Issue 
A&A
Volume 619, November 2018



Article Number  A117  
Number of page(s)  8  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201732242  
Published online  15 November 2018 
A dedicated sourceposition transformation package: pySPT^{⋆}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
email: owertz@alumni.ulg.ac.be
Received:
6
November
2017
Accepted:
23
January
2018
Modern timedelay cosmography aims to infer the cosmological parameters with a competitive precision from observing a multiply imaged quasar. The success of this technique relies upon a robust modeling of the lens mass distribution. Unfortunately strong degeneracies between density profiles that lead to almost the same lensing observables may bias precise estimates of the Hubble constant. The source position transformation (SPT), which covers the wellknown masssheet transformation (MST) as a special case, defines a new framework to investigate these degeneracies. In this paper, we present pySPT, a python package dedicated to the SPT. We describe how it can be used to evaluate the impact of the SPT on lensing observables. We review most of its capabilities and elaborate on key features that we used in a companion paper regarding SPT and time delays. The pySPT program also comes with a subpackage dedicated to simple lens modeling. This can be used to generate lensing related quantities for a wide variety of lens models independent of any SPT analysis. As a first practical application, we present a correction to the first estimate of the impact on time delays of the SPT, which has been experimentally found in a previous work between a softened power law and composite (baryons + dark matter) lenses. We find that the large deviations previously predicted have been overestimated because of a minor bug in the public lens modeling code lensmodel (v1.99), which is now fixed. We conclude that the predictions for the Hubble constant deviate by ∼7%, first and foremost as a consequence of an MST. The latest version of pySPT is available on Github, a software development platform, along with some tutorials to describe in detail how making the best use of pySPT.
Key words: cosmological parameters / gravitational lensing: strong
© ESO 2018
1. Introduction
For about a decade, the modern timedelay cosmography, namely the cosmological parameter inferences from time delay measurements in strong gravitational lensing, have been achieved with an increasingly competitive precision (for a recent review, Treu & Marshall 2016; see and references therein). A crucial step of this technique relies upon a robust characterization of the gravitational potential, which produces the multiply imaged configuration of a background bright quasar (see, e.g., Keeton 2003; Fassnacht et al. 2006; Suyu et al. 2010, 2013; Wong et al. 2011; Wong et al. 2017; Birrer et al. 2016). This gravitational potential is essentially produced by a main deflector but also by any mass distributions lying along the line of sight to the source (Seljak 1994; BarKana 1996). Unfortunately, modeling the main lens mass distribution faces a major hurdle, namely the existence of degeneracies between plausible lens density profiles. In fact, significant freedom exists in choosing lens models that produce the same image configurations but predict different products of time delays and Hubble constant, H_{0} Δt (Schneider & Sluse 2013). Thus, these degeneracies translate into systematic errors that propagate to H_{0}.
Now wellknown to the lensing community, the first lensing invariance to have been pointed out is the masssheet degeneracy (MSD; Falco et al. 1985). A dimensionless surface mass density κ(θ) and all the modified κ_{λ}(θ) under the masssheet transformation (MST) defined as
along with the corresponding unobservable source rescaling β → λ β, lead to identical lensing observables with the exception of time delays, which transform as in Δt → λ Δt. This pure mathematical degeneracy has nothing to do with the gravitational perturbations caused by external masses along the line of sight. Whereas various solutions have been already proposed to reduce the impact of the MST on timedelay cosmography (see, e.g., Sect. 3 in Treu & Marshall 2016; and references therein), none of them succeed in unambiguously breaking the MSD. The source position transformation (SPT), a more general class of degeneracies which includes the MST as a special case, has been introduced in Schneider & Sluse (2014) and carried forward in Unruh et al. (2017). The SPT defines a new mathematical framework that includes degeneracies that have been neither described nor considered in timedelay cosmography before. The first estimation of its impact on time delays is given in Schneider & Sluse (2014) where the authors show experimentally that predictions for H_{0} can deviate by ∼20% when comparing a softened power law (SPL) with a composite (baryons + dark matter) lens models. However, this value has been overestimated owing to a bug in the software used by the author. We address this issue in Sect. 5 and show that the predictions for H_{0} can deviate by ∼7% instead.
Recently, a detailed analysis of how the SPT may affect the timedelay cosmography has been presented in Wertz et al. (2018). To address this question, we started by developing a flexible numerical framework that encompasses welltested and efficient implementations of most of the analytical results published in Schneider & Sluse (2014) and Unruh et al. (2017). Numerous additional features were then added, giving rise to pySPT, an easytouse and welldocumented python package dedicated to the SPT. We used pySPT to draw the conclusions presented in Wertz et al. (2018). Thus, our package is released as open source, making our results easy to reproduce. Beyond that, we also hope that it will be useful to the timedelay cosmography community to quantify the systematic errors that are introduced by the SPT when inferring H_{0}.
This paper is organized as follows. We outline the basic principles of SPT in Sect. 2. Section 3 gives an overall description of the package design and features, and details are provided in Sect. 4. In Sect. 5, we present the corrected version of some results presented in Schneider & Sluse (2014). We summarize our findings in Sect. 6.
2. Principle of source position transformation
This section focuses on the basic idea underlying the SPT and its mathematical framework. A detailed discussion is provided in Schneider & Sluse (2014) and Unruh et al. (2017). In this paper we adopt the standard gravitational lensing notation (see Schneider 2006).
The relative lensed image positions θ_{i}(θ_{1}) of a background pointlike source located at the unobservable position β constitute the lensing observables that we measure with the highest accuracy and precision. As a typical example, just a few mas can be achieved with deep HST observations. When n images are observed, the mapping θ_{i}(θ_{1}) only provides the constraints
where α(θ) corresponds to the deflection law caused by a foreground surface mass density κ(θ), the socalled lens. The SPT addresses the question of whether we define an alternative deflection law, denoted as , which preserves the mapping θ_{i}(θ_{1}) for a unique source? If such a deflection law exists, the alternative source position differs in general from β. Furthermore, it defines the new lens mapping , which leads to
An SPT consists in a global transformation of the source plane formally defined by a onetoone mapping , unrelated to any physical contribution such as the external convergence (Schneider & Sluse 2013). To preserve the mapping θ_{i}(θ_{1}), the alternative deflection law is thus written
where in the first step we used Eq. (3) and in the last step we inserted the original lens equation. As defined, the deflection laws α(θ) and yield exactly the same image positions of the source β and , respectively.
Because is in general not a curlfree field, it cannot be expressed as the gradient of a deflection potential caused by a mass distribution κ̂. Provided its curl component is sufficiently small, Unruh et al. (2017) have established that one can find a curlfree deflection law ᾶ = ∇ψ̃ that is similar to in a finite region. In other words, ᾶ yields the same sets of multiple images up to the astrometric accuracy ε_{acc} of current observations. Two image configurations are considered indistinguishable when they satisfy
for all images θ of the source β, and corresponding images θ̃ of the source with . In Unruh et al. (2017), the adopted similarity criterion is written
in a finite region of the lens plane denoted as 𝒰 where multiple images occur. Even though this criterion cannot actually guarantee Δθ< ε_{acc} over 𝒰, there exists in general a subregion 𝒰′⊂𝒰 that includes image configurations for which Eq. (5) is satisfied (see Sect. 4.1 in Wertz et al. 2018). Thus, an SPT leads to an alternative lens profile κ̂, which gives rise to the SPTtransformed deflection law , then its curlfree counterpart ᾶ is defined based upon the criterion Eq. (6), and may lead to indistinguishable image configurations produced by α. To conclude this section, we note that the deflection law ᾶ is produced by a surface mass density κ̃ := ∇ · ᾶ/2, which equals κ̂ by construction.
3. Package overview
The pySPT program was developed in python (Watters et al. 1996) and only relies on packages included in the python standard library^{1} and the proven opensource libraries numpy (Van Der Walt et al. 2011), scipy (Jones et al. 2001), and matplotlib (Hunter 2007). The code design and development follow effective practices for scientific computing such as using test cases and a profiler to identify bugs and bottlenecks, keeping an effective collaboration thanks to a version control system (git), and providing an extensive documentation (Wilson et al. 2014). This open source code is available on GitHub^{2} and comes along with a clear description on how to install it. To make the best use out of pySPT, a quick start guide and several tutorials are also provided in the form of Jupyter notebooks. Benefiting from the python objectoriented paradigm, the structure of pySPT is highly modular, which avoids code clones and makes it less sensitive to bug propagation.
The code is composed of several modules built from various classes and is organized in a dozen subpackages. Its core is separated into three main subpackages, referred to as lensing, sourcemapping, and spt The reason for that is simple. The mapping presented in Sect. 2 defines a onetoone connection between the source plane and its SPTtransformed counterpart. Through the lens equation β = θ − α(θ), the deflection law (arising from a lens model) characterizes the link between the source plane and the image plane. Thus, only together with α(θ) does give rise to the alternative deflection law defined in Eq. (4). As a result, the most significant pySPT subpackages are lensing to basically deal with α(θ), sourcemapping to define , and spt to describe (and all the SPTtransformed lensing quantities). In the remainder of this section, we describe briefly their main functionalities^{3}.
The subpackage lensing shares a lot of functionalities with other lensingdedicated softwares such as gravlens (Keeton 2001b). From its class Model, we generated a lens model object that allowed us to perform a wide range of basic lensing calculations. Strictly speaking, this part of the code is not related to the SPT and it can be used independent of any SPT analysis. Nevertheless, we decided to implement this part of the code for a simple reason: as a builtin feature of pySPT, the class Model provides to its instances the adequate structure that is required to match the spt requirements. This makes lensing more convenient to use than a third party code. The subpackage sourcemapping is used to define an SPT . Several forms of are implemented, such as particular radial stretchings presented in Schneider & Sluse (2014). Moreover, the code is designed to accept any userdefined SPT. With sourcemapping also comes functionalities to characterize the mapping , such as testing whether it is onetoone over a given region. More details are given in Sect. 4.2. The subpackage spt is the heart of pySPT. It is designed to work alongside with sourcemapping and lensing to provide all the basic SPTtransformed quantities, such as . One can also derive the SPTtransformed time delays Δt̂ and Δt̃, and most of the quantities presented in Schneider & Sluse (2014), Unruh et al. (2017), and Wertz et al. (2018). This makes the results presented in these papers straightforward to reproduce. A detailed description of the tools provided by spt is given in Sect. 4.3.
The pySPT code also includes several subpackages dedicated to specific tasks. Based on the package multiprocessing of the standard library, multiproc provides an efficient tool to parallelize functions and methods in a straightforward way. As such, most of the pySPT features support parallel computing to leverage multiple processors fully. The grid feature helps us to create different types of mesh grids. These are of practical interest for generating maps of lensing quantities in a particular region. To calculate ᾶ efficiently, we followed Unruh et al. (2017), which suggests using a Riemann mapping to handle a pole numerically. Thus, pySPT contains the subpackage integrals that includes functionalities to deal with and to illustrate conformal mappings.
4. Analyzing the impact of the SPT with pySPT
4.1. Subpackage lensing to deal with lens models
To analyze how the SPT may alter lensing observables in a quantitative way, we needed to select a model that characterizes a mass distribution. A range of lensing observables were then generated and compared against those produced by an SPTtransformed version of the original model. The subpackage lensing relies on the class Model of which one of its instances defines a lens model and provides efficient tools to compute a wide range of lensing quantities. Most of the mass profiles described in Keeton (2001a) are available in the subpackage catalog, together with a complete documentation. We note that, to our knowledge, no analytical expression of the deflection potential for the generalized pseudoNFW can be found in the literature. Hence, we derived this expression and provide the result in Appendix A. The class Model is also designed to accept a userdefined lens model, as long as some conventions are respected. At least, the deflection potential ψ(θ) must be provided as either a python function or a C shared library. Those that are not defined among α, κ, and ∂α/∂θ are then computed from numerical approximation at the expense of a more intensive usage of computer resources. Thanks to operator overloading^{4}, arbitrarily complicated composite models can be obtained by simply adding several Model instances. For example, this feature is used to generate quadrupole models, namely the combination of an axisymmetric matter distribution plus some external shear. To go a step further, the axisymmetric part may itself be composed of different components, whereas additional contributions can be included at arbitrary positions.
The computational methods implemented with the class Model follow the prescription of Keeton (2001b). In particular, the socalled tiling algorithm (with adaptive grid) is used to solve the lens equation and locate the critical curves. Other fundamental lensing quantities are available, such as the caustics, basic image properties (e.g., image type, magnification factor, parity, oddnumber, and magnification theorem checks), Fermat potential τ(θ), and time delays between image pairs Δt(θ_{i}, θ_{j})≡Δt_{ij}. In Fig. 1, we illustrate the capabilities of the subpackage lensing by showing several lensing quantities produced by a complex mass distribution. The workflow for generating data used in this figure is as follows^{5}:
Fig. 1. Example illustrating some capabilities of the subpackage lensing. Left panel: eight lens profiles with different parameters are combined to generate a complex mass distribution. The total surface mass density is shown in tones of gray and dashed curves highlight few particular isodensity contours. The colored thick curves show the critical curves and the filled and empty markers locate the lensed image positions of two different sources shown in the right panel. The inverted triangles correspond to images of type I (maximum of τ), the diamonds to type II (saddle point of τ), and triangles to type III (maximum of τ). The size of the markers is logproportional to the magnification of the images. Right panel: the colored lines show the caustics and the two squares locate the sources. The filled (resp. empty) square has seven (resp. nine) images, all shown in the left panel. The axis scale is arcseconds but unit of θ_{E} can also be used. 

