Unified lensing and kinematic analysis for any elliptical mass profile

TL;DR

This paper introduces a fast, general method for calculating gravitational lensing effects for any elliptical mass profile by decomposing it into Gaussian components, solving a long-standing computational challenge.

Contribution

It presents a novel, efficient technique to compute lensing quantities for elliptical profiles using Gaussian decomposition, overcoming a three-decade-old computational hurdle.

Findings

Method accurately computes lensing deflection and magnification for elliptical profiles.

Decomposition into Gaussian components simplifies complex lensing calculations.

The approach is computationally efficient and broadly applicable.

Abstract

We demonstrate an efficient method to compute the strong-gravitational-lensing deflection angle and magnification for any elliptical surface-density profile. This method solves a numerical hurdle in lens modelling that has lacked a general solution for nearly three decades. The hurdle emerges because it is prohibitive to derive analytic expressions of the lensing quantities for most elliptical mass profiles. In our method, we first decompose an elliptical mass profile into Gaussian components. We introduce an integral transform that provides us with a fast and accurate algorithm for the Gaussian decomposition. We derive analytic expressions of the lensing quantities for a Gaussian component. As a result, we can compute these quantities for the total mass profile by adding up the contributions from the individual components. This lensing analysis self-consistently completes the kinematic description in terms of Gaussian components presented by Cappellari (2008). Our method is general without extra computational burden unlike other methods currently in use.

Equations

$$\Sigma (x,y)\approx \sum _{j=0}^J{{\Sigma }_0}_j\exp \left(-\frac{q^2{x}^2+y^2}{2{\sigma }_j^2}\right),$$

$$F(x)\approx \sum _{n=0}^N{A}_n\exp \left(-\frac{x^2}{2{\sigma }_n^2}\right),$$

$$F(x)\equiv {\int }_0^{\infty }\frac{f(\sigma )}{\sqrt{2\text{uppi}}\sigma }\exp \left(-\frac{x^2}{2{\sigma }^2}\right)\mathrm{d}\sigma .$$

$$f(\sigma )=\frac{1}{\mathrm{i}{\sigma }^2}\sqrt{\frac{2}{\text{uppi}}}{\int }_{C}zF(z)\exp \left(\frac{z^2}{2{\sigma }^2}\right)\mathrm{d}z.$$

$$f(\sigma )\approx \sum _{n=0}^{2P}{\eta }_n\mathrm{Re}\left[F\left(\sigma {\chi }_n\right)\right].$$

$$\begin{array}{rl}& {\chi }_n={\left[\frac{2P\log (10)}{3}+2\text{uppi}\mathrm{i}n\right]}^{1/2},\\ & {\eta }_n=(-1{)}^{n}2\sqrt{2\text{uppi}}10^{P/3}{\xi }_p,\\ & {\xi }_0=\frac{1}{2},{\xi }_n=1,1\le n\le P,{\xi }_{2P}=\frac{1}{2^P},\\ & {\xi }_{2P-n}={\xi }_{2P-n+1}+2^{-P}\left(\genfrac{}{}{0ex}{}{P}{n}\right),0<n<P.\end{array}$$

$$\begin{array}{rl}F(x)& =\frac{1}{\sqrt{2\text{uppi}}}{\int }_0^{\infty }f(\sigma )\exp \left(-\frac{x^2}{2{\sigma }^2}\right)\mathrm{d}(\log \sigma )\\ \Rightarrow F(x)& =\underset{N\to \infty }{\lim }\sum _{n=1}^N\frac{f({\sigma }_n)}{\sqrt{2\text{uppi}}}\exp \left(-\frac{x^2}{2{\sigma }_n^2}\right)\Delta (\log \sigma {)}_n\\ \Rightarrow F(x)& \approx \sum _{n=1}^N{A}_n\exp \left(-\frac{x^2}{2{\sigma }_n^2}\right).\end{array}$$

$$\Sigma_{\rm NFW}(R)=\frac{\rho_{\rm s}r_{\rm s}}{(R/r_{\rm s})^{2}-1}\times\left\{\begin{array}[]{ll}1-\frac{\sec^{-1}\left(R/r_{\rm s}\right)}{\sqrt{(R/r_{\rm s})^{2}-1}}&(R>r_{\rm s}),\\ \frac{1}{3}\left[(R/r_{\rm s})^{2}-1\right]&(R=r_{\rm s}),\\ 1-\frac{{\rm sech}^{-1}{(R/r_{\rm s})}}{\sqrt{1-(R/r_{\rm s})^{2}}}&(R<r_{\rm s})\end{array}\right.$$

$${\Sigma }_{\text{Seˊrsic}}(R)={\Sigma }_{\mathrm{eff}}\exp \left[-b_n\left\{(R/R_{\mathrm{eff}}{)}^{1/n_{\text{Seˊrsic}}}-1\right\}\right]$$

$$\begin{array}{r}F(x,y)={\int }_0^{\infty }\mathrm{d}{\sigma }_1{\int }_0^{\infty }\mathrm{d}{\sigma }_2\frac{f({\sigma }_1,{\sigma }_2)}{2\text{uppi}{\sigma }_1{\sigma }_2}\exp \left(-\frac{x^2}{2{\sigma }_1^2}-\frac{y^2}{2{\sigma }_2^2}\right).\end{array}$$

$$\begin{array}{rl}F(x,y)& =F\left(R(q)\right)={\int }_0^{\infty }F\left(R(\varrho )\right)\delta (\varrho -q)\mathrm{d}\varrho \\ & ={\int }_0^{\infty }\mathrm{d}\varrho \delta (\varrho -q){\int }_0^{\infty }\mathrm{d}\sigma \frac{f_y(\sigma )}{\sqrt{2\text{uppi}}\sigma }\exp \left(-\frac{R(\varrho {)}^2}{2{\sigma }^2}\right),\end{array}$$

$$\begin{array}{rl}F(x,y)& =\frac{1}{2\text{uppi}}{\int }_0^{\infty }\mathrm{d}{\sigma }_1{\int }_0^{\infty }\mathrm{d}{\sigma }_2\frac{\sqrt{2\text{uppi}}q\delta ({\sigma }_2/{\sigma }_1-q)f_y({\sigma }_2)}{{\sigma }_1{\sigma }_2}\\ & \times \exp \left(-\frac{x^2}{2{\sigma }_1^2}-\frac{y^2}{2{\sigma }_2^2}\right).\end{array}$$

$$f({\sigma }_1,{\sigma }_2)=\sqrt{2\text{uppi}}q\delta \left(\frac{{\sigma }_2}{{\sigma }_1}-q\right)f_y({\sigma }_2).$$

$$F(x,y)\approx \sum _{n=1}^N{A}_n\exp \left(-\frac{q^2{x}^2+y^2}{2{\sigma }_n^2}\right).$$

$$\Sigma (x,y)\approx \sum _{j=1}^J{{\Sigma }_0}_j\exp \left(-\frac{q_j^2{x}^2+y^2}{2{\sigma }_j^2}\right),$$

$$I(x,y)\approx \sum _{k=1}^K{I_0}_k\exp \left(-\frac{q_k^2{x}^2+y^2}{2{\sigma }_k^2}\right).$$

$$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \iota & -\sin \iota \\ 0& \sin \iota & \cos \iota \end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right).$$