Open with DEXTER 
The subpackage lensing also includes C shared libraries that implement ψ, α, κ, ∂α/∂θ for all the lens models. As shown in Sect. 4.3, the curlfree deflection angle ᾶ is obtained by means of line and surface integrals of functions that involve (see Eq. (4)) and (Unruh et al. 2017). Thus, the original deflection angle α is evaluated as many times as the number of iterations required by the solver to converge. Even though this procedure is efficient when only one ᾶ is evaluated, high performance becomes critical when ᾶ needs to be evaluated on a dense grid. Thanks to the foreign function interface module ctypes^{6}, the use of the C shared libraries speeds up significantly the execution of pure python code. The pySPT code comes along with a tutorial in the form of a Jupyter notebook that describes in details how to deal with C shared libraries.
4.2. Defining with the subpackage sourcemapping
After we defined the lens model, we chose an SPT by defining a source mapping . To each position β of the source plane, this mapping associates a new and unique position in the source plane. The socalled radial stretching is simply defined as
where f is called the deformation function. With defined this way, always lies on the line passing through 0 and β_{j}. The most simple case of radial stretching corresponds to a constant deformation function, f(β) = λ − 1 with λ ∈ ℝ, which leads to the wellknown MST, . In Wertz et al. (2018), we focused most of the work on the radial stretching Eq. (7) where the deformation function f(β) is defined as the lowest order expansion of more general functions
where is the Einstein angular radius, and f is an even function of β to preserve the symmetry (Schneider & Sluse 2014). When f_{2} = 0, this case simplifies to a pure MST with λ = 1 + f_{0}. In Table 1, we provide a list of the radial stretchings implemented in sourcemapping. The rationale behind the choice of these particular deformation functions are motivated in Schneider & Sluse (2014). The first column refers to the id used to identify which source mapping one wants to select. A specific example is given below. Aside from defining an SPT, sourcemapping also includes a few functionalities for characterizing the mapping . In particular, one can test whether the source mapping is one to one over a given region, compute the Jacobi matrix , and provide its decomposition ℬ(β)=B_{1}ℐ + B_{2}Γ(η), which is useful to evaluate the amplitude of the curl component of (see Sect. 4.1 in Schneider & Sluse 2014). As an example, the code below illustrates how to define and characterize an MST with λ = 1(≡1 + f_{0}): Similarly, we can first define the deformation function f that characterizes an MST with λ as unique argument, and pass it to sourcemapping as an argument. For efficiency, the first derivative f′ of the deformation function with respect to β can also be passed^{7}. This process illustrates how we can work with a userdefined SPT.
List of radial stretchings implemented in lensmapping.
4.3. Deriving SPTtransformed quantities with spt
In the two previous sections, we have shown how to generate lensing observables produced by a given lens model and how to define an SPT. We present the subpackage spt, which provides the tools required to determine the SPTtransformed quantities .
For any SPT and lens model, the alternative deflection law is implemented as defined in Eq. (4), and the Jacobi matrix of the alternative lens equation given in Eq. (3) as . The determination of perfectly illustrates how spt works along with lensing and sourcemapping: ℬ(β) is obtained with sourcemapping, β(θ) and 𝒜(θ) are obtained with lensing, and spt combines the results to derive . In the axisymmetric case, is a curlfree field and there exists a deflection potential such that where θ = θ. Thus, is obtained as follows:
up to a constant independent of θ. Because the integrand in Eq. (9) depends on the SPT and the lens model choices, no general analytical solution can be derived. The integral is therefore computed numerically using the python module integrate.quad^{8} from scipy (Jones et al. 2001). Under a radial stretching, the axisymmetric mass profile κ(θ) transforms into , where β = β (Schneider & Sluse 2014). Otherwise, the more general form is valid regardless of the lens model symmetry.
When the axisymmetry assumption for the original lens model is dropped, spt also provides the physically meaningful ᾶ and ψ̃. The analytical expressions implemented in SPT are slightly simplified versions of those first presented in Unruh et al. (2017). The deflection potential ψ̃ evaluated at the position θ in the lens plane is explicitly written
where the region 𝒰 is a disk of radius R, 〈ψ̃〉 is the average of ψ̃ on 𝒰, ds the line element of the boundary curve ∂𝒰,
and
In Appendix B, we explicitly show that both versions are fully equivalent. The corresponding simplified version of the deflection angle ᾶ can be derived by obtaining the gradient of H_{1} and H_{2} with respect to θ, which is written
To deal with the pole ϑ = θ in the first term of Eqs. (11) and (13), we used a Riemann mapping as described in Appendix A in Unruh et al. (2017) and implemented in the subpackage integrals. This mathematical trick makes the previous integrals easier to solve by mapping 𝒰 onto the unit disk in the complex plane, such as the pole is moved at the origin^{9}z = 0. Additional care is however needed in the vicinity of θ = 0, where the gradient of κ(θ) may vary significantly and produce a second sharp peak of the integrand. This peak may even be a new pole when κ is singular at the origin. Thus, a physically meaningful lens model should always be favored. In the left panel in Fig. 2, we illustrate the integrand^{10} of the first term of the integral over 𝒰 defined in Eq. (10). For this illustrative example, we choose an nonsingular isothermal sphere (NIS, θ_{c} = 0.1θ_{E}) lens model plus external shear (γ_{p} = 0.1) transformed by the radial stretching with (f_{0}, f_{2})=(0, 0.55) and θ = ( − 0.5, 0.25) θ_{E}. The negative peak comes from the pole ϑ = θ while the central peak is caused by κ(ϑ=0). The right panel in Fig. 2 illustrates the integrand after applying the Riemann mapping. The peaks have moved and the pole ϑ = θ now lies at the origin z = 0, as expected. To deal with the second peak, the subpackage integrals splits the integration domain to place it on a boundary and to ensure the gradient of the integrand is smooth in the sub\domains.
Fig. 2. Representative example that illustrates how the integration domain 𝒰 is mapped onto the unit disk in the complex plane under a Riemann mapping. Left panel: the color coding depicts the integrand for all ϑ ∈ 𝒰. It shows a pole at ϑ = θ with θ = ( − 0.5, 0.25) θ_{E} and a secondary peak at the origin caused by κ(ϑ=0). The lens model is an NIS (θ_{c} = 0.1θ_{E}) plus external shear (γ_{p} = 0.1) transformed by the radial stretching with (f_{0}, f_{2})=(0, 0.55). The red circle delimits 𝒰 with radius R. Right panel: integrand after applying the Riemann mapping described in Appendix A in Unruh et al. (2017). The pole ϑ = θ now lies at the origin of the unit (blue) circle of the complex plane. The polar grid (gray lines) allows us to visualize how the Riemann mapping acts on 𝒰. For this figure, we used the subpackage integrals that addresses both the pole and second peak. 