$$\rho (R^{\prime },z^{\prime })=\sum _{j=1}^J\frac{{q_j^{\prime }}^2{{\Sigma }_0}_j}{\sqrt{2\text{uppi}}{\sigma }_j{q}_j}\exp \left(-\frac{{q_j^{\prime }}^2{R^{\prime }}^2+{z^{\prime }}^2}{2{\sigma }_j^2}\right),$$

$$l(R^{\prime },z^{\prime })=\sum _{k=1}^K\frac{{q_k^{\prime }}^2{I_0}_k}{\sqrt{2\text{uppi}}{\sigma }_k{q}_k}\exp \left(-\frac{{q_k^{\prime }}^2{R^{\prime }}^2+{z^{\prime }}^2}{2{\sigma }_k^2}\right).$$

$$q^{\prime }=\frac{\sqrt{q^2-{\cos }^2\iota }}{\sin \iota }.$$

$$\begin{array}{rl}& \frac{bl\overline{v_{z^{\prime }}^2}-l\overline{v_{{\varphi }^{\prime }}^2}}{R^{\prime }}+\frac{\partial \left(bl\overline{v_{z^{\prime }}^2}\right)}{\partial R^{\prime }}=-l\frac{\partial \Phi }{\partial R^{\prime }},\\ & \frac{\partial \left(l\overline{v_{z^{\prime }}^2}\right)}{\partial z^{\prime }}=-l\frac{\partial \Phi }{\partial z^{\prime }}.\end{array}$$

$${\beta }_{z^{\prime }}(R^{\prime },z^{\prime })\equiv 1-\frac{\overline{v_{z^{\prime }}^2}}{\overline{v_{R^{\prime }}^2}}\approx 1-\frac{{\sum }_k{l}_k}{{\sum }_k{b}_k{l}_k}$$

$$\begin{array}{rl}\overline{v_{\mathrm{los}}^2}(x,y)=& 2\sqrt{\text{uppi}}G{\int }_0^1\sum _{j=1}^J\sum _{k=1}^K\frac{{q_j^{\prime }}^3{q_k^{\prime }}^2{{\Sigma }_0}_j{I_0}_k{u}^2}{{\sigma }_j{q}_j{\sigma }_k{q}_k}\\ & \times \frac{{\sigma }_k^2\left({\cos }^2\iota +b_k{\sin }^2\iota \right)+\mathcal{D}{x}^2{\sin }^2\iota }{\left(1-\mathcal{C}{u}^2\right)\sqrt{\left(\mathcal{A}+\mathcal{B}{\cos }^2\iota \right)\left[1-\left(1-{q_j^{\prime }}^2\right)u^2\right]}}\\ & \times \exp \left(-\mathcal{A}\left[x^2+\frac{(\mathcal{A}+\mathcal{B})y^2}{\mathcal{A}+\mathcal{B}{\cos }^2\iota }\right]\right)\mathrm{d}u,\end{array}$$

$$\begin{array}{rl}\mathcal{A}& =\frac{1}{2}\left(\frac{u^2{q_j^{\prime }}^2}{{\sigma }_j^2}+\frac{{q_k^{\prime }}^2}{{\sigma }_k^2}\right),\\ \mathcal{B}& =\frac{1}{2}\left\{\frac{1-{q_k^{\prime }}^2}{{\sigma }_k^2}+\frac{{q_j^{\prime }}^2\left(1-{q_j^{\prime }}^2\right)u^4}{{\sigma }_j^2\left[1-\left(1-{q_j^{\prime }}^2\right)u^2\right]}\right\},\\ \mathcal{C}& =1-{q_j^{\prime }}^2-\frac{{q_j^{\prime }}^2{\sigma }_k^2}{{\sigma }_j^2},\\ \mathcal{D}& =1-b_k{q_k^{\prime }}^2-\left[(1-b_k)\mathcal{C}+(1-{q_j^{\prime }}^2)b_k\right]u^2\end{array}$$

$$\rho (r)=\sum _{j=1}^J\frac{{{\Sigma }_0}_j}{\sqrt{2\text{uppi}}{\sigma }_j{q}_j}\exp \left(-\frac{r^2}{2{\sigma }_j^2}\right),$$

$$l(r)=\sum _{k=1}^K\frac{{I_0}_k}{\sqrt{2\text{uppi}}{\sigma }_k{q}_k}\exp \left(-\frac{r^2}{2{\sigma }_k^2}\right).$$

$$M(r)=\sum _{j=1}^J\frac{2\text{uppi}{\sigma }_j^2{{\Sigma }_0}_j}{q_j}\left[\mathrm{erf}\left(\frac{r}{\sqrt{2}{\sigma }_j}\right)-\sqrt{\frac{2}{\text{uppi}}}\frac{r}{{\sigma }_j}\exp \left(-\frac{r^2}{2{\sigma }_j^2}\right)\right],$$

$$\frac{\mathrm{d}\left(l\overline{v_r^2}\right)}{\mathrm{d}r}+\frac{2\beta l\overline{v_r^2}}{r}=-l\frac{\mathrm{d}\Phi }{\mathrm{d}r},$$

$$\beta (r)=1-\overline{v_{\theta }^2}/\overline{v_r^2}.$$

$$\overline{v_{\mathrm{los}}^2}(x,y)=\frac{2G}{I(x,y)}{\int }_{\sqrt{x^2+y^2}}^{\infty }{\mathcal{K}}_{\beta }\left(\frac{r}{\sqrt{x^2+y^2}}\right)l(r)M(r)\frac{\mathrm{d}r}{r}$$

Full text

††pubyear: 2019††pagerange: Unified lensing and kinematic analysis for any elliptical mass profile–References

1 Introduction: precise computation made easy

Gravitational lensing (hereafter, lensing) has versatile applications in astrophysics and cosmology. Lensing is the effect when light bends while passing by a massive object. If two galaxies sit along the same line-of-sight of an observer, the background galaxy appears multiple times due to lensing. This system is called a galaxy-scale strong-lensing system (hereafter, lens). Lenses are useful to measure the Hubble constant $H_{0}$, dark matter subhalo mass function, dust characteristics in galaxies, mass of super-massive black holes, stellar initial mass function, etc. (e.g., Falco et al., 1999; Peng et al., 2006; Suyu et al., 2013; Vegetti et al., 2014; Schechter et al., 2014).

In several applications, stellar kinematics play a complementary role to lensing. Lensing-only observables suffer from mass-sheet degeneracy (MSD) – if we appropriately rescale a mass profile after adding an infinite mass-sheet on top of it, then all the lensing observables stay invariant except the time delay (Falco et al., 1985; Schneider & Sluse, 2014). When a measured quantity depends on the mass profile, the MSD adds uncertainty to the measurement. Kinematics help break this degeneracy. Lensing observables probe the projected mass; kinematic observables probe the three-dimensional potential. Thus, the lensing–kinematics combination tightly constrains the mass profile, and enables us to robustly measure astrophysical and cosmological quantities (e.g., Treu & Koopmans, 2004; Barnabè et al., 2011).

In practice, it is difficult to compute lensing and kinematic quantities for elliptical mass profiles, which are required to describe the most common kind of lenses, i.e., elliptical galaxies. We need to integrate over density profiles while fitting a model to either the lensing or the kinematic data. For example, in lensing, the deflection angle is an integral over the surface density profile; in kinematics, the line-of-sight velocity dispersion is a double integral over three-dimensional mass and light profiles. For most elliptical density profiles, we are unable to express these integrals with elementary or special functions. Special functions are numerically well-studied, hence fast algorithms to compute them are usually available. Thus, these functions can be numerically convenient for evaluating model-predicted observables in large numbers (e.g., $\sim 10^{6}$) when sampling the lens model posterior with Markov chain Monte Carlo (MCMC) methods, or when only searching for the best-fit model even. Otherwise, it would be inefficient to numerically compute lensing and kinematic integrals for general elliptical profiles.

Usually, the numerical difficulties are circumvented through simplifying assumptions or approximations. We now outline common assumptions and approximations in kinematic and lensing analyses, noting that these can limit the accuracy and precision of the inferences.

Either axisymmetry or spherical symmetry is usually assumed for kinematic analysis. If we start with a surface density profile for lensing analysis, we need to deproject this profile along the line of sight to compute the kinematics. This deprojection has an infinite degeneracy (Contopoulos, 1956). Therefore, it is necessary to choose a line-of-sight symmetry when deprojecting, for which either axisymmetry or spherical symmetry is often a convenient choice.

In lensing analysis, spherical symmetry is rarely sufficient and we need to consider ellipticity to achieve the required precision. All the lensing quantities are related to the deflection potential or its derivatives. The time delay depends on the deflection potential difference. The gradient of the deflection potential gives the deflection angle. The Hessian of the deflection potential relates to the surface density, the shear, and the magnification. To efficiently compute these quantities for the elliptical case, we can find the following three approximations in the literature:

Ellipticity in deflection potential: The gradient and the Hessian of an elliptical deflection potential can be easily computed through numerical differentiation (e.g., Kovner, 1987; Golse & Kneib, 2002). However, this solution is not general, as the surface density becomes dumbbell-shaped for an elliptical deflection potential with axis ratio $q\lesssim 0.6$ (Fig. 1, Kassiola & Kovner 1993). This oddly shaped surface-density is unphysical. 2. 2.

Elliptical power-law profile: We can efficiently compute the lensing quantities for the elliptical power-law profile using numerical approximation or analytical expressions (for the isothermal case, Kormann et al. 1994; and for the general case, Barkana 1998; Tessore & Metcalf 2015). This profile can be sufficient to use in statistical studies that do not require detailed modelling of individual lenses (e.g., Koopmans et al., 2009; Sonnenfeld et al., 2013). However, adopting a power-law profile artificially breaks the MSD and it could potentially bias the $H_{0}$ measurement (e.g., Schneider & Sluse, 2013; Sonnenfeld, 2018). Therefore, we need to explore different, physically-motivated mass models, such as a composite model that explicitly accounts for the luminous and the dark components (Suyu et al., 2014; Yıldırım et al., 2019). 3. 3.

Chameleon profile: The Sérsic profile well describes the surface brightness of a galaxy (Sérsic, 1968). The Chameleon profile approximates the Sérsic profile within a few per cent (Dutton et al., 2011). We can efficiently compute lensing quantities for the Chameleon profile using analytic expressions. However, this profile only describes the baryonic component. The precise lensing analysis of elliptical Navarro–Frenk–White (NFW) profile for the dark component still lacks a general solution (Navarro et al., 1997).

In a nutshell, these approximations for elliptical lensing analysis are only applicable in restricted regimes as described above.

In this paper, we present a general method to precisely compute the gradient and the Hessian of the deflection potential for any elliptical surface-density profile. The method follows a “divide and conquer” strategy. We can approximately divide, or decompose, an elliptical profile into Gaussian components (e.g., Bendinelli, 1991). For this Gaussian decomposition, we devise a fast and accurate algorithm by introducing an integral transform. For each Gaussian component, we derive analytic expressions of the gradient and the Hessian of the deflection potential. For the deprojected Gaussian component, Cappellari (2008) derives the line-of-sight velocity dispersion. We can combine the computed quantities from each Gaussian component back together to obtain these quantities for the total density profile. In this way, the lensing and kinematic descriptions are self-consistently unified. At the same time, this method is general, as we can apply it to lensing with any elliptical surface-density and to kinematics with either axisymmetry or spherical symmetry. Our method is more efficient than numerical integration to compute lensing quantities for an elliptical surface-density profile.

We organize this paper as follows. In Section 2, we motivate the “divide and conquer” strategy behind our method and introduce an integral transform that provides a fast algorithm for decomposing any elliptical surface-density profile into Gaussian components. In Section 3, we summarize the kinematic description of the Gaussian components from Cappellari (2008). In Section 4, we derive the gradient and the Hessian of the deflection potential for elliptical Gaussian surface-density. Next in Section 5, we demonstrate a proof-of-concept for our method using simulated data. Then, we summarize the paper in Section 6. Additionally in Appendix A, we prove some fundamental theorems for the integral transform introduced in Section 2. This integral transform with a Gaussian kernel is the continuous case of Gaussian decomposition. The theorems in Appendix A establish the existence, uniqueness, and invertibility of this transform.

2 Decomposing an elliptical profile into Gaussian components

To make the computation tractable, we aim to decompose an elliptical surface-density profile into simpler functions. This function should be simple enough so that we can both

express the deflection angle in terms of elementary or special functions, and 2. 2.

easily deproject it into three-dimension and compute the enclosed mass for kinematic analysis.

The Gaussian function meets both of these criteria. We validate criterion (i) in Section 4.2, where we express the deflection angle for an elliptical Gaussian surface-density profile with the complex error function. For criterion (ii), deprojecting a two-dimensional Gaussian into three-dimension is straightforward, as the Abel inversion of a two-dimensional Gaussian is a three-dimensional Gaussian. The enclosed mass for a three-dimensional Gaussian has the form of the error function, which we can efficiently compute without integrating numerically. Given these points, we approximately decompose an elliptical surface-density profile as

[TABLE]

where ${\Sigma_{0}}_{j}$ is the amplitude of the $j$-th Gaussian, and all the components have a common axis ratio $q$. Similar decomposition into Gaussian components has been used in the literature to fit the surface brightness profile of galaxies [called as the multi-Gaussian expansion (MGE) by Emsellem et al. (1994), and the mixture-of-Gaussians by Hogg & Lang (2013)]. The lensing and kinematic quantities of our interest – namely the gradient and the Hessian of the deflection potential, and the line-of-sight velocity dispersion – follow the principle of superposition. As a result, we can compute these quantities separately for each Gaussian component and then add them together to recover these quantities for the total surface-density profile (see Fig. 2).

We describe the kinematics and lensing analyses for the Gaussian components in Sections 3 and 4, but first, we need a fast method to decompose an elliptical profile into Gaussian components. Cappellari (2002) presents a method that uses non-linear optimization. However, this non-linear optimization method is computationally too expensive to implement within MCMC. Although, this method has been implemented to compute the kinematic observable while sampling from the lens model posterior (e.g., Birrer et al., 2019). Computing the kinematic observable can involve at least one numerical integration, which is the main bottleneck in the efficiency, not the non-linear optimization method for Gaussian decomposition. To make our lens-modelling method efficient, we require a (i) general, (ii) precise, and (iii) fast technique to decompose a function into Gaussians. In Sections 2.1 and 2.2, we provide a technique that satisfies these three requirements.