Open with DEXTER 
With and ᾶ, we can evaluate the SPT validity criterion adopted in Unruh et al. (2017) and recalled in Eq. (6). With the same lens model and SPT adopted for Figs. 2 and 3 shows the map Δα(θ) over a circular grid θ ≤2 θ_{E} in the lens plane. It is worth noting that this figure is similar to the map Δα(θ) illustrated in the Fig. 7 in Unruh et al. (2017) while ᾶ was obtained from two different and independent approaches. Unruh et al. (2017) first calculated ψ̃ by solving numerically a Neumann problem thanks to a successive overrelaxation method (Press et al. 1992) on a square grid of width 4 θ_{E}; then they derived ᾶ from ψ̃ using a secondorder accurate finite differencing scheme (see their Sect. 3.2 for a detailed overview). This iterative process necessarily requires systematically calculating ψ̃ over the whole square grid. Conversely, in Fig. 3, ᾶ is obtained directly from the explicit Eq. (13) for each position on a circular sampling grid^{11}. Thus, although the similarity between the two figures confirms the consistency of the two approaches, the semianalytical approach implemented in pySPT yields ᾶ at a particular position.
Fig. 3. Map of Δ α(θ) over a circular grid θ≤2 θ_{E} for f_{2} = 0.55, θ_{c} = 0.1 θ_{E} and γ_{p} = 0.1. We set the radius R of the circular region 𝒰 in such a way that the area of 𝒰 is equal to the area of the square grid used in the pure numerical approach, i.e., . This figure is similar to the Fig. 7 in Unruh et al. (2017) even though it is based on a different approach (see the text for more details). This figure has been obtained with the subpackage spt in less than five minutes for a grid of about 2 × 10^{4} positions. 