2.1 An integral transform for fast Gaussian decomposition

Now, we introduce an integral transform with a Gaussian kernel. Using this transform, we obtain an algorithm to efficiently decompose an elliptical surface-density profile into Gaussian components.

We start with the simple one-dimensional case of the integral transform. We aim to approximate a function $F(x)$ as a sum of Gaussian components as

[TABLE]

where $A_{n}$ and $\sigma_{n}$ are respectively the amplitude and the standard deviation of the $n$-th Gaussian component. We can convert this discrete summation into a continuous integral by taking $N\to\infty$. Accordingly, we define the following integral transform:

[TABLE]

Here, the amplitude $A_{n}$ is converted into a function $f(\sigma)/\sqrt{2\uppi}\sigma$. We call $F(x)$ as the transform of $f(\sigma)$. We prove three fundamental properties of this integral transform in Appendix A. These three properties tell us that

this integral transform exists for most mass and light profiles of practical use, 2. 2.

the transform is unique for these functions, and 3. 3.

the integral transform is invertible.

We call $f(\sigma)$ as the inverse transform of $F(x)$. The inverse transform is given by

[TABLE]

Here, $\mathrm{i}$ is the imaginary unit as $\mathrm{i}=\sqrt{-1}$. Also, we have extended $F(x)$ to some region on the complex plane and wrote it as $F(z)$, where $z$ is a complex variable. The contour $C$ for the integral lies within the region where $F(z)$ is defined (for details, see Appendix A and Fig. 7). In Section 2.1.1, we provide an algorithm that does not require $C$ to be explicitly specified for computing the inverse transform $f(\sigma)$.

We can use the inverse transform to decompose a function into Gaussian components and the forward transform to recover the original function by combining the Gaussian components. We first provide an efficient algorithm to compute the inverse transform from equation (2.4), then we discuss a method for computing the forward transform from equation (2.3).

2.1.1 Computing the inverse transform

The integral transform in equation (2.3) can be converted into a Laplace transform by suitable change of variables (Remark A.9). Therefore, we can use any of the several algorithms available for inverse Laplace transform by appropriately changing the variables (for a simple overview of the algorithms, see Abate & Whitt, 2006). In this paper, we modify the Euler algorithm to approximate equation (2.4) as

[TABLE]

(Abate et al., 2000). Here, the weights $\eta_{n}$ and nodes $\chi_{n}$ can be complex-valued and they are independent of $f(\sigma)$. The weights and the nodes are given by

[TABLE]

We can precompute the weights and the nodes just once before the MCMC sampling. In that way, computing them does not add any extra burden in computing the likelihood. The precision of the inverse transform is $\sim\mathcal{O}(10^{-0.6P})$ (Abate & Whitt, 2006). Therefore, the value of $P$ can be appropriately chosen to achieve a required precision. Note that the decimal precision of the machine sets an effective upper limit for $P$. For example, the precision will not improve with increasing $P$ when $P\gtrsim 12$ for 32-bit floating point number, and when $P\gtrsim 27$ for 64-bit floating point number. Thus, equation (2.5) gives a straightforward, fast, and precise algorithm to compute the inverse transform.

2.1.2 Computing the forward transform

Let us approximate the forward transform integral such that we can recover $F(x)$ from only a finite number of $f(\sigma)$ values computed at fixed $\sigma$’s. This finite number should be on the order of tens to keep the lensing analysis manageable, as we have to compute lensing quantities for each Gaussian component individually. We write equation (2.3) as

[TABLE]

We have recovered the form of equation (2.2) by taking logarithmically spaced $\sigma_{n}$. Here, the amplitudes are $A_{n}=w_{n}f(\sigma_{n})\Delta(\log\sigma)_{n}/\sqrt{2\uppi}$. The weights $w_{n}$ depend on the choice of the numerical integration method. We use the trapezoidal method with weights $w_{1}=0.5$, $w_{n}=1$ for $1<n<N$, $w_{N}=0.5$, as the trapezoidal method is highly efficient to numerically compute integrals of this form (Goodwin, 1949). As a result, equation (2.7) efficiently recovers $F(x)$ from only a finite number of $f(\sigma)$ values.

As an example, we demonstrate the integral transform method to decompose the NFW profile and the Sérsic profile into Gaussian components (Fig. 3). The two-dimensional projected NFW profile is given by

[TABLE]

(Bartelmann, 1996). Here, $\rho_{\rm s}$ is the three-dimensional density normalization, and $r_{\rm s}$ is the scale radius. The Sérsic profile is given by

[TABLE]

(Sérsic, 1968). Here, the normalizing factor $b_{n}$ ensures that half of the total projected mass is contained within the effective radius $R_{\rm eff}$. We only need to compute as many $f(\sigma)$ values as the number of Gaussian components. We can appropriately choose this number to achieve the required precision for approximating the original function within a given range of $R$. In this example, we set $r_{\rm s}=5R_{\rm eff}$ and assume 1 per cent Poisson noise at $R=R_{\rm eff}$. Then, we can approximate both the projected NFW profile and the Sérsic function within the noise level with only 15 Gaussian components in the range $0.1R_{\rm eff}\leq R\leq 10R_{\rm eff}$. The standard deviations $\sigma_{n}$ of the 15 Gaussians are logarithmically spaced between $0.005r_{\rm s}$ and $50r_{\rm s}$ for the NFW profile, and between $0.02R_{\rm eff}$ and $15R_{\rm eff}$ for the Sérsic profile. Thus using the integral transform method, we can decompose a function into Gaussian components within any required precision by appropriately choosing the component number $N$.

2.1.3 The integral transform method is more efficient than the MGE method.

We compare our method to decompose a one-dimensional function into Gaussian components with the MGE method (Fig. 4). We use both methods to decompose the Sérsic function into 15 Gaussian components. The MGE method approximates the Sérsic function within the noise level up to $\sim 6R_{\rm eff}$, whereas our method approximates the Sérsic function within the noise level up to 10$R_{\rm eff}$ with the same number of components. Albeit, we can increase the number of Gaussian components in the MGE method to reach the desired precision within a given radius. In our method, the precision of the decomposition can be affected by both $P$ in equation (2.5) and the number of Gaussians $N$ in equation (2.7). However, if $P$ is appropriately chosen so that $10^{-0.6P}$ is sufficiently (e.g., by a factor of $10^{-2}$–$10^{-4}$) smaller than the required precision, then the precision predominantly depends on $N$. For lensing and kinematic analyses, increasing the number of Gaussian components $N$ introduces more computational burden than increasing $P$. Therefore, it is advisable to first choose a sufficiently large $P$ and then adjust the number of Gaussians to achieve the required precision. Note that the real power of our method is in its efficiency. A python implementation of our method is $\sim 10^{3}$ times faster than the MGE method to decompose a one-dimensional function into Gaussian components with similar or better precision.

2.2 Decomposing a two-dimensional elliptical profile with the one-dimensional transform

So far, we have discussed the one-dimensional case of the integral transform; now we show that the one-dimensional transform is sufficient to decompose a two-dimensional elliptical profile. We can extend the one-dimensional integral transform from equation (2.3) into a two-dimensional integral transform for a function $f(\sigma_{1},\sigma_{2})$ as

[TABLE]

If $F(x,y)$ is elliptically symmetric, then we can express it as $F(R)$ in terms of the elliptical radius $R=\sqrt{q^{2}x^{2}+y^{2}}$ and axis ratio $q$. Then, we can write

[TABLE]

where $f_{y}(\sigma)$ is the inverse transform of $F(0,y)$. If we make the change of variables $\sigma_{1}=\sigma/\varrho$, $\sigma_{2}=\sigma$, this integral becomes

[TABLE]

Because of the uniqueness property (Theorem A.7), comparing equations (2.10) and (2.12) we can write

[TABLE]

Therefore, for an elliptically symmetric function, it is sufficient to numerically compute the one-dimensional inverse transform $f_{y}(\sigma_{2})$ along the $y$-axis. As a result, we can express a two-dimensional elliptical function as a sum of elliptical Gaussian components as

[TABLE]

By now, we have shown that the integral transform method meets all of our three requirements for decomposing a function into Gaussian components:

generality: the method applies to most mass and light profiles of practical use, 2. 2.

precision: the method achieves any required precision over a given range by appropriately choosing the number of Gaussian components, and 3. 3.

efficiency: the method runs approximately $\sim 10^{3}$ times faster than the previously available method.

With these three requirements met, our lensing analysis method has cleared the first hurdle to be feasible in practice. In Section 4.2, we show that we can efficiently compute the gradient and the Hessian of the deflection potential for an elliptical Gaussian surface-density. With that, our method also clears the final hurdle to be efficient. Next in Section 3, we summarize the kinematic analysis for the Gaussian components from Cappellari (2008). Then in Section 4, we describe the lensing analysis for the Gaussian components and complete the unification of lensing and kinematic descriptions.

3 Kinematics of Gaussian components

Cappellari (2008) presents the Jeans anisotropic modelling of kinematics for a mass profile decomposed into Gaussian components. We summarize the analysis here to complete our unified framework. The kinematic observable is the luminosity-weighted, line-of-sight velocity dispersion. The velocity dispersion can be an integrated measurement within a single aperture or it can be spatially resolved on the plane of the sky. To compute this quantity for a combination of mass and light profiles, we need to solve the Jeans equations. We can decompose the surface mass-density profile into Gaussian components as

[TABLE]

and decompose the surface brightness profile into Gaussian components as

[TABLE]

Here, we use different subscript letters to make the context of the Gaussian decomposition clear: we use the subscript $j$ for a component of the mass profile and the subscript $k$ for a component of the light profile. We have also allowed different ellipticity for each Gaussian component represented by $q_{j}$ or $q_{k}$. This is the most general case, for example, when structures with different ellipticities constitute the total mass or light distribution. We first need to deproject these two-dimensional profiles into three-dimension as the kinematic quantities depend on the three-dimensional distributions of mass and light. We assume axisymmetry or spherical symmetry for the deprojected three-dimensional structure to circumvent the infinite degeneracy in deprojection. First, we provide the kinematic analysis for the axisymmetric case in Section 3.1; then we do the same for the simpler case of spherical symmetry in Section 3.2.

3.1 Axisymmetric case

For an axisymmetric system, the cylindrical coordinates $(R^{\prime},z^{\prime},\phi^{\prime})$ are the most suitable to express the Jeans equations. We use the prime symbol to denote the coordinates in the system’s symmetry-frame, where the $z^{\prime}$-axis aligns with the axis of symmetry. We assign $(x,y,z)$ coordinates to the sky frame, where the $z$-axis aligns with the line of sight and the $x$-axis aligns with the projected major axis. If the galaxy is inclined by an angle $\iota$, then the $(x^{\prime},y^{\prime},z^{\prime})$ coordinates in the symmetry frame relate to the the $(x,y,z)$ coordinates in the sky frame as

[TABLE]

In the symmetry frame, equation (3.1) deprojects into the mass density profile $\rho$ as

[TABLE]

and equation (3.2) deprojects into the light density $l$ as

[TABLE]

Here, the intrinsic axis ratio $q^{\prime}$ relates to the projected axis ratio $q$ as

[TABLE]

We can first solve the Jeans equations for these mass and light density profiles in the symmetry frame to get the intrinsic velocity dispersions and then integrate along the line of sight to obtain the line-of-sight velocity dispersion.

The Jeans equations to solve for an axisymmetric system are

[TABLE]

Here, the gravitational potential $\Phi$ relates to the three-dimensional mass density $\rho$ by $\nabla^{2}\Phi=\rho$, and $b$ represents the anisotropy as in $\overline{v_{R^{\prime}}^{2}}=b\overline{v_{z^{\prime}}^{2}}$. We can let $b_{k}$ for each luminous Gaussian component have different values to approximate the luminosity-weighted anisotropy parameter as

[TABLE]

(Binney & Mamon, 1982; Cappellari, 2008). The line-of-sight second velocity moment for total mass profile obtained from solving the Jeans equations is given by

[TABLE]

where $G$ is the gravitational constant and

[TABLE]

(Cappellari, 2008). The line-of-sight velocity dispersion $\sigma_{\rm los}$ relates to the second velocity moment by $\overline{v^{2}_{\rm los}}=v_{\rm mean}^{2}+\sigma^{2}_{\rm los}$, where $v_{\rm mean}$ is the stellar mean velocity.

3.2 Spherical case

If we assume the system is spherically symmetric, then the spherical coordinates $(r,\phi,\theta)$ are the most suitable to express the Jeans equations. In this coordinate system, the mass density profile deprojected from equation (3.1) takes the form

[TABLE]

and the light density profile deprojected from equation (3.2) turns into

[TABLE]

The projected axis ratio $q$ shows up in these equations to keep the total mass and luminosity conserved. We can express the three-dimensional enclosed mass for this density profile as

[TABLE]

where $\mathrm{erf}\ (x)$ is the error function. The spherical Jeans equation is

[TABLE]

where $\beta(r)$ is the anisotropy parameter given by

[TABLE]

Spherical symmetry imposes that $\overline{v_{\theta}^{2}}=\overline{v_{\phi}^{2}}$. By solving the Jeans equation for the spherically symmetric case, we can obtain the line-of-sight second velocity moment as

[TABLE]

(Mamon & Łokas, 2005). Here, the function $\mathcal{K}_{\beta}(\upsilon)$ depends on the form of the anisotropy parameter $\beta(r)$. For the isotropic case with $\beta=0$, the function $\mathcal{K}_{\beta}$ shapes into

[TABLE]

For the Osipkov-Merritt parameterization $\beta(r)=r^{2}/(r^{2}+r_{\rm ani}^{2})$, where $r_{\rm ani}$ is a scale radius, the function $\mathcal{K}_{\beta}$ takes the form

[TABLE]

with $\upsilon_{\rm ani}=r_{\rm ani}/\sqrt{x^{2}+y^{2}}$ (Osipkov, 1979; Merritt, 1985a, b). See equation (A16) of Mamon & Łokas (2005) for the form of $\mathcal{K}_{\beta}$ corresponding to other parameterizations of $\beta(r)$. When assuming spherical symmetry is sufficient, we can use equation (3.16) to compute the line-of-sight velocity dispersion in a much simpler way than the axisymmetric case [cf. equation (3.9)].

The kinematic description of an elliptical mass distribution by decomposing it into Gaussian components is thus well developed in the literature. In the next section, we unify the lensing description with the kinematic description under the same framework.