Open with DEXTER 
In Wertz et al. (2018), we analyzed the impact of the SPT on time delays in detail. To achieve this, we compared, for a given lens model, the time delays Δt_{ij} between image pairs (θ_{i}, θ_{j}) of a source β with the time delays Δt̃_{ij} between the image pairs (θ̃_{i}, θ̃_{j}) of the modified source under an SPT. The images θ̃ satisfy the lens equation . By construction of ᾶ, we expect θ̃_{i} to be close to θ, at least in a subregion 𝒰′⊂𝒰 (see Sect. 4 in Wertz et al. 2018). These images θ̃ can be obtained easily thanks to the subpackage lensing. Because its main class Model is designed to accept a userdefined lens model, we benefited from all the tools provided by lensing to characterize the SPTtransformed lens model (see Sect. 4.1). The code below illustrates how to generate the SPTtransformed lensing quantities thanks to the use of the subpackages lensing, sourcemapping, and spt.
5. Impact of the SPT on time delays: Empirical estimation
The first empirical estimation of the impact of the SPT on time delays was presented in Schneider & Sluse (2013, 2014). These authors showed that a quadrupole model composed of a SPL profile predicts the same lensed image positions (with a 0.004 arcsec accuracy) as a composite fiducial model (Hernquist profile + generalized NFW + external shear) for a set of sources and β, respectively^{12}. We want to stress that the SPL model is not an SPTgenerated model from the composite fiducial model but they represent two different models for which the nature of the degeneracy can be approximated by an SPT. Thus, the set of sources was obtained independent of β by fitting the lensed image positions produced by the fiducial model. The top panel of the Fig. 4 in Schneider & Sluse (2014) represents as a function of β. For sources located in a disk of radius 0.7 arcsec, the connection between and β is slightly anisotropic and roughly resembles a radial stretching of the form Eq. (8) with f_{0} = −0.068 and f_{2} ≈ 0.012. The bottom panel of the Fig. 4 in Schneider & Sluse (2014) shows that were almost constant in the double image regime and not conserved in the quadruple image regime. These authors also found that was never smaller than 1.2, reaching a mean of 1.45 for the quadruple image configuration of a source located at β ≈0.18 arcsec. Thus, this figure shows that the degeneracy between the SPT and fiducial models can affect the inferred value of H_{0} by an average of 20% to up to 45% for particular image configurations.
Fig. 4. Time delay ratios of image pairs between the composite fiducial model and its SPTtransformed counterpart under a radial stretching with (f_{0}, f_{2})=(−0.068, 0.012). Top panel: Δt̂/Δt normalized by the ratio are very close to 1, even in the quadruple image regime, which disagree with the results obtained in Schneider & Sluse (2014). Bottom panel: The impact of the SPT (cleaned from the pure MST with λ = 1 + f_{0} = 0.932) is around only a few tenths of percent. The error bars illustrates that the time delay ratios are not conserved in the quadruple image regime, as first suggested in Schneider & Sluse (2014). 