4 Lensing by Gaussian components

In this section, we present the lensing analysis for an elliptical surface-density profile decomposed into Gaussian components. In Section 2, we introduced an integral transform that efficiently decomposes an elliptical surface-density profile into Gaussian components as

[TABLE]

We can compute a lensing quantity for each individual Gaussian component, and then linearly add the contributions from all the components to obtain the total lensing quantity. For example, if $\mn@boldsymbol{\alpha}_{j}(x,y)$ is the deflection at position $(x,y)$ for the $j$-th Gaussian component, then the total deflection is simply given by $\mn@boldsymbol{\alpha}(x,y)=\sum_{j=1}^{J}\mn@boldsymbol{\alpha}_{j}(x,y)$. Therefore, it is sufficient to analyze the lensing properties of one elliptical Gaussian component. We use the complex formulation of lensing to solve the deflection integral for an elliptical Gaussian surface-density profile. Below, we first lay out the complex formalism of lensing in Section 4.1; then we study the lensing properties of an elliptical Gaussian surface-density profile in Section 4.2.

4.1 Complex formulation of lensing

The strong lensing effect is usually described using the vector formulation on the two-dimensional image plane. We first define the lensing quantities in the familiar vector formulation, then we translate them to the complex formulation. The convergence $\kappa$ is a dimensionless surface-density defined as $\kappa\equiv\Sigma/\Sigma_{\rm crit}$, where the critical density $\Sigma_{\rm crit}$ is given by

[TABLE]

Here, $c$ is the speed of light. The three angular diameter distances are $D_{\rm d}$: between the observer and the deflector, $D_{\rm s}$: between the observer and the source, and $D_{\rm ds}$: between the deflector and the source. The convergence $\kappa$ relates to the vector deflection angle ${\mn@boldsymbol{\alpha}}$ as $\kappa=\nabla\cdot\mn@boldsymbol{\alpha}/2$. The deflection angle $\mn@boldsymbol{\alpha}$ is the gradient of the deflection potential as $\mn@boldsymbol{\alpha}=\nabla\psi$, thus the convergence $\kappa$ relates to the deflection potential $\psi$ as $\kappa=\nabla^{2}\psi/2$. The Hessian of the deflection potential is

[TABLE]

The convergence $\kappa$, the shear parameters $(\gamma_{1},\ \gamma_{2})$, and the magnification $\mu$ relate to the Hessian, since we can express them as

[TABLE]

where $\bm{\mathsf{I}}$ is the identity matrix. Therefore, if we start with a convergence $\kappa$ and derive the deflection $\mn@boldsymbol{\alpha}$ and the shear parameters $(\gamma_{1},\ \gamma_{2})$, then we can obtain the gradient and the Hessian of the deflection potential from them.

Now we reformulate the lensing quantities on the complex plane. Following Bourassa et al. (1973), we can express the deflection vector $\mn@boldsymbol{\alpha}$ as a complex quantity

[TABLE]

where the complex quantity $z=x+\mathrm{i}y$ corresponds to the position vector $\mathbb{r}=(x,\ y)$. We can define a complex deflection potential $\psi(z)$ with its real part equal to the usual deflection potential (Schramm, 1990). Then, the complex deflection angle is the Wirtinger derivative of the deflection potential as

[TABLE]

We can express the convergence $\kappa$ as

[TABLE]

Furthermore, the complex shear $\gamma\equiv\gamma_{1}+\mathrm{i}\gamma_{2}$ satisfies the relation

[TABLE]

Using this complex formulation, we analyze the lensing properties of an elliptical Gaussian convergence next in Section 4.2.

4.2 Lensing by elliptical Gaussian convergence

We derive the deflection angle and shear for the elliptical Gaussian convergence

[TABLE]

where $R=\sqrt{q^{2}x^{2}+y^{2}}$ is the elliptical radius. Using the complex formulation, the deflection angle for the elliptical convergence can be obtained from

[TABLE]

where $\mathrm{sgn}(z)\equiv\sqrt{z^{2}}/z$ is the complex sign function, and $R(z)=\sqrt{q^{2}x^{2}+y^{2}}$ is the semi-minor axis length for the ellipse with axis-ratio $q$ that goes through the point $z=x+iy$ (Bourassa & Kantowski, 1975; Bray, 1984). With changes of variables $s=1/2\sigma^{2}$, $t=(1-q^{2})/q^{2}z^{2}$, $\tau=\sqrt{1-t\zeta^{2}}$, we can express equation (4.10) as

[TABLE]

where $\mathrm{erfi}\ (z)\equiv-\mathrm{i}\ \mathrm{erf}\ (\mathrm{i}z)$ and we have defined the function

[TABLE]

We obtain the complex conjugate of the complex shear from equation (4.8) as

[TABLE]

Both the deflection angle and the shear contain the function $\varsigma(z;q)$. This function relates to the Faddeeva function $w_{\rm F}(z)$. First, we write the function $\varsigma(z;q)$ as

[TABLE]

where $\varpi(z;q)=\exp\left(-z^{2}\right)\ \rm{erfi}\left(qx+iy/q\right)$. We can express $\varpi(z;q)$ using the Faddeeva function $w_{\rm F}(z)$ as

[TABLE]

Thus, we can compute the deflection angle and the shear using the Faddeeva function (Fig. 5). Faddeeva function is a well-studied special function for its various applications in physics, for example, in radiative transfer and in plasma physics (e.g., Armstrong, 1967; Jiménez-Domínguez et al., 1989). We can readily compute $w_{\rm F}(z)$ in python using the scipy.special.wofz function. For some other popular programming languages, code-packages to compute this function are available at the web-address http://ab-initio.mit.edu/Faddeeva. In this paper, we use the algorithm outlined by Zaghloul (2017) to compute $w_{\rm F}(z)$ with relative error less than $4\times 10^{-5}$ over the whole complex plane. We state this algorithm in Appendix B. A python implementation of this algorithm is about twice as fast as the function provided by scipy. As a result, we can efficiently compute the gradient and Hessian of the deflection potential for an elliptical Gaussian convergence using equations (4.11) and (4.13).

Now, we turn our attention to computing the deflection potential. The deflection potential $\psi(z)$ is given by

[TABLE]

where we set $\psi(0)=0$. Often times we are interested in the potential difference between two points $z_{1}$ and $z_{2}$ given by

[TABLE]

This integral is independent of the choice of a contour. We have to carry out this integral numerically. However, the number of times we need to compute it in most applications, e.g., for computing time delays, is much fewer than that for $\alpha(z)$. We can also numerically solve the Poisson equation $\nabla^{2}\psi=2\kappa$ using the Fourier transform of the deflection potential $\hat{\psi}\equiv\mathcal{F}[\psi]$ (van de Ven et al., 2009). This equation turns into $k^{2}\hat{\psi}=2\hat{\kappa}$ in the Fourier domain. The solution of the Poisson equation is then $\psi=\mathcal{F}^{-1}(2\hat{\kappa}/k^{2})$. We can analytically compute the forward Fourier transform because the convergence has the Gaussian form. Then, we need to compute only the inverse transform numerically. Although obtaining the deflection potential necessitates a numerical integration or a numerical Fourier transform, we can keep the computational burden under control in most applications by computing this quantity only for a feasible number of models sampled from the lens model posterior.

5 Adding it all together: proof of concept

In this section, we demonstrate the feasibility of our method to model lenses. We first simulate synthetic data of a mock lens and then model the lens using our method. We use the publicly available lens-modelling software lenstronomy to simulate the synthetic data and perform the model-fitting (Birrer & Amara, 2018). We added extra modules to lenstronomy to implement the lensing analysis presented in Section 4.2.

For the mock strong lensing system, we adopt an elliptical NFW deflection potential for the dark component and an elliptical Chameleon convergence for the luminous component. We take realistic scale sizes and normalizations for these profiles. We choose the Chameleon profile for two reasons:

we can analytically simulate the data with ellipticity in the convergence, and 2. 2.

we know the Sérsic-profile parameters that approximates the chosen Chameleon profile a priori, so we can check the fidelity of our method.

We parameterize the scaling of the NFW profile with two parameters: scale radius $r_{\rm s}$ and the deflection angle $\alpha_{\rm s}$ at $r_{\rm s}$. For a spherical NFW profile given by

[TABLE]

the normalization $\rho_{s}$ relates to $\alpha_{s}$ as

[TABLE]

(Meneghetti et al., 2003). The elliptical Chameleon convergence is given by

[TABLE]

(Suyu et al., 2014). We also add external shear to the mass profiles. Therefore, our fiducial lens mass profile has three components in total: elliptical NFW deflection potential, elliptical Chameleon convergence, and external shear.

We simulate data for this fiducial lens system with image quality similar to the Hubble Space Telescope (HST) Wide-Field Camera 3 imaging in the F160W filter (see top-left panel of Fig. 6). We adopt 0.08 arcsec for the pixel size, 2197 s for exposure time, and a realistic point spread function (PSF) to achieve data quality similar to the lens sample presented by Shajib et al. (2019).

We fit the synthetic data with a model composed of elliptical NFW deflection potential, elliptical Sérsic convergence profile decomposed into Gaussians, and external shear. Note that we take ellipticity in the deflection potential for the NFW profile due to a design restriction of lenstronomy. However, we can also extend an elliptical NFW convergence into Gaussians for lensing analysis in principle. For now, this limitation does not affect the point of this exercise to show that lensing analysis with Gaussian components is feasible. We take the PSF as known during model-fitting for simplicity. We also separately model the lens with the fiducial mass profiles for comparison. In both cases, the parameters for the light and the luminous mass profiles are joint except for the amplitudes. We fit the model to the data by using the MCMC method. For every sample point in the parameter space, we first decompose the elliptical Sérsic profile into Gaussian components using equations (2.5) and (2.7). Similar to the example in Section 2.1.2, we take 15 Gaussian components with logarithmically spaced $\sigma$’s between $0.02R_{\rm eff}$ and $15R_{\rm eff}$. These 15 Gaussian components approximate the Sérsic function well within the noise level for $0.1R_{\rm eff}\leq R\leq 10R_{\rm eff}$ (Fig. 3). We compute the gradient and the Hessian of the deflection potential for each Gaussian component. Finally, we add the contributions from all the individual components together to obtain these quantities for the total mass profile. These total quantities are used to compute the likelihood in the MCMC method for fitting the model to the data.

Our method fits the synthetic data very well (see the ‘Normalized Residuals’ plot in Fig. 6). The fiducial Sérsic profile parameters are also recovered with reasonable to high accuracies at the same time (Table 1). The total runtime is only approximately three times longer than using the fiducial model with the Chameleon profile. This loss in efficiency is a reasonable tradeoff for generality. Thus, we have demonstrated a feasible implementation of our lens-modelling method.

6 Conclusion: precision is feasible.

In this paper, we present a general method for precise lensing analysis of any elliptical convergence profile. Our method follows a “divide and conquer” strategy. In our method, we first decompose an elliptical convergence profile into Gaussian components as

[TABLE]

We then compute lensing quantities, e.g., the gradient and the Hessian of the deflection potential, for each Gaussian component. Finally, we add the lensing quantities from individual Gaussian components together to obtain these quantities for the total surface-density profile. Moreover, we can straightforwardly deproject a Gaussian component to obtain its corresponding three-dimensional density profile assuming either axisymmetry or spherical symmetry. Then, we can also compute the kinematic properties, such as the line-of-sight velocity dispersion, for each Gaussian component (Cappellari, 2008). We can then add the velocity dispersions from individual Gaussians together to obtain the total line-of-sight velocity dispersion. In this way, we self-consistently unify the lensing and kinematic descriptions of any elliptical mass profile.

We introduce an integral transform with a Gaussian kernel that leads us to a general, precise, and fast algorithm for decomposing a surface density profile into Gaussian components. Without such an algorithm, decomposing into Gaussians would end up as a bottleneck in the lens modelling efficiency. We obtain the algorithm by first inverting the integral transform as

[TABLE]

Although this is an integral, we provide a straightforward formula to compute $f(\sigma)$. The computed values of $f(\sigma)$ then quantify the amplitudes $\kappa_{0j}$ of the Gaussian components in equation (6.1). As a result, this integral transform fulfills the three requirements for a decomposition algorithm to be (i) general, (ii) precise, and (iii) fast. To be specific, this decomposition algorithm is $\sim 10^{3}$ times faster than the MGE algorithm from Cappellari (2002). Consequently, our lensing analysis requires the same order of CPU time as other methods currently in use to model a lens with a composite mass profile. Thus, the integral transform enables the lens modelling with the Gaussian components to be efficient and, in turn, makes our unified framework for lensing and kinematic analysis of an elliptical mass profile feasible.

Our method enables precise lens modelling with an elliptical mass profile for several astrophysical applications. Specifically, our method gives an efficient method to model composite mass profiles with separate components for the baryon and the dark matter. For example, the usual choices for these components are the Sérsic and the NFW profiles; both are computationally difficult to directly implement in lens modelling for the elliptical case. Our method makes both of these profiles computationally tractable while achieving the required precision. Thus, our method will be useful in applications where a composite mass profile is essential for lens modelling, for example, in detecting dark-matter substructure, in measuring the Hubble constant, and in testing massive elliptical-galaxy formation theories (e.g., Vegetti et al., 2012; Wong et al., 2017; Nightingale et al., 2019).

Acknowledgements

AJS thanks Simon Birrer and Shouman Das for helpful discussions. AJS additionally thanks the anonymous referee for very useful suggestions that improved this paper. AJS also expresses gratitude to Adriano Agnello, Simon Birrer, Xuheng Ding, Xinnan Du, Abhimat K. Gautam, Briley Lewis, Michael Topping, and Tommaso Treu for providing feedbacks that greatly improved the writing and the presentation of this paper. AJS acknowledges support by National Aeronautics and Space Administration (NASA) through Space Telescope Science Institute grant HST-GO-15320.

This research made use of lenstronomy (Birrer & Amara, 2018), numpy (Oliphant, 2015), scipy (Jones et al., 2001), jupyter (Kluyver et al., 2016), and matplotlib (Hunter, 2007).

Appendix A Properties of the Integral transform with a Gaussian Kernel

In this appendix, we prove some fundamental properties of the integral transform with a Gaussian kernel. In three theorems, we prove that

the integral transform exists for a function with certain characteristics, 2. 2.

the transform is unique for a continuous function, and 3. 3.

the transform is invertible.

First, we define the integral transform.

Definition A.1