Open with DEXTER 
As a first application of pySPT, we compare these time delay ratios with those obtained when we transform the fiducial model under a radial stretching Eq. (8) with (f_{0}, f_{2})=(−0.068, 0.012). As shown in Fig. 4, we find that the impact of the SPT on time delays is much smaller than predicted in Fig. 4 in Schneider & Sluse (2014). For instance, they found that a source located at β=0.7 arcsec leads to Δt̂/Δt ≈ 1.169 (between the two outer images), whereas we find Δt̂/Δt ≈ 0.935, knowing that the major contribution comes from an MST with λ = 1 + f_{0} = 0.932. Moreover, the small anisotropic feature of the empirical source mapping alone cannot explain this tension.
The discrepancy between our prediction and the results obtained in Schneider & Sluse (2014) are triggered by a minor bug in the public lens modeling code lensmodel (v1.99; Keeton 2001a) that the authors used to compute the time delays. We spotted this bug when we compared the outputs produced with our package pySPT and lensmodel for the deflection angle and deflection potential for the SPL model. Denoting the logarithm slope of the SPL model as a, we found that and , showing adifferent normalization factor between . This extra normalization factor propagates into the code, leading to a biased value of the time delay Δt̂_{lensmodel}. It is worth mentioning that, for isothermal profiles (a = 1), this extra normalization factor a has no impact on the lensing quantities computed with lensmodel. This might explain why this minor bug has remained unnoticed so far. Nevertheless, the latter has been fixed and a corrected version of lensmodel has been immediately released by Chuck Keeton. To obtain the Fig. 4 in Schneider & Sluse (2014), the SPL model was characterized by the logarithmic threedimensional slope γ′=2.24, which is linked to a by the relation a = 3 − γ′, hence a = 0.76 ≠ 1, which finally explains the discrepancy mentioned before.
In Fig. 5, we show the corrected version of the Fig. 4 that we produced with our package pySPT. We confirm that the exact same result can now be obtained from the corrected version of lensmodel. The normalized time delay ratios plotted against β (top panel) are very close to 1. Thus, the time delay ratios (Δt̂/Δt) (bottom panel) closely resembles the source ratios represented in top panel in Fig. 4. This behavior is well understood and is described in detail in the companion paper Wertz et al. (2018). The impact of the SPT (separated from the MST with 1 + f_{0} = 0.932) on time delays now reaches only around 0.32% for β=0.7 arcsec, which is obviously much smaller than the 20% previously found in Schneider & Sluse (2014). The corrected version now fully agrees with the conclusions drawn from Fig. 4. In addition, this confirms that the degeneracy between the SPT and the fiducial models mimics an SPT, as first established in Unruh et al. (2017).
Fig. 5. Time delay ratios of image pairs between the composite fiducial model and quadrupole composed of a SPL. This figure constitutes the corrected version of the Fig. 4 published in Sect. 4.3 in SS14. The hat lensing quantities correspond to the SPL model while standard notation is used for the composite fiducial model. 