Define an integral transform $\mathcal{T}$ that takes a function $f(\sigma):\mathbb{R}_{\geq 0}\to\mathbb{R}$ to a function $F(z):\mathbb{C}\to\mathbb{C}$ as

[TABLE]

Next, we define a transformable function in the context of this paper.

Definition A.2

A function $f(\sigma):\mathbb{R}_{\geq 0}\to\mathbb{R}$ is said to be transformable, if it satisfies the following conditions:

the function $f(\sigma)$ is piecewise continuous, 2. 2.

the function $f(\sigma)=\mathcal{O}(\exp(c/2\sigma^{2}))$ as $\sigma\to 0$, where $c\in\mathbb{R}_{\geq 0}$, 3. 3.

the function $f(\sigma)=\mathcal{O}(\sigma^{\lambda})$ with $\lambda<0$ as $\sigma\to\infty$.

We refer to these three conditions as the transformability conditions. The namesake for the transformable function is made clear next in Theorem A.3.

Theorem A.3** (Existence)**

If $f(\sigma)$ is transformable, then its integral transform $F(z)$ exists in the region of convergence (ROC) ${\rm Re}(z^{2})>c$.

Proof

Divide the integral in equation (A.1) as

[TABLE]

where $0<a<b<\infty$.

The integral $\mathcal{I}_{2}$ converges, as the integrand is piecewise continuous in $[a,b]$ according to transformability condition (i) in Definition A.2. 2. 2.

According to transformability condition (ii), there exists $M_{1}\in\mathbb{R}_{>0}$ such that $f(\sigma)\leq M_{1}\exp(c/2\sigma^{2})$ for $\sigma\leq a$. Then using Jensen’s inequality, we have

[TABLE]

Therefore, the integral $\mathcal{I}_{1}$ converges in the region ${\rm Re}(z^{2})=x^{2}-y^{2}>c$, where $z=x+\mathrm{i}y$, $x\in\mathbb{R},y\in\mathbb{R}$. 3. 3.

According to transformability condition (iii), there exists $M_{2}\in\mathbb{R}_{>0}$ such that $f(\sigma)\leq M_{2}\sigma^{\lambda}$ for $\sigma\geq b$. Then, we have

[TABLE]

Here, we applied the inequality $\left\lvert\exp(-z^{2}/2\sigma^{2})\right\rvert\leq 1$ for ${\rm Re}(z^{2})>c\geq 0$. As a result, the integral $\mathcal{I}_{3}$ converges.

Therefore, the transform $F(z)$ exists in the ROC ${\rm Re}(z^{2})>c$.

Fig. 7 shows the ROC for the integral in equation (A.1). We can extend the ROC by the following two corollaries.

Corollary A.4

If a transformable function $f(\sigma)$ additionally satisfies the condition $f(\sigma)=\mathcal{O}(\sigma^{\beta}\exp(c/2\sigma^{2}))$ with $\beta\geq 1$ as $\sigma\to 0$, then the integral in equation (A.1) converges in the ROC ${\rm Re}(z^{2})\geq c$.

Proof

According to the additional condition, there exists $M_{3}\in\mathbb{R}_{>0}$ such that $f(\sigma)\leq M_{3}\sigma^{\beta}\exp(c/2\sigma^{2})$ for $\sigma\leq a$. Then, we can rewrite equation (A.3) as

[TABLE]

For ${\rm Re}(z^{2})=c$, this becomes

[TABLE]

Therefore, the ROC for the integral in equation (A.1) extends to ${\rm Re}(z^{2})\geq c$.

Corollary A.5

If a transformable function $f(\sigma)$ additionally satisfies $f(\sigma)=\mathcal{O}(\sigma^{\beta})$ with $\beta\geq 1$ as $\sigma\to 0$, then the integral in equation (A.1) converges in the ROC ${\rm Re}(z^{2})\geq 0$.

Proof

The proof is trivial with the substitution $c=0$ in Corollary A.4.

Next we prove the uniqueness theorem for the transform. First, we state a well-known proof of the following lemma for completeness.

Lemma A.6

If $f(x)$ is continuous in $[0,1]$, and $\int_{0}^{1}x^{n}f(x)\ \mathrm{d}x=0$ for $n=0,1,2,\dots$, then $f(x)=0$.

Proof

From the Weierstrass approximation theorem, for any $\epsilon>0$, there exists a polynomial $P_{\epsilon}(x)$ such that $\left\lvert f(x)-P_{\epsilon}(x)\right\rvert<\epsilon$ for all $x\in[0,1]$. The hypothesis implies that $\int_{0}^{1}P_{\epsilon}(x)f(x)\ \mathrm{d}x=0$. By taking the limit $\epsilon\to 0$, this equation becomes $\int_{0}^{1}f(x)f(x)\ \mathrm{d}x=0$. As $f(x)^{2}\geq 0$, we have $f(x)=0$.

Theorem A.7** (Uniqueness)**

If $f(\sigma)$ and $g(\sigma)$ are continuous, and $\mathcal{T}[f](z)=\mathcal{T}[g](z)$ for all $z$ in the ROC, then $f(\sigma)=g(\sigma)$.

Proof

Due to linearity, it is sufficient to prove that if $\mathcal{T}[f](z)=0$, then $f(\sigma)=0$. Take $d$ such that the contour ${\rm Re}(z^{2})=d$ lies in the ROC. By making the change of variables $s=\exp(-1/2\sigma^{2})$, for $z^{2}=d+n+1$ with $n=0,1,2,\dots$, we have

[TABLE]

This integral exists as $s\to 0$, because

[TABLE]

Therefore, according to Lemma A.6, we have $f(\sigma)=0$.

Theorem A.8** (Inversion)**

If $F(z)$ is the transform of $f(\sigma)$, then $f(\sigma)$ is given by the inverse transform

[TABLE]

where the contour $C$ is the hyperbola ${\rm Re}(z^{2})=d$ such that $C$ lies in ROC of $F(z)$.

Proof

Write equation (A.1) for $z^{2}=d$ as

[TABLE]

With the change of variables $p=1/\sigma^{2}$, this equation transforms into

[TABLE]

where

[TABLE]

Define a new function

[TABLE]

Theorem A.3 implies that $\int_{-\infty}^{\infty}\left\lvert h(p)\right\rvert\ \mathrm{d}p<\infty$, thus $h(p)$ belongs to the Lebesgue space $L^{1}(\mathbb{R})$. Therefore, we can take the Fourier transform of $h(p)$ as

[TABLE]

where we used equation (A.11) for the substitution in the last line. Now, take the inverse Fourier transform of $\hat{h}(\nu)$ as

[TABLE]

Here, we used the substitution of variable $z^{2}=d+2\mathrm{i}\nu$ in the last line, which transforms the integral path to the hyperbolic contour $C$ given by ${\rm Re}(z^{2})=d$.

Remark A.9

Equation (A.11) has the form of a Laplace transform. Therefore, the integral transform with a Gaussian kernel for a transformable function can be converted into a Laplace transform by suitable change of variables.

Appendix B Efficient algorithm to compute the Faddeeva function

In this appendix, we state an efficient algorithm to compute the Faddeeva function $w_{\rm F}(z)$ (Algorithm 1, for details see Zaghloul, 2017). The relative error of this algorithm is less than $4\times 10^{-5}$ over the whole complex plane.