Open with DEXTER 
6. Conclusions
In this paper we have presented the pySPT package for analyzing the SPT. The pySPT program relies on several subpackages, the most important of which are lensing to deal with lens model, sourcemapping to define an SPT, and SPT to provide the SPTtransformed lensing quantities defined in previous papers. The subpackage lensing is particular in a sense that it can be used independent of any SPT analysis. To some extent, it somewhat offers a python alternative to the public lens modeling code lensmodel with which it shares a lot of functionalities. The pySPT program implements functionalities for generating lensing quantities produced by SPTtransformed lens models. Thanks to its modularity, pySPT is also designed to accept both userdefined lens model and SPT. In such a case, pySPT constitutes a user friendly interface to deal with the SPT.
As a first application, we used pySPT to explore how a radial stretching may affect the time delay measurements for a fiducial model composed of a Hernquist profile + generalized NFW + external shear. We found that the impact of the SPT is much smaller than first suggested in Schneider & Sluse (2014). We addressed the tension between these results by spotting a minor bug in the public lens modeling code lensmodel that was used by the authors. This bug resulted in biased values of the deflection potential for the nonisothermal SPL model, leading to an overestimated impact of the SPT on the time delays. Using the subpackage lensing, we produced a corrected version of the Fig. 4 published in Schneider & Sluse (2014), which now fully agrees with the results presented in this paper. As a result, the impact of the SPT on timedelay cosmography might not be as crucial as initially suspected. We address this question in details in the companion paper Wertz et al. (2018).
With the next version of pySPT, we plan to include stateoftheart lens modeling capabilities. For example, combining stellar dynamics data obtained from spectroscopy of the lens galaxy with lensing measurements has become a standard practice within the strong lensing community. Furthermore, we still do not have a clear answer to the question of how the kinematic information of a mass distribution is affected under an SPT. Schneider & Sluse (2013) showed that the fiducial and SPL models discussed in Sect. 5 could not be satisfactorily distinguished thanks to the measurement of the stellar velocity dispersion σP with a typical 10% uncertainty. This thus suggests that the use of σP may be of limited help for breaking the SPT. In a future work, we aim to address this open question with the use of pySPT. In this context, we plan to update the software with a new subpackage fully dedicated to the determination of the stellar velocity dispersion associated with an SPTmodified mass profile. Finally, pySPT is provided not only as an ultimate tool for the lens modeling community, but primarily as an attractive choice to identify the possible degeneracies from which timedelay cosmography may suffer. Nonetheless we hope that its permanent development will attract more users and will extend its purpose to more than just dealing with the SPT.
https://github.com/owertz/pySPT. For git users, pySPT can be cloned directly from the source code repository via the following bash command git clone recursive https://github.com/owertz/pySPT
We adopt the naming conventions advocated in the Python Enhancement Proposals 8 (PEP8). In particular, subpackages have alllowercase names while class names use the socalled CapWords convention. The PEP8 is accessible at https://www.python.org/dev/peps/pep0008/
A Jupyter notebook dedicated to this figure is available in the pySPT repository on GitHub. In addition to details about the workflow shown above, it includes the code we used to plot Fig. 1.
ctypes is included in the python standard library: https://docs.python.org/2.7/library/ctypes.html
For axisymmetric profile κ(θ)=κ(θ), the analytical expression for involves f′ (see Eq. (16) in Schneider & Sluse 2014).
This package provides an interface to QUADPACK (Piessens et al. 1983) whose routines use the adaptive quadrature method to approximate integrals.
We may have chosen a square or another shape for the sampling grid. The choice of a disk is motivated by the fact that the region in which multiple images occur is typically θ≤2 θ_{E}. Furthermore, the radius of the circular sampling grid over which we evaluate Δ α, which is depicted in Fig. 3, is not the radius R of 𝒰. For each position on the grid, solving the Eq. (13) requires to define R. For consistency, we must adopt the same R for each evaluation of ᾶ, which implies that R must be chosen (at least equal or) larger than the radius of the circular sampling grid.
We follow the same notation as in Schneider & Sluse (2014), i.e., we denote the lensing quantities associated with the SPL model with a hat, e.g., Δt̂ for the time delay, whereas no hat is used for the composite fiducial model.
Acknowledgments
We would like to thank Dominique Sluse and Chuck Keeton for valuable discussions that allowed us to spot and fix the small bug in lensmodel. We are very grateful to the anonymous referee for his comments and suggestions that contributed to improving the quality of pySPT and this paper. This work was supported by the Humboldt Research Fellowship for Postdoctoral Researchers.
References
 BarKana, R. 1996, ApJ, 468, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmol. Astropart. Phys., 8, 020 [NASA ADS] [CrossRef] [Google Scholar]
 Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fassnacht, C. D., Gal, R. R., Lubin, L. M., et al. 2006, ApJ, 642, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, Online; accessed 20170807. [Google Scholar]
 Keeton, C. R. 2001a, ArXiv eprints [arXiv:astroph/0102341]. [Google Scholar]
 Keeton, C. R. 2001b, ArXiv eprints [arXiv:astroph/0102340] [Google Scholar]
 Keeton, C. R. 2003, ApJ, 584, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, J. A., Kochanek, C. S., & Keeton, C. R. 2001, ApJ, 558, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Piessens, R., de DonckerKapenga, E., Überhuber, C., & Kahaner, D. 1983, QUADPACK: A Subroutine Package for Automatic Integration (Berlin, Heidelberg: SpringerVerlag), 0179 [Google Scholar]
 Press, W. H., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C 2nd edn.: The Art of Scientific Computing (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Schneider, P. 2006, in SaasFee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro, eds. G. Meylan, P. Jetzer, P. North, C. S. Kochanek, & J. Wambsganss, 1 [Google Scholar]
 Schneider, P., & Sluse, D. 2013, A&A, 559, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Sluse, D. 2014, A&A, 564, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seljak, U. 1994, ApJ, 436, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., & Marshall, P. J. 2016, A&ARv, 24, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Unruh, S., Schneider, P., & Sluse, D. 2017, A&A, 601, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [CrossRef] [Google Scholar]
 Watters, A., Ahlstrom, J. C., & Rossum, G. V. 1996, Internet Programming with Python (New York: Henry Holt and Co., Inc.) [Google Scholar]
 Wertz, O., Orthen, B., & Schneider, P. 2018, A&A, 617, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wilson, G., Aruliah, D. A., Titus Brown, C., et al. 2014, PLOS Biol., 12, e1001745 [CrossRef] [Google Scholar]
 Wong, K. C., Keeton, C. R., Williams, K. A., Momcheva, I. G., & Zabludoff, A. I. 2011, ApJ, 726, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Generalized pseudoNFW
To our knowledge, no analytical expression of the deflection potential ψ(θ) for the generalized pseudoNFW model has ever been published in the literature. For practical purposes, we present such an analytical expression, which has been implemented into pySPT. We first recall that the spherical density distribution ρ(r) of the generalized pseudoNFW model is defined as (see the Eq. (1) in Muñoz et al. 2001, with n = 3)
where ρs is a characteristic density, rs the scale radius, and γ the logarithmic slope of the density profile at small radius. Up to an additive constant, the axisymmetric deflection potential ψ(θ) is given by
where κ_{s} = ρ_{s} r_{s}/Σ_{cr}, the term Γ^{*}(υ)=Γ(υ)/Γ(υ − 1/2) is a particular combination of the gamma function Γ,
where is a particular Gauss hypergeometric function, and the dilogarithm Li_{2}(z) can be defined by the series
Appendix B: Proof of the relation (10)
As a preamble, the notation adopted here differs from that used in Unruh et al. (2017). We use θ as a position in the lens plane and ϑ as the corresponding integration variable for θ. The inverse is partially used in Unruh et al. (2017), in particular in the Sect. 3.3 where ψ̃ and ᾶ are derived.
Starting with Eq. (18) in Unruh et al. (2017), the deflection potential ψ̃ evaluated at the position θ is given by
where a solution for the Green function H is analytically known when 𝒰 is a disk of radius R
Firstly, we note that ϑ=R for all ϑ on the boundary ∂𝒰, which implies
Thus, the two logarithmterms in H(θ; ϑ) are equal when we consider the line integral.
Secondly, the term −θ^{2}/(4πR ^{2}) in H(θ; ϑ) does not depend on ϑ, and therefore contributes neither to the integral over 𝒰 nor to the line integral. Therefore, Eq. (B.1) contains the term
where the equality holds because of and we made use of Gauß divergence theorem. As a result, the term −θ^{2}/(4πR ^{2}) in H(θ; ϑ) does not contribute to ψ̃.
Finally, combining Eqs. (B.1), (B.3), and (B.4) leads to the definition of ψ̃ given in Eq. (10). We note that the same reasoning holds for ᾶ, leading to the Eq. (13).
All Tables
All Figures
Fig. 1. Example illustrating some capabilities of the subpackage lensing. Left panel: eight lens profiles with different parameters are combined to generate a complex mass distribution. The total surface mass density is shown in tones of gray and dashed curves highlight few particular isodensity contours. The colored thick curves show the critical curves and the filled and empty markers locate the lensed image positions of two different sources shown in the right panel. The inverted triangles correspond to images of type I (maximum of τ), the diamonds to type II (saddle point of τ), and triangles to type III (maximum of τ). The size of the markers is logproportional to the magnification of the images. Right panel: the colored lines show the caustics and the two squares locate the sources. The filled (resp. empty) square has seven (resp. nine) images, all shown in the left panel. The axis scale is arcseconds but unit of θ_{E} can also be used. 

Open with DEXTER  
In the text 
Fig. 2. Representative example that illustrates how the integration domain 𝒰 is mapped onto the unit disk in the complex plane under a Riemann mapping. Left panel: the color coding depicts the integrand for all ϑ ∈ 𝒰. It shows a pole at ϑ = θ with θ = ( − 0.5, 0.25) θ_{E} and a secondary peak at the origin caused by κ(ϑ=0). The lens model is an NIS (θ_{c} = 0.1θ_{E}) plus external shear (γ_{p} = 0.1) transformed by the radial stretching with (f_{0}, f_{2})=(0, 0.55). The red circle delimits 𝒰 with radius R. Right panel: integrand after applying the Riemann mapping described in Appendix A in Unruh et al. (2017). The pole ϑ = θ now lies at the origin of the unit (blue) circle of the complex plane. The polar grid (gray lines) allows us to visualize how the Riemann mapping acts on 𝒰. For this figure, we used the subpackage integrals that addresses both the pole and second peak. 

Open with DEXTER  
In the text 
Fig. 3. Map of Δ α(θ) over a circular grid θ≤2 θ_{E} for f_{2} = 0.55, θ_{c} = 0.1 θ_{E} and γ_{p} = 0.1. We set the radius R of the circular region 𝒰 in such a way that the area of 𝒰 is equal to the area of the square grid used in the pure numerical approach, i.e., . This figure is similar to the Fig. 7 in Unruh et al. (2017) even though it is based on a different approach (see the text for more details). This figure has been obtained with the subpackage spt in less than five minutes for a grid of about 2 × 10^{4} positions. 

Open with DEXTER  
In the text 
Fig. 4. Time delay ratios of image pairs between the composite fiducial model and its SPTtransformed counterpart under a radial stretching with (f_{0}, f_{2})=(−0.068, 0.012). Top panel: Δt̂/Δt normalized by the ratio are very close to 1, even in the quadruple image regime, which disagree with the results obtained in Schneider & Sluse (2014). Bottom panel: The impact of the SPT (cleaned from the pure MST with λ = 1 + f_{0} = 0.932) is around only a few tenths of percent. The error bars illustrates that the time delay ratios are not conserved in the quadruple image regime, as first suggested in Schneider & Sluse (2014). 

Open with DEXTER  
In the text 
Fig. 5. Time delay ratios of image pairs between the composite fiducial model and quadrupole composed of a SPL. This figure constitutes the corrected version of the Fig. 4 published in Sect. 4.3 in SS14. The hat lensing quantities correspond to the SPL model while standard notation is used for the composite fiducial model. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.