Probing the small-scale structure in strongly lensed systems via transdimensional inference

TL;DR

This paper introduces a probabilistic cataloging approach for analyzing strongly lensed systems, enabling the detection and characterization of small-scale dark matter subhalos below traditional detection thresholds, thereby advancing our understanding of cosmic structure.

Contribution

It presents a novel transdimensional Bayesian framework for modeling lens systems that captures uncertainties and exploits all data information, improving subhalo mass function constraints.

Findings

Probabilistic cataloging robustly characterizes subhalo populations.

Method can detect subhalos below traditional detection thresholds.

Simulations show potential for future space telescopes to constrain dark matter structures.

Abstract

Strong lensing is a sensitive probe of the small-scale density fluctuations in the Universe. We implement a novel approach to modeling strongly lensed systems using probabilistic cataloging, which is a transdimensional, hierarchical, and Bayesian framework to sample from a metamodel (union of models with different dimensionality) consistent with observed photon count maps. Probabilistic cataloging allows us to robustly characterize modeling covariances within and across lens models with different numbers of subhalos. Unlike traditional cataloging of subhalos, it does not require model subhalos to improve the goodness of fit above the detection threshold. Instead, it allows the exploitation of all information contained in the photon count maps, for instance, when constraining the subhalo mass function. We further show that, by not including these small subhalos in the lens model, fixed-dimensional inference methods can significantly mismodel the data. Using a simulated Hubble Space Telescope (HST) dataset, we show that the subhalo mass function can be probed even when many subhalos in the sample catalogs are individually below the detection threshold and would be absent in a traditional catalog. With the planned Wide Field Infrared Space Telescope (WFIRST), simultaneous probabilistic cataloging of dark subhalos in high-resolution, deep strong lens images has the potential to constrain the subhalo mass function at even lower masses.

Tables

Table 1. Table 1: Parameter List of the PCAT Lens Metamodel

name	prior	min/mean	max/std	true	unit	explanation
${\mu }_{\mathrm{sub}}$	uniform	0	100	-		mean number of subhalos
$\beta$	Gaussian	1.9	0.5	1.9		slope of the deflection strength distribution
${\sigma }_{\text{psf}}$	Gaussian	0.087	0.01	0.087	arcsec	radius of the Airy disk that represents the PSF
$A$	$\log$-uniform	${10}^{-8}$	${10}^{-6}$	$2\times {10}^{-7}$	erg/cm2/s/Å/sr	normalization of the isotropic background emission
${\theta }_{1,\mathrm{src}}$	uniform	-2	2	RC	arcsec	horizontal coordinate of the background source
${\theta }_{2,\mathrm{src}}$	uniform	-2	2	RC	arcsec	vertical coordinate of the background source
$I_{\mathrm{src}}$	$\log$-uniform	${10}^{-20}$	${10}^{-15}$	${10}^{-18}$	erg/cm2/s/Å	flux of the background source
${\theta }_{e,\mathrm{src}}$	$\log$-uniform	0.1	2	0.5	arcsec	half-light radius of the background source
${ϵ}_{\mathrm{src}}$	uniform	0	0.3	R		ellipticity of the background source
${\varphi }_{\mathrm{src}}$	uniform	0	2$\pi$	R	radian	azimuthal orientation of the background source
${\theta }_{1,\mathrm{hst}}$	uniform	-2	2	RC	arcsec	horizontal coordinate of the host galaxy
${\theta }_{2,\mathrm{hst}}$	uniform	-2	2	RC	arcsec	vertical coordinate of the host galaxy
$I_{\mathrm{hst}}$	$\log$-uniform	${10}^{-20}$	${10}^{-15}$	${10}^{-16}$	erg/cm2/s/Å	flux of the host galaxy
${\theta }_{e,\mathrm{hst}}$	$\log$-uniform	0.1	2	1	arcsec	scale size of the Sérsic profile of the host galaxy
${\theta }_{E,\mathrm{hst}}$	$\log$-uniform	0.5	2	1.5	arcsec	Einstein radius of the host galaxy
${ϵ}_{\mathrm{hst}}$	uniform	0	0.5	R		ellipticity of the host galaxy
${\varphi }_{\mathrm{hst}}$	uniform	0	2$\pi$	R	radian	azimuthal angle of the ellipticity of the host galaxy
${\gamma }_{\mathrm{ext}}$	uniform	0	0.3	R		amplitude of the external shear
${\varphi }_{\mathrm{ext}}$	uniform	0	2$\pi$	R	radian	azimuthal angle of the external shear
$N_{\mathrm{sub}}$	Poisson	${\mu }_{\mathrm{sub}}$	$\sqrt{{\mu }_{\mathrm{sub}}}$	25		number of subhalos in a model
$\overrightarrow{{\theta }_1}$	uniform	-2	2	R	arcsec	horizontal coordinate of the $n$th subhalo
$\overrightarrow{{\theta }_2}$	uniform	-2	2	R	arcsec	vertical coordinate of the $n$th subhalo
$\overrightarrow{{\alpha }_s}$	power-law	0.01	1	R	arcsec	deflection strength of the $n$th subhalo
$\overrightarrow{{\theta }_s}$	uniform	0	0.1	R	arcsec	projected scale radius of the $n$th subhalo
$\overrightarrow{{\theta }_c}$	uniform	0	2	R	arcsec	projected cutoff radius of the $n$th subhalo

Equations

$${\tilde{f}}_{\mathrm{src}}({\theta }_1,{\theta }_2)=f_{\mathrm{src}}({\theta }_1-{\alpha }_{{\theta }_1},{\theta }_2-{\alpha }_{{\theta }_2}),$$

$$\overrightarrow{\stackrel{~}{\theta }}=\vec{\theta }-\vec{\alpha }(\vec{\theta }),$$

$$\vec{\alpha }=\frac{2}{c^2}{\vec{\nabla }}_{\theta }\int \text{differential}\chi \frac{{\chi }_{\mathrm{src}}-\chi }{\chi {\chi }_{\mathrm{src}}}\Psi (\vec{r}),$$

$$\vec{\alpha }=\frac{2}{c^2}\frac{D_{\mathrm{LS}}}{D_{\mathrm{S}}{D}_{\mathrm{L}}}{\vec{\nabla }}_{\theta }\int \text{differential}z\Psi (\vec{r}),$$

$$\vec{\alpha }={\vec{\alpha }}_{\mathrm{ext}}+{\vec{\alpha }}_{\mathrm{hst}}+\sum _{n=1}^N{\vec{\alpha }}_n.$$

$${\rho }_{\mathrm{hst}}(r)=\frac{{\sigma }_{\mathrm{hst}}^2}{2\pi G{r}^2},$$

$$\vec{\alpha}_{\rm{hst}}=\frac{\theta_{\rm{E,hst}}q_{\rm{hst}}}{\sqrt{1-q_{\rm{hst}}^{2}}}\Bigg{(}\tan^{-1}\Big{(}\frac{\sqrt{1-q_{\rm{hst}}^{2}}\theta_{1,\rm{hst}}^{\prime}}{\omega}\Big{)}\hat{\theta}_{1,\rm{hst}}^{\prime}+\\ \tanh^{-1}\Big{(}\frac{\sqrt{1-q_{\rm{hst}}^{2}}\theta_{2,\rm{hst}}^{\prime}}{\omega}\Big{)}\hat{\theta}_{2,\rm{hst}}^{\prime}\Bigg{)},$$

$${\rho }_n(r)=\frac{M_{0,n}}{4\pi r_{\mathrm{s},\text{n}}^3}\frac{1}{(r/r_{\mathrm{s},\text{n}})(1+r/r_{\mathrm{s},\text{n}}{)}^2}\frac{{\tau }_n^2}{(r/r_{\mathrm{s},\text{n}}{)}^2+{\tau }_n^2},$$

$$M_{\mathrm{c},\text{n}}=M_{0,n}\frac{{\tau }_n^2}{({\tau }_n^2+1{)}^2}\left[({\tau }_n^2-1)\ln {\tau }_n+\pi {\tau }_n-(1+{\tau }_n^2)\right].$$

$${\vec{\alpha }}_{\mathrm{ext}}={\gamma }_{\mathrm{ext}}\left(\begin{array}{cc}\cos 2{\varphi }_{\mathrm{ext}}& +\sin 2{\varphi }_{\mathrm{ext}}\\ \sin 2{\varphi }_{\mathrm{ext}}& -\cos 2{\varphi }_{\mathrm{ext}}\end{array}\right)\vec{\theta }.$$

$$\begin{array}{rl}(1-{\kappa }_{\mathrm{ext}})& \to \lambda {\kappa }_{\mathrm{ext}},\\ {\gamma }_{\mathrm{ext}}& \to \lambda {\gamma }_{\mathrm{ext}}.\end{array}$$

$$\displaystyle f_{\rm{gal}}(\theta)=f_{\rm{e,gal}}\exp\Bigg(-b_{n}\Big{(}(\theta/\theta_{\rm{e,gal}})^{1/n}-1\Big{)}\Bigg{missing}),$$

$$b_n=2n-\frac{1}{3}+\frac{4}{405n}+\frac{46}{25515n^2}$$

$$f_{\text{m}}=\mathcal{P}\ast (f_{\mathrm{bac}}+f_{\mathrm{hst}}+{\tilde{f}}_{\mathrm{src}}).$$

$$\begin{array}{rl}\ln P(\mathcal{D}|\mathcal{M})& =\sum _{ij}\ln P(k_{ij}^{\mathrm{d}}|k_{ij}^{\mathrm{m}})\\ & =\sum _{ij}{k}_{ij}^{\mathrm{d}}\ln k_{ij}^{\mathrm{m}}-k_{ij}^{\mathrm{m}}-\ln k_{ij}^{\mathrm{d}}!,\end{array}$$

$$\begin{gathered} k_{ij}^{\mathrm{m}}= \\ \iiint f_{\text{m}}(E,{\theta }_1,{\theta }_2)ϵ(E,{\theta }_1,{\theta }_2)T(E)\text{differential}E\text{differential}{\theta }_1\text{differential}{\theta }_2, \end{gathered}$$

$${\tilde{k}}_m({\theta }_1,{\theta }_2)={\bar{k}}_m({\theta }_1,{\theta }_2)+{\vec{\alpha }}_n\cdot \vec{\nabla }{\bar{k}}_m({\theta }_1,{\theta }_2)+\mathcal{O}({\vec{\alpha }}_n^2).$$

$$f_{\mathrm{E},\mathrm{sub}}\equiv \frac{{\sum }_n{\int }_{\mathrm{ann}}{\kappa }_n{\text{differential}}^2\theta }{{\int }_{\mathrm{ann}}{\kappa }_{\mathrm{hst}}{\text{differential}}^2\theta }.$$

$$\displaystyle\begin{split}\alpha_{n}&\equiv\theta\frac{\bar{\Sigma}_{n}(\theta)}{\Sigma_{\rm{crit}}}=\frac{M_{0,n}}{\pi D_{\rm{L}}^{2}\Sigma_{\rm{crit}}\theta_{{\rm s},n}}\frac{1}{\theta^{\prime}}\Big{(}\ln\frac{\theta^{\prime}}{2}+F(\theta^{\prime})\Big{)}\\ &=\frac{\alpha_{{\rm s},n}}{\theta^{\prime}}\Big{(}\ln\frac{\theta^{\prime}}{2}+F(\theta^{\prime})\Big{)},\end{split}$$

$${\Sigma }_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{S}}}{D_{\mathrm{LS}}{D}_{\mathrm{L}}}$$

$$F({\theta }^{\prime })=\{\begin{array}{ll}\frac{\arccos (1/{\theta }^{\prime })}{\sqrt{{{\theta }^{\prime }}^2-1}}& \text{if θ′>1}\\ 1& \text{if θ′=1}\\ \frac{\text{arccosh}(1/{\theta }^{\prime })}{\sqrt{1-{{\theta }^{\prime }}^2}}& \text{if θ′<1.}\end{array}$$

$${\alpha }_n$$

$$\displaystyle L(\theta^{\prime})=\ln\Bigg(\frac{\theta^{\prime}}{\tau_{n}+\sqrt{\tau_{n}^{2}+{\theta^{\prime}}^{2}}}\Bigg{missing}).$$

Full text

gravitational lensing: strong, methods: data analysis, methods: statistical, cosmology: dark matter

††slugcomment: Submitted to ApJ††software: Probabilistic Cataloger (PCAT) (Daylan et al., 2017)

1 Introduction

The concordance $\Lambda$ Cold Dark Matter ($\Lambda$CDM) model of cosmology is currently the most accurate description of matter density fluctuations in the Universe (Davis et al., 1981; Blumenthal et al., 1982, 1984; Davis et al., 1985). Through the collapse and bottom-up hierarchical merger of dark matter halos, it explains how structure could form in less than a billion years after the Big Bang. At the largest scales probed by the Cosmic Microwave Background (CMB), the six-parameter model agrees with precise observations of CMB anisotropies (Ade et al., 2016). Nevertheless, discrepancies between observations and $\Lambda$CDM predictions have been claimed to remain at subgalactic scales (see, e.g. Bull et al., 2016; Del Popolo & Delliou, 2016). For instance, the number of observed satellites in the Local Group ($\sim$ 50) falls short of the much larger numbers of low-mass subhalos (> 1000) predicted by $N$-body simulations of $\Lambda$CDM (Klypin et al., 1999; Moore et al., 1999). Furthermore, the projected density profiles of the observed subhalos are rather shallow in the central regions, compared to the much cuspier profiles predicted by $\Lambda$CDM (see, e.g. Flores & Primack, 1994; Moore, 1994; Burkert, 1995; Gentile et al., 2004; Walker & Penarrubia, 2011; Oh et al., 2015).

These apparent discrepancies can be resolved in several ways. First, the inclusion of baryonic feedback processes within numerical simulations (Brooks et al., 2013; Sawala et al., 2016; Wetzel et al., 2016; Despali & Vegetti, 2017; Sawala et al., 2017; Creasey et al., 2017; Garrison-Kimmel et al., 2017; Strigari et al., 2007; Nierenberg et al., 2016; Amorisco et al., 2014) has been shown to significantly reduce the tension between predictions and observations. Alternatively, the discrepancy could be due to a suppression of the small-scale matter power spectrum prior to structure formation. This may arise from the free streaming of dark matter (DM) particles after they kinetically decouple from the rest of the cosmic plasma, which partially erases inhomogeneities at the $\sim$ kpc comoving scale (Bond & Szalay, 1983; Dalcanton & Hogan, 2001; Bode et al., 2001; Boyanovsky et al., 2008; Boyanovsky & Wu, 2011). Another venue to resolve the tension is to modify the microphysics of DM by, for instance, postulating that DM particles interact with each other (Spergel & Steinhardt, 2000; Yoshida et al., 2000; Davé et al., 2001; Colín et al., 2002), or with a relativistic species (see, e.g. Carlson et al., 1992; Boehm et al., 2002; Ackerman et al., 2009; Feng et al., 2009; Kaplan et al., 2010; van den Aarssen et al., 2012; Cyr-Racine & Sigurdson, 2013; Chu & Dasgupta, 2014). Both possibilities result in a modification to the internal density profile of DM halos, while the latter could also suppress the number of small-mass subhalos (see, e.g. Buckley et al., 2014; Schewtschenko et al., 2015; Vogelsberger et al., 2016). Last, the problem could indicate a suppression of the primordial power spectrum following inflation (Kamionkowski & Liddle, 2000; Minor & Kaplinghat, 2015). In any case, the possibility of resolving the potential small-scale crisis of $\Lambda$CDM and exploring the particle nature of dark matter makes a strong case for observationally probing the subgalactic structure of DM.

Small-scale DM structure can, in principle, be probed by studying light emission from faint and ultra-faint dwarf galaxies in the Local Group and beyond. However, star formation in these low-mass subhalos is suppressed due to reionization, interaction with their host, or self-quenching (e.g. Fitts et al., 2016; Read et al., 2017). Consequently most of the subhalos with virial halo mass $\lesssim 10^{8}M_{\odot}$ should be mostly dark, making their gravitational lensing (Einstein, 1915) the only observable probe of their properties. Gravitational lensing thus provides a unique window into the smallest self-bound DM structure populating our Universe. In particular, galaxy-scale strong lensing systems in which a background source (e.g. galaxy or quasar) is multiply imaged by a foreground massive galaxy provide a key way to probe the DM distribution on subgalactic scales.

Indeed, the observation of flux-ratio anomalies in strongly lensed quasars (Mao & Schneider, 1998; Metcalf & Madau, 2001; Chiba, 2002; Metcalf & Zhao, 2002; Dalal & Kochanek, 2002; Keeton et al., 2003; Metcalf & Amara, 2012; Fadely & Keeton, 2011; Nierenberg et al., 2014; gilman_strong_2017; Nierenberg et al., 2017, 2017; Gilman et al., 2017; Hsueh et al., 2017) has been used to study the abundance of substructure within lens galaxies. In particular, Dalal & Kochanek (2002) used the observation of seven radio-loud lensed quasars to put a constraint on the substructure mass fraction of $0.006<f_{\rm sub}<0.07$ (90% confidence level). More recently, direct gravitational imaging of mass substructure (Koopmans, 2005; Vegetti & Koopmans, 2009; Hezaveh et al., 2013) lying close in projection to magnified arcs and Einstein rings has opened a new avenue to study the small-scale distribution of DM. This method works by fitting observed strong lenses with a smooth mass distribution, i.e. without substructure, and then adding pixel-based corrections Vegetti et al. (2010, 2012) or a number of linearized subhalos Hezaveh et al. (2016) to the mass model, with the requirement that the additional degrees of freedom improve the goodness of fit above some stringent detection threshold, e.g., greater than $\sim 10\sigma$. For instance, Vegetti et al. (2010, 2012) detected substructure with masses above $10^{8}M_{\odot}$ at very high statistical significance using optical images of two strong lens systems. Similarly, a $10^{9}M_{\odot}$ substructure was detected by Hezaveh et al. (2016) in an interferometric data set at a statistical significance of $6.9\sigma$ in the strong lens SDP.81 (ALMA Partnership et al., 2015). Constraints on the mass fraction in substructure and on the subhalo mass function can, in principle, be placed using these detections along with nondetections in other lens systems (e.g. Vegetti & Koopmans, 2009; Vegetti et al., 2014; Li et al., 2016; Hezaveh et al., 2016; Li et al., 2017).

In practice, turning a direct mass substructure detection into constraints on the subhalo mass function and their internal properties (e.g. concentration and truncation) is highly nontrivial. Indeed, degeneracies in the lensing model can blur the interpretation of a direct substructure detection (e.g. Minor et al., 2016). For instance, covariances between the smooth lens model, the structure of the source, and the subhalo population can lead to biases and underestimated error bars in the recovered substructure properties if they are not explored properly. Further, the nonlinearity of the lensing magnification implies that the likelihood improvement resulting from the addition of a subhalo to the lens mass model is not entirely independent of the presence of other subhalos within the lens. Although less likely, it remains a concern that two less massive subhalos serendipitously appearing in projection along the same line-of-sight could be misinterpreted as one more massive subhalo. Properly taking into account these degeneracies is important in order to obtain fully realistic error bars on subhalo population parameters (later identified as hyperparameters).

With the exception of several studies (e.g. Birrer et al., 2015; Hezaveh et al., 2016; Cyr-Racine et al., 2016; Birrer et al., 2017), which aim to infer subhalos below the detection threshold statistically, current analyses (e.g. Vegetti et al., 2014; Li et al., 2016; Hezaveh et al., 2016) place constraints on the subhalo population a posteriori using highly statistically significant detections (or nondetections) of mass substructure within lens galaxies. In general, this requires writing down an approximate likelihood for the abundance of detected substructure, given some subhalo population parameters characterizing their mass function and spatial distribution. We identify two issues with this approach. First, by exclusively using an extremely reduced version of the data (i.e. the position and mass of directly detected substructure), these analyses cannot take into account possible lensing degeneracies in which, for instance, the gravitationally imaged mass substructure actually corresponds to two close-by physical subhalos. Another way to phrase this problem is to say that correlations between models having different numbers of subhalos are neglected, which could potentially bias the recovered subhalo mass function constraints. Second, the use of highly reduced data neglects, by construction, marginal subhalo detections that, while not constituting statistically significant direct detections of substructure, can still provide statistical constraints on the subhalo population hyperparameters.

In this work, we develop a Bayesian data analysis framework that directly uses the lensed images as input and addresses the above issues by performing transdimensional sampling over the union of lens models with different numbers of physical subhalos (hereafter, the metamodel). Specifically, the lens metamodel refers to the union of mass models with a common macro (or smooth) lens model, and different numbers of subhalos, each with their own properties. The approach applies probabilistic cataloging (Hogg & Lang, 2010; Brewer et al., 2013; Daylan et al., 2017; Brewer et al., 2015; Portillo et al., 2017) to images of strongly lensed systems. In contrast with standard gravitational imaging, our technique outputs an ensemble of probability-weighted subhalo catalogs providing good fits to the strongly lensed images. We use a version of the recently developed Markov chain Monte Carlo (MCMC) sampler, . PCAT is a Poisson mixture sampler that provides a transdimensional, hierarchical, and Bayesian inference framework to sample from the posterior probability distribution of a metamodel consistent with Poisson-distributed count data. In the context of the lensed images, we assume that the count data is collected by a photon-counting experiment such as a charged coupled device (CCD), but our results could, in principle, incorporate other types of auxiliary data (e.g., time delays of multiple images when working with quasar background sources (Suyu et al., 2017), or spectroscopic data of the foreground galaxy to constrain its mass budget (Mason et al., 2017)), as well as informative priors on the relevant metamodel parameters.

Our approach shares its transdimensional nature with that of (Brewer et al., 2015), which implemented an independent transdimensional framework to infer substructure in strong lens systems. However, in addition to differences concerning substructure modeling and the sampling method used, our code and analysis were developed independently of theirs, and hence provide an important cross-check of the method and results. It has become customary in cosmological analyses to have at least two publicly available independent codes to check for numerical accuracy and find bugs. While both codes are public, we understand that this approach is not yet ubiquitous in the field of galaxy-scale strong lensing (for instance, it has been difficult thus far to reproduce/check the results of Vegetti et al. (2014) and Hezaveh et al. (2016)), but we hope that our work can help make substructure lensing more reproducible. Given the uniqueness of the constraints it can provide on low-mass substructures and the impact that these detections can have on dark matter physics, we believe it is extremely important to have multiple cross-checks to validate any results. On a more technical level, our approach differs from that of (Brewer et al., 2015) in two major ways:

• •

Subhalo lens modeling. We represent subhalos with NFW profiles with a certain mass, scale, and cutoff radius, whereas (Brewer et al., 2015) uses the so-called blobs, i.e., deflection as a power-law in the subhalo-centric radius. We think that choosing a physically motivated basis function set to represent the deflection field due to subhalos is essential for correctly marginalizing them out, i.e., propagating their (rather large) uncertainties to the parameter of interests such as the subhalo mass function.

• •

The sampling method and the sampler used to sample from the posterior probability distribution of the (different) lens models given the data. (Brewer et al., 2015) uses Diffusive Nested Sampling (DNS) as implemented in RJObject. We use a homebrew reversible-jump MCMC sampler with Metropolis-type within-model proposals.

Our work is also significantly different than that of (Fadely & Keeton, 2011). In their analysis of the lens HE0435, they either consider models with a fixed number (1 to 3) of subhalos, thereby not taking into account the covariance between models with different numbers of subhalos, or they consider large populations of subhalos that are Poisson-distributed across the lens plane. Our approach of letting the data drive the actual number of subhalos is a significant improvement over this previous work.

After thoroughly sampling from the lens metamodel, one can marginalize over the subhalo catalogs to infer the underlying lens mass model and population characteristics of the subhalos, such as their mass function. The latter is achieved via hierarchical modeling, where the probability distributions of the subhalo properties, e.g., the subhalo mass function, are parameterized by hyperparameters. When the subhalos in the metamodel are marginalized out, the posterior distribution of the hyperparameters constrains the population characteristics of the subhalos, even though most of the sampled subhalos are below the detection threshold.

We base our inference directly on the observed photon count maps, as opposed to secondary, potentially biased estimates, such as astrometric perturbations or relative magnifications of the multiple images of the background source. This allows one to exploit all the information contained in the photon count maps regarding our metamodel, and to fully propagate uncertainties in the observed data to the parameter of interest, e.g., the subhalo mass function.

To illustrate the features of probabilistic cataloging, we focus here on simple smooth lens models and source configurations, and present results for mock data sets; doing so allows one to study potential biases and systematics. Of course, analyses of real data would require more complex lens models – and especially a more advanced lensing source reconstruction (e.g. Tagore & Keeton, 2014). While such complications are likely to make sampling from the metamodel more numerically onerous, they will not change the details of the probabilistic cataloging approach described below.

The rest of the paper is structured as follows. In Section 2, we discuss how we model lensed images. We then give an overview of probabilistic cataloging in Section 3. We present our results using a mock data set in Section 4 and assess the inference performance of probabilistic cataloging under mismodeling in Sections 5 and 6. We discuss our results and conclude in Section 7.

2 Lens metamodel

In this section, we outline our lens metamodel and how we forward-model a strong lens image. The lensed image, $\tilde{f}_{\rm{src}}$, can be written as the source brightness profile $f_{\rm{src}}$, evaluated on the image plane,

[TABLE]

where $\vec{\alpha}=(\alpha_{\theta_{1}},\alpha_{\theta_{2}})$ is the deflection vector at the image position, $\vec{\theta}=(\theta_{1},\theta_{2})$. Here, $\theta_{1}$ and $\theta_{2}$ denote the horizontal and vertical coordinates on the two-dimensional image plane, and $\vec{\tilde{\theta}}$ is the position on the source plane, such that

[TABLE]

which is the lens equation. Furthermore, the deflection field $\vec{\alpha}$ can be written as the gradient of the Newtonian gravitational potential integrated along the line-of-sight,

[TABLE]

where $\chi$ is the comoving distance along the line-of-sight, $\chi_{\rm{src}}$ is the comoving distance to the source, and $\Psi(\vec{r})$ is the gravitational potential at the three-dimensional position vector, $\vec{r}$.

As in most lensing works, we make use of the thin lens approximation, whereby we assume that the thickness of the lensing halo is much smaller than other distance scales in the problem. This reduces Equation (3) to

[TABLE]

where $z$ is the line-of-sight coordinate, $D_{\rm{L}}$ and $D_{\rm{S}}$ are the angular diameter distances from the observer to the lens and source planes, respectively, and $D_{\rm{LS}}$ is the angular diameter distance from the lens to the source plane. Analytic expressions for the line-of-sight integral of various foreground mass components can then be used to calculate the deflection field (see, e.g. Keeton, 2001).

2.1 Deflection field

Because of the linearity of the Poisson equation, the deflection field, $\vec{\alpha}$, can be expressed as the superposition of deflections due to several components. In a typical strong lens system, the dominant contribution comes from the smooth mass distribution of the main foreground galaxy (hereafter called the host halo), $\vec{\alpha}_{\rm{hst}}$ and the external shear due to the large-scale structure (LSS), $\vec{\alpha}_{\rm{ext}}$. Subhalos perturb this deflection field at the percent level, with individual contributions to the deflection field, $\vec{\alpha}_{n}$, where the subscript $n=1,2,..,N$ is the index of a subhalo. As a result, the total deflection field is given by

[TABLE]

Let us now describe each component separately.

Smooth halo

Cosmological $N$-body simulations of $\Lambda$CDM (Genel et al., 2014) suggest that the equilibrium distribution of collisionless dark matter particles is approximately spherically symmetric about the dynamical center of mass of the self-gravitating halo of dark matter particles. Furthermore, the resulting virialized halos have a universal radial profile that is well-described by a broken power-law in the radial halo-centric distance, $r$, and is known as the Navarro Frenk White (NFW) profile (Dubinski & Carlberg, 1991; Navarro et al., 1996). However, the host halo also possesses baryonic matter, which actually dominates the mass budget in the bulge. Taking the mass budget of the host halo to be the sum of the baryonic and DM components, it has been found that the overall mass density profile of early-type galaxies (ETGs) can be fitted well by singular isothermal ellipsoids (SIEs), which is known as the bulge-halo conspiracy (Treu, 2010) because neither dark nor baryonic matter density distribution is individually isothermal.

We therefore take the mass density profile of the host halo as

[TABLE]

where $\sigma_{\rm{hst}}$ is the (constant) velocity dispersion and $G$ is the gravitational constant. In general, galaxies also possess some ellipticity due to the existence of a baryonic disc or as artifacts of recent gravitational interactions. By taking into account the ellipticity of the isothermal gas, and assuming that the contours of equal surface mass density are in the form of concentric ellipses, the deflection is given by (Keeton, 2001; Kassiola & Kovner, 1993),

[TABLE]

where $\omega\equiv\sqrt{q_{\rm{hst}}^{2}{\theta_{1}^{\prime}}^{2}+{\theta_{2}^{\prime}}^{2}}$, $q_{\rm{hst}}\equiv 1-\epsilon_{\rm{hst}}$ is the minor to major axis ratio, and $\epsilon_{\rm{hst}}$ is the ellipticity. Note that the primed quantities, i.e, $\theta_{1,\rm{hst}}^{\prime}$, $\theta_{2,\rm{hst}}^{\prime}$, $\hat{\theta}_{1,\rm{hst}}^{\prime}$, and $\hat{\theta}_{2,\rm{hst}}^{\prime}$, denote the angular coordinates and the associated unit vectors along the minor and major axes of the elliptical host halo, respectively. If we use $\phi_{\rm hst}$ to denote the angle between the major axis of the ellipse and the horizontal axis of the image, $\theta_{1}$, then the primed coordinates are obtained by translating the image by $(\theta_{1,\rm{hst}},\theta_{2,\rm{hst}}$) and then rotating it by $\phi_{\rm{hst}}$.

Furthermore, the lensing strength of the isothermal deflector is parameterized using the projected Einstein radius, $\theta_{E,\rm{hst}}$, i.e., the radius of the circular image that would be produced if the lens were spherical and exactly aligned with the background light source. In a more general sense that holds for non-spherical lenses, we define $\theta_{E,\rm{hst}}$ as the radius of a circle that has an area equal to the area inside the critical curve.

Note that many observed strong lenses require departures from an isothermal model, where the preferred inner log-slope of the three-dimensional density distribution can be significantly different from -2. However, because the method we present is independent of the details of the lens model, we assume a relatively simple mass model for the smooth halo, leaving further refinements to future work.

Subhalos

$\Lambda$CDM predicts that, with cosmological time, more density peaks collapse to form new self-gravitating halos, while halos that formed earlier are accreted into more massive ones. Because hierarchical structure formation is a continuous and ongoing process, at any given time some fraction of matter is expected to be tied in halos that have been accreted into more massive halos. In the rest of the paper, we refer to these dark, gravitationally bound satellites inside the virial extent of the host as subhalos.

Subhalos grow by accretion and lose mass via tidal stripping, collisions, and evaporation. The relative efficiency of these processes determines the level of substructure in a host halo of mass $M_{\rm{hst}}$ and redshift $z_{\rm{hst}}$. $N$-body simulations can predict the abundance and properties of these subhalos, although the finite mass resolution of simulations can preclude robust conclusions (van den Bosch, 2017). Nevertheless, it is expected that subhalos have a truncated mass profile due to tidal stripping by the host galaxy. Therefore, we take the three-dimensional subhalo density profile of the $n$th subhalo, $\rho_{n}(r)$, to be the NFW profile with a power-law truncation as (Baltz et al., 2009):

[TABLE]

where $M_{\rm{c},\textit{n}}$ is the truncated mass of the same subhalo. In this parameterization most of the mass of a subhalo is contained inside the scale radius, $r_{\rm{s},\textit{n}}$, which would be the subhalo-centric distance at which the density profile steepens from a power-law with index -1 to a power-law with index -3 if it was not tidally truncated. Moreover, the central mass concentration is encoded in $\tau_{n}=r_{\rm{c},\textit{n}}/r_{\rm{s},\textit{n}}$ similarly to the concentration parameter $c=r_{\rm{vir},\textit{n}}/r_{\rm{s},\textit{n}}$ in the non-truncated case, where $r_{\rm{vir},\textit{n}}$, $r_{\rm{c},\textit{n}}$ and $r_{\rm{s},\textit{n}}$ are the virial, tidal cutoff, and scale radii of a subhalo. $M_{0,n}$ is again a mass scale relevant to the $n$th subhalo. Subhalos with large $\tau_{n}$ represent rather old, highly concentrated subhalos that orbit the host at large halocentric distances. In contrast, those with low $\tau_{n}$ typically correspond to subhalos that are close to the center and are thus being tidally stripped by the host. Note that we use the subscript $n$ when referring to a parameter of the $n$th subhalo when the host halo has the same type of parameter.

It is also expected that subhalos near the dynamical center of the host halo, and thus closest to disruption, should be tidally elongated. Nevertheless, because current data sets do not have the signal-to-noise ratio (SNR) necessary to constrain departures from spherical symmetry, we neglect any triaxiality in the subhalo shapes and assume that they are spherically symmetric. Furthermore, because virialized halos are generally found to have concentrations that reflect the average background density of the Universe at the time of their collapse, earlier halos are expected to have higher concentrations. Upon accretion into the host halo, their high concentrations can help them remain intact despite strong tidal stretching and evaporation by the host halo. However, motivated by recent $N$-body simulations (Springel et al., 2008), we neglect substructure inside subhalos and assume that the subhalos of the host galaxy have smooth mass distributions.

Although the mass of a non-truncated NFW profile is logarithmically divergent with halocentric radius, tidal truncation makes the total mass $M_{\rm{c},\textit{n}}$ finite (Baltz et al., 2009) such that

[TABLE]

This finite mass is plotted, in Figure 1 for a subhalo in units of $M_{0,n}$. It changes by about an order of magnitude in the relevant region of $1\lesssim\tau_{n}\lesssim 10$. Below $\tau_{n}\sim 1$, the subhalo is likely to be completely disrupted by evaporation and tidal stripping. We report truncated masses when referring to subhalos.

Our simulated images are approximately centered at the dynamical center of the host halo and have a size comparable to its Einstein radius, which is much smaller than its virial radius. Because we work on the (projected) lens plane, subhalos in the field of view have a large spread in their three-dimensional distance from the dynamical center of the host. As a result, there should be a large spread in their deflection profile. We crudely take this into account by drawing the projected scale and cutoff radii of the mock subhalos from a uniform distribution. The limits of these distributions are given in Table 1. The resulting ratios between the projected cutoff and scale radii, $\tau_{n}$, fall between 3 and 25.

We emphasize that these mock subhalos are meant to provide a test bed for introducing probabilistic cataloging of subhalos and do not fully reflect their rich dynamics. Furthermore, although simulations and analytical treatments of $\Lambda$CDM can predict the scale and cutoff radii of a given subhalo, if provided its mass and distance from the center of the host galaxy (Metcalf & Madau, 2001), which implies that predictions can be folded into the inference as a joint prior distribution on the projected scale and cutoff radii, the fact that strong lens images are not informative about the three-dimensional distance of a subhalo from the center of its host still means there is a degeneracy in this description.

We parameterize the deflection field due to the $n$th subhalo, $\alpha_{n}$, using the angular distances corresponding to the scale and cutoff radii, $\theta_{\rm{s},\textit{n}}=r_{\rm{s},\textit{n}}/D_{\rm{L}}$ and $\theta_{\rm{c},\textit{n}}=r_{\rm{c},\textit{n}}/D_{\rm{L}}$, and the normalization of the deflection profile, $\alpha_{\rm{s},\textit{n}}$, which we refer to as the deflection strength. The functional form of this deflection profile is presented in Appendix A.

The deflection as a function of radial distance in three dimensions from the dynamical center of mass of the subhalo is proportional to the integrated mass up to that radius. Therefore, it initially rises because of the shallow log-slope of the NFW profile, turns over at the scale radius, and further steepens beyond the cutoff radius. The projected deflection profile of a subhalo is given in Figure 2.

$\Lambda$CDM predicts that the variance of the matter density fluctuations in the Universe monotonically increases as one averages the density field inside successively smaller spheres. Therefore, the number of collapsed halos (and hence of subhalos) increases steeply as a power-law toward low masses with a log-slope of $\sim-1.9$ (Springel et al., 2008). When $\Lambda$CDM is modified at small scales via self-interactions, free-streaming at kinetic decoupling, or by the Heisenberg uncertainty in the case of wave-like dark matter (Hu et al., 2000), this is no longer the case and one expects a mass scale below which the subhalo mass function is suppressed.

Because subhalo deflection strengths are proportional to mass, we assume $\alpha_{\rm{s},\textit{n}}$ are also distributed as a power-law between some minimum and maximum value. However, note that the lower limit of this prior does not necessarily correspond to a physical cutoff in the mass distribution. Rather, it is a computational convenience, because probabilistic cataloging of very low-mass subhalos results in very little information gain – at the expense of significantly increased computational cost. Furthermore, the pixel size and the point-spread function (PSF) of the photon-collecting device sets an angular distance scale below which deflections cannot significantly affect the image. This scale also depends on the level of background emission and the flux of the lensed light source. For our simulated images of the Wide Field Camera 3 / Ultraviolet Visible (WFC3/UVIS) camera on board the Hubble Space Telescope (HST), whose pixel size is 0.04 arcsec, this scale roughly corresponds to $\sim$ 0.01 arcsec. We therefore choose the minimum allowed value of $\alpha_{\rm{s},\textit{n}}$ as 0.01 arcsec, in order to allow low-significance subhalos into the lens metamodel. We also choose the maximum of $\alpha_{\rm{s},\textit{n}}$ to be 1 arcsec. However, this maximum value does not affect our prior on $\alpha_{\rm{s},\textit{n}}$ because most of the prior volume is contained at small $\alpha_{\rm{s},\textit{n}}$ due to the steepness of the power-law. The normalization of this power-law is then given by $\mu_{\rm{sub}}$.

In the hierarchical modeling approach of probabilistic cataloging, the negative of the log-slope of the $\alpha_{\rm{s},\textit{n}}$ distribution, $\beta$, is also allowed to vary. In this inference scheme, the sampler visits all prior configurations consistent with the data and generates a posterior distribution over $\beta$. Therefore, although $\beta$ parameterizes the prior, it is treated just like any other parameter in the model, and is referred to as a hyperparameter. It is the prior on this hyperparameter that specifies our prior belief on the distribution of $\alpha_{\rm{s},\textit{n}}$. We set this prior such that $\beta$ itself is Gaussian-distributed with a mean and standard deviation of 1.9 and 0.5, respectively. The mean reflects the value motivated by $\Lambda$CDM, whereas the spread accounts for any deviations and allows the sampler to visit subhalo configurations that deviate from a power-law with index 1.9.

The deflection strength, $\alpha_{\rm{s},\textit{n}}$, and projected scale and cutoff radii, $\theta_{\rm{s},\textit{n}}$ and $\theta_{\rm{c},\textit{n}}$, of a subhalo uniquely determine its truncated mass, $M_{\rm{c},\textit{n}}$. Therefore, given the above priors, there is no hard boundary on the prior distribution of the subhalo truncated masses. However, most of the prior volume falls in the interval $10^{8}-10^{9}M_{\odot}$ for the fiducial host and source redshifts of $z_{\rm{hst}}=0.2$ and $z_{\rm{src}}=1$. This interval changes as a function of angular diameter distances to the host and source planes. To put our prior mass interval into cosmological context, Figure 3 shows the lower bound of the prior on the truncated mass of a subhalo, if it had $\theta_{\rm{s},\textit{n}}=0.05$ arcsec and $\theta_{\rm{c},\textit{n}}=1$ arcsec.

Finally, we adopt a spatially uniform prior on the projected surface density of subhalos. The prior structure of our lens metamodel is summarized in Table 1.

External shear

Lensing galaxies are usually found in galaxy groups or clusters whose overall contribution adds angular structure to the deflection field in the vicinity of the lensing galaxy. Even in the case of an isolated galaxy, matter distribution in the foreground of the lensing galaxy can introduce additional multipole terms to the deflection field. We parameterize this external deflection by means of a reduced shear field $\gamma_{\rm{ext}}$ oriented at an angle $\phi_{\rm{ext}}$ with respect to the longitudinal axis (Keeton et al., 1996). The field is traceless and invariant under a rotation of 180∘.

[TABLE]

Note that we only model the reduced shear, $\gamma_{\rm{ext}}=\gamma_{\rm{ext}}^{\prime}/(1-\kappa_{\rm{ext}})$, where $\kappa_{\rm{ext}}$ represents a potentially underlying, spatially uniform convergence in the field of interest. This is because adding a constant mass sheet to the model, $(1-\kappa_{\rm{ext}})$, reproduces the same lensed image by simultaneously rescaling the source position and flux or the deflector mass (Falco et al., 1985). These degenerate mass models are linked by the transformations

[TABLE]

An independent measurement of the deflector mass, e.g., using stellar kinematics, could constrain $\lambda$ (Schneider & Sluse, 2013; Bradač et al., 2004). Alternatively, because the magnification field and light travel-time change with $\lambda$, knowledge of the non-lensed fluxes of the background sources or having multiple images of the same system taken at different times could break this degeneracy. Specifically, the knowledge of time delays between multiple images of a background quasar and the time-delay distance allows one to constrain the gravitational potential of the main deflector, independent of the astrometric information that we use in this work. Because we only rely on a single photometric exposure of the strong lens system, we simply determine all convergences up to the transformation in Equation 12 and do not attempt to lift the degeneracy.

We note that the inclusion of only a quadrupolar shear field is motivated by the assumption that the lensing galaxy is isolated and gravitationally relaxed. If there are nearby galaxies, higher-order multipoles would be needed.

2.2 Light emission from the source

In general, the background light source in a strong lens system could be: a quasar; a blue, young galaxy with bright Lyman-$\alpha$ emission; or a red ETG similar to the lensing galaxy we consider in this work. When the background light source is a quasar, i.e., a point source without any spatial feature, the resulting multiple images are not extended and can be used to probe subhalos only over a small region on the lens plane depending on the mass of the subhalo. We therefore instead consider galaxy-galaxy lenses, where light from a background galaxy is strongly deflected by a giant elliptical in the foreground. In these systems, multiple images are more extended and a higher fraction of the pixels are informative on the subhalo parameters.

We assume that the emission from both the foreground and the background galaxies have Sérsic profiles, i.e., the surface brightness is given by (Sérsic, 1963)

[TABLE]

where the subscript gal refers to both lensing (foreground) and lensed (background) galaxies, i.e., gal={hst,src}. In this relation, $b_{n}$ is a coefficient that depends on the index $n$ (Ciotti & Bertin, 1999),

[TABLE]

which controls the level at which the inner and outer slopes anti-correlate. Furthermore, $\theta_{\rm{e,gal}}$ is the projected distance within which half of the total emission is contained and $f_{\rm{e,gal}}$ is the surface brightness at this radius. We assume $n=4$, which corresponds to the de Vaucouleurs profile (de Vaucouleurs, 1948). For this profile, the surface brightness at the half-light radius, $f_{\rm{e,gal}}$, is related to the flux of the galaxy, $I_{\rm{gal}}$, by $I_{\rm{gal}}\approx 7.2\pi\theta_{\rm{e},gal}^{2}f_{\rm{e,gal}}$.

Modeling the source plane emission with the surface brightness of a single galaxy can be a large simplification compared to real lenses. In particular, emission or absorption regions in the source galaxy can significantly change the appearance of the resulting arcs and arclets. However, in this work we aim to construct a simple test bed to compare probabilistic cataloging to mainstream analysis tools of strong lensing. Therefore, we leave a pixel-based representation of the source plane to future work.

2.3 Light emission from the foreground galaxy

We use the same surface brightness profile in Equation 13 to represent emission from the host galaxy, but with a different half-light radius and surface brightness, $\theta_{\rm{e,hst}}$, and $f_{\rm{e,hst}}$.

Because the scattering cross section between dark and baryonic matter is not necessarily zero, in principle there can be a spatial offset between the surface brightness and mass density of a host galaxy that has not been gravitationally relaxed (Shu et al., 2016). We do not model such unrelaxed systems, and assume that the smooth mass distribution and the light emission of the host galaxy are co-centric.

Furthermore, we do not subtract an estimate of foreground galaxy or isotropic emission from the observed photon-count map. This is because subtracting a best-fit emission component from the photon count map is equivalent to placing a delta function prior at the associated value, or equivalently, assuming that the parameter has vanishing covariance with all other parameters. Instead, we sample from the joint posterior probability distribution of all parameters including those of the host galaxy, and later marginalize over the nuisance parameters.

2.4 Background emission

Any CCD image includes a spatially uniform contribution from instrumental and astronomical backgrounds. We model this by adding an isotropic emission template with a floating normalization to the predicted images. This template also absorbs any true sky background in the vicinity of the strong lensed system. Note that the shot noise is built into the model by evaluating the Poisson probability of realizing the number of counts in the observed image, given the number of photons predicted by the model.

2.5 Instrumental PSF

Due to atmospheric distortions or imperfections in the optics, an optical photon-collecting device spreads an input delta function to a characteristic shape on the image plane, known as the PSF. In general, this shape can be modeled as a linear combination of Gaussians or (when the PSF has heavy tails) Lorentzians. A frequently used distribution is that of Student’s t-distribution (also known as the Moffat profile), which can account for both the Gaussian core and the heavy tails of the PSF shape.

In this paper, we simulate an HST image. In the absence of an intervening atmosphere and within the limit of distortion-free optics, we assume that the PSF is instead limited by diffraction, although the real HST PSF has diffraction spikes due to the support vanes of the secondary mirror. Given the circular aperture of the HST, we model the PSF in band $i$ as an Airy pattern whose first zero-crossing occurs at a radius of $\sigma_{\textrm{psf},i}$.

In general, uninformative priors can be placed on the parameters that characterize the PSF shape, e.g., $\sigma_{\textrm{psf},i}$. This would allow PSF parameters to be inferred along with the lens metamodel. However, in the absence of a bright point source in the image, the PSF cannot be constrained-well. In fact, for a given telescope, the shape of the PSF is usually well-constrained using independent calibration data. Therefore, PCAT assumes Gaussian priors on the PSF parameters, whose means and standard deviations are provided by the external calibration. This allows exploitation of previous calibration while also still allowing propagation of uncertainties due to potential covariances between the PSF and model parameters.

2.6 Model image

Given the above ingredients, we express the model image (prior to the PSF convolution) as the sum of

a spatially uniform template to model detector background and isotropic sky emission, $f_{\rm{bac}}$, 2. 2.

emission from the host galaxy, $f_{\rm{hst}}$, 3. 3.

the gravitationally lensed emission of a background light source due to a foreground host galaxy and a variable number of its subhalos, $\tilde{f}_{\rm{src}}$.

The total model emission, $f_{\textrm{m}}$, is then obtained by convolving this image with the PSF kernel, $\mathcal{P}$, such that

[TABLE]

2.7 Model indicator

Our lens metamodel is a union of lens models, each with a fixed number of subhalos. This implies that the number of subhalos, $N_{\rm{sub}}$, which is also known as the model indicator, is a discrete parameter of the metamodel. Through transdimensional proposals that respect detailed balance, it gets updated as the Markov chain state evolves, generating a posterior distribution over $N_{\rm{sub}}$. We assume that $N_{\rm{sub}}$ is Poisson-distributed with a mean of $\mu_{\rm{sub}}$ and place a uniform prior on the distribution of $\mu_{\rm{sub}}$. Therefore, in principle, $N_{\rm{sub}}$ can be any non-negative integer. However, we place an upper limit of 100 on it, for computational convenience, and check that this upper limit is never reached. Implementation of transdimensional proposals in PCAT is discussed in Section 3.

2.8 Summary of the lens metamodel

Table 1 lists the metamodel parameters. The group of rows at the top contains the hyperparameters, which parameterize the conditional distribution of subhalo properties. These subhalo parameters are shown in the group of rows at the bottom of the list. The remaining groups of rows are the PSF, background, and lens parameters, from top to bottom. The lens parameters are further divided into three subgroups, separately showing the foreground host lens, background light source, and external shear parameters.

The third and fourth columns indicate the minimum and maximum when the associated distribution is uniform, log-uniform or power-law, and show the mean and standard deviation when the distribution is Gaussian. The fifth column indicates the parameters of the metamodel from which our mock data sets are drawn. Generation of mock data sets will be discussed in Section 4.1.

The hierarchical nature of probabilistic cataloging can be seen in Figure 4, where we show the probabilistic graphical model associated with the lens model. The nodes in the graph represent random variables with a certain probability density, and directed lines indicate the conditional relation between random variables. Green nodes correspond to subhalo parameters, whose multiplicity, $N_{\rm{sub}}$, is another (discrete) parameter shown with blue. Red nodes denote the hyperparameters, which control the distribution of subhalo parameters such as multiplicity and distribution of Einstein radii. For instance, the number of subhalos, $N_{\rm{sub}}$, is assumed to be a Poisson realization of an underlying mean number of subhalos, $\mu_{\rm{sub}}$. Similarly, the deflection strengths, $\alpha_{\rm{s},\textit{n}}$, are distributed as a power-law with slope $\beta$.

The yellow nodes are the fixed-dimensional parameters, which affect the model image. Here, $\vec{A}$ collectively represents the set of parameters that control the background emission. In this work, we assume that the background is spatially uniform and data (hence the model) is in a single energy band. Therefore, $\vec{A}$ is reduced to a single parameter, $A$, which is the amplitude of spatially uniform emission in units of erg/cm2/s/sr. We collectively denote all the macro lens parameters with $\vec{\chi}$. We use $\vec{\eta}$ to denote the set of parameters that characterize the PSF of the instrument

3 Probabilistic Cataloging

In this section, we give an overview of probabilistic cataloging and refer the reader to (Brewer et al., 2013; Daylan et al., 2017; Portillo et al., 2017) for details. Probabilistic cataloging is a transdimensional sampling scheme, where an a priori unknown number of elements, each with certain features, are inferred from an observed data set, e.g., a series of images, through evaluation of the Poisson likelihood of a model given the data. We use the python2.7 implementation of the PCAT framework, which is publicly available at https://github.com/tdaylan/pcat. PCAT is a parallelized Poisson mixture sampler originally designed to sample from a gamma-ray emission metamodel consistent with the Fermi-LAT data. In this work, we extend PCAT to evaluate lens models.

Given some Poisson-distributed count data, PCAT proposes reversible jumps (Green, 1995) across members of the metamodel in the form of births, deaths, splits, and merges of transdimensional elements, which represent dark, light-deflecting subhalos in this context. Reversibility of the transdimensional proposals ensure that detailed balance is respected across the metamodel. The resulting chain contains fair samples from the posterior distribution of the metamodel given the data, which can be used to compare models or constrain their parameters.

When sampling from the posterior probability density of the metamodel, we assume that a photon-counting experiment, e.g., a CCD, measures impinging photon counts in spatial pixels and spectral bands. We denote the observed number of photons in energy band $i$ and pixel $j$ by $k_{ij}^{\rm d}$. Given an observed image, ideally our aim would ideally be to infer the underlying true model, $k_{ij}^{\rm true}$, such that $k_{ij}^{\rm d}$ is a Poisson realization of $k_{ij}^{\rm true}$. In practice, however, the true model can be highly complex. Therefore, we approximate it using a parametric metamodel, $k_{ij}^{\rm m}$, that predicts the mean photon count in energy band $i$ and pixel $j$. The observed photon count map can therefore be written as a Poisson realization of the photon count map predicted by the approximate model, yielding the log-likelihood

[TABLE]

where $\mathcal{D}$ represents the set of observed photon count maps and $\mathcal{M}$ denotes those predicted by the metamodel. The photon count map predicted by the metamodel in energy band $i$ and pixel $j$, $k_{ij}^{\rm m}$, is then the model flux convolved with the exposure of the detector toward a particular direction on the sky, ($\theta_{1},\theta_{2}$), at energy $E$, $\epsilon(E,\theta_{1},\theta_{2})$, and the transmission efficiency of the optical filter, $T(E)$,

[TABLE]

where $f_{\textrm{m}}(E,\theta_{1},\theta_{2})$ is the image predicted by the metamodel, i.e., number of photons per area, time, energy, and solid angle. In other words, we project the prediction of the lens metamodel onto the space where data are measured. In this space of pixels and energy bands, measurement uncertainties are simply Poissonian, which makes error propagation straightforward. Once $k_{ij}^{\rm m}$ is calculated over all pixels and energy bands, the Poisson log-likelihood of observing the image, given the data, is summed over all pixels and energy bands as in Equation (16), which yields the test statistics that can be used to compare models. Note that the term $\ln k_{ij}^{\rm d}$ is fixed for a given data set, and is neglected when comparing the likelihood of different models for the same data set.

Because inference is based directly on the observed image and involves a fairly large model space, probabilistic cataloging is a time-consuming method. However, parallel computing and fast evaluation of the log-likelihood significantly reduces the required execution times. For the results presented in this paper, 10 chains with $2\times 10^{6}$ samples each were collected in parallel over 8 hr (72 CPU-hours). By evaluating the likelihood in a compiled language, the execution speed can be increased by a factor of $\sim 5$ in the future to allow sampling simultaneously from the catalog space of multiple strong lens systems in the future.

4 Results

In order to validate the inference performance of PCAT, we run it on mock (simulated) data. We generate these mock count maps as Poisson realizations of the photon count maps predicted by forward-modeling a metamodel, which we refer to as the true metamodel. Because we know the parameters of the true metamodel, we can compare the posterior probability distributions sampled by PCAT with the underlying true parameters. This approach yields full control over systematic errors by allowing us, for instance, to fit a data set that has been drawn from a different metamodel.

4.1 Mock Data Set

When simulating an image, we assume that the data are taken using the WFC3/UVIS detector on the Hubble Space Telescope (HST). This choice is only intended to provide a concrete example, because probabilistic cataloging can be applied to any photometric data set such as ground-based, high-resolution images that use adaptive optics. We further assume that the F814W filter is used to collect photons. For this HST band, we multiply differential flux predicted by the metamodel, i.e., in units of erg cm*-2* s*-1* Å*-1* sr*-1*, by $6.1\times 10^{18}$ erg*-1* cm2 s Å(Ryan et al., 2016), as well as the pixel area in order to obtain the number of photons predicted per pixel. We assume that the photometric data are collected in a single exposure and that no drizzling is applied. Furthermore, we fix the observation time to 1000 seconds, which is roughly equivalent to half an orbit of HST. Choosing the mock galaxy brightnesses low enough allows one to ensure that the CCD does not sature over the selected exposure time. Finally, defining SNR as the ratio of the lensed surface brightness with the square root of the total surface brightness, our nominal mock data set has a maximum per-pixel SNR of $9$, which quickly decays away from the multiple images of the background source.

We sample the parameters of the true metamodel randomly from its prior, unless stated otherwise in Table 1. In particular, we choose the minimum of the true distribution of $\alpha_{{\rm s},n}$ to be 0.003 arcsec, and fix the number of true subhalos to 25. These two choices are made so that the resulting true subhalo mass fraction, averaged within 0.1 arcsec of the critical curve (see the discussion on the subhalo mass fraction in Section 4.2) is $\sim 2\%$ (Gao et al., 2011). Furthermore, we assign deterministic values to some parameters as indicated under the true column of the table. The letter combination RC under this column indicates that the associated parameters (galaxy coordinates) are drawn from Gaussians at the center of the image with a standard deviation equal to the pixel size (0.04 arcsec). Likewise, R denotes that the parameter is drawn randomly from the prior. Once the true metamodel parameters are determined, the predicted image is calculated and a Poisson realization of the map is drawn to obtain the mock photon count map.

Note that the true subhalos in this mock data set are mostly low-significance subhalos. Out of the 25 true subhalos, only 3 improve the goodness of fit, $\Delta\log P(D|M)$, by more than 35 when added to the best-fitting model, meaning that the rest are below formal detection ($5\sigma$ for five degrees of freedom).

4.2 Nominal Results

We begin by showing, in Figure 5 one realization of the mock photon count map (repeated across panels) and six fair samples from the ensemble of subhalo catalogs, represented with blue pluses. The true catalog of subhalos is indicated by the green markers.

The multiple images in the simulated data set are formed by a background galaxy near the fold caustic of an elliptical foreground galaxy. The green right triangle and the blue left triangle show the position of the true and model background source, respectively. Similarly, the green diamond and the blue filled square indicate the positions of the true and model host galaxy on the lens plane, respectively. The positions of the true and model galaxies are overlapping in all samples, although the positions of the blue markers have unnoticeable variations across samples due to the rather small uncertainties in the posterior distributions of the source and host galaxy positions.

Macro lens model

In this work, we focus on the inference of subhalo properties in the lens metamodel. Nevertheless, subhalos have only a perturbative effect on the deflection field, which is otherwise dominated by the host galaxy. Therefore, any bias in the host galaxy parameters, as well as other macro lens parameters – such as those of the background galaxy and the external shear – can potentially leak into constraints on the subhalo properties. We therefore first discuss the posterior distribution of the macro lens model parameters.

In general, we find that the best-constrained macro lens parameters are the positions of the host and source galaxies as well as the Einstein radius of the host galaxy. They generally have posterior 68$\%$ credible intervals that are $6\times 10^{-4}$ arcsec, $5\times 10^{-3}$ arcsec, and $0.01$ arcsec wide, respectively. However, we also find that the Einstein radius of the host halo correlates with its ellipticity, the amplitude and orientation angle of which are constrained inside 68% credible intervals that are 0.01 and 2∘ wide, respectively. This anti-correlation can be attributed to the fact that the critical curve is constrained more strongly along one axis of symmetry than the other.

As an example, Figure 6 shows the distribution of posterior samples on the horizontal position of the host halo. The green vertical line indicates the true value of the parameter, whereas the red vertical line shows the value of the parameter in the maximum likelihood sample. The dashed and dotted-dashed vertical lines show the 95% and 68% uncertainties around the median.

The host and source galaxy fluxes, on the other hand, are constrained at the $\sim 2\%$ and $\sim 10\%$ level, respectively. Specifically, the fact that we impose the background galaxy to have a Sérsic profile of known index reduces the uncertainty on it. If the source plane emission had been allowed to vary from pixel to pixel, subject to some regularization, the uncertainties on the background emission would be even larger. Nevertheless, because we are mainly interested in modeling the foreground mass distribution, the background emission is eventually marginalized out, which reflects the uncertainty of the marginal posterior probability distribution of the subhalo properties. The choice of Sérsic profile also introduces a strong covariance between the luminosity and half-light radius of the galaxies, which can have correlation coefficients as large as 0.9.

The external shear is the least strongly constrained ingredient of the macro lens model, with typical uncertainties of $\sim 10\%$ on the magnitude of the shear field and $\sim 20^{\circ}$ error on its direction. The large uncertainty in the alignment is likely due to the degeneracy of the quadrupole shear term with the overall deflection field due to subhalos. Similarly, the posterior probability distribution over the lens metamodel is marginalized over the shear field and its uncertainties are propagated to the uncertainties on the subhalo properties.

Finally, we obtain $\sim 0.1\%$ uncertainties on the isotropic background emission and recover our prior uncertainty on the PSF width in the posterior distribution. The latter allows one to marginalize over uncertainties that are due to our lack of perfect knowledge of the PSF, and is not intended to gain information about the PSF using the image.

Subhalo mass function

We present the posterior truncated mass distribution of the subhalos in Figure 7. The figure shows the histogram of the truncated masses of the true subhalos in green. Overplotted with blue error bars, we show the median histogram of the truncated masses of the subhalos in the posterior samples. These blue points with error bars are obtained by first calculating the histogram of the truncated masses for each posterior sample from the metamodel. The 16th, 50th, and 84th percentiles in each bin are then plotted as the lower cap, central point, and the upper cap of the error bar, respectively.

It may initially seem surprising that even the highest-mass bins that contain a true subhalo may have a median posterior at zero, or that the posterior uncertainties go above zero when there is no true subhalo. This is a natural result of accounting for transdimensional (across-model) covariances in the problem. When births and deaths of subhalos are allowed, a true subhalo can sometimes be fitted with multiple and lower-mass subhalos. Similarly, multiple true subhalos can be fitted with a single more massive model subhalo. It is important to emphasize that such transdimensional covariances exist even for true subhalos with high statistical significance. In other words, even though the log-likelihood difference between two models with and without a given subhalo may be well above 35, it may be as low as a few when comparing two models where the subhalo exists as a single, high-mass clump and as two, smaller-mass clumps. Furthermore, it is also true that a high-mass subhalo can be far from the multiple images, causing the posterior subhalo mass function to be uninformed of its existence.

We note that the roll-off at small masses is due to the choice of the minimum deflection strength allowed in the model, $\alpha_{\rm{s,min}}$. If this value is decreased, the posterior probability density extends to lower masses. Therefore, the constraints derived on the subhalo mass function must be quoted for a given $\alpha_{\rm{s,min}}$. Given our choice of the mock data set and fiducial redshifts, our posterior subhalo mass function becomes prior-driven below $\sim 10^{8}M_{\odot}$. This is illustrated in Figure 8, where we plot the truncated mass and significance of subhalos in a different mock data set with a higher number of subhalos.

As for the deflection profile parameters $\theta_{{\rm s},n}$ and $\theta_{{\rm c},n}$, we generally find that the projected cutoff radius is better constrained than the projected scale radius. In fact, WFC3’s pixel size of 0.04 arcsec is comparable to the projected scale radius of an NFW subhalo with $r_{{\rm s},n}=200$ pc at the fiducial redshift of $z=0.2$. This implies that current optical photometry is insensitive to the scale radii of most subhalos. However, more massive subhalos close to the multiple images can still be distinguished by their projected scale radii, which motivates our choice to make $\theta_{{\rm s},n}$ subject to inference.

Note that the mass of a subhalo does not have a one-to-one correlation with its statistical significance. In a photometric inference problem, how the flux of a given point source compares to the background emission inside the full width at half maximum (FWHM) of the PSF, largely determines its statistical significance (Portillo et al., 2017). Therefore, we expect that isolated, bright light sources will be, on average, to be more statistically significant than dimmer ones. However, how strongly a subhalo can be constrained in the lensing problem depends on its mass, as well as how its deflection field is oriented with respect to the lensed emission. In other words, because the deflection profile of a subhalo drops beyond its projected scale radius, a subhalo needs to be close to an already lensed emission on the image plane as well as have a sufficiently large mass, in order to have non-negligible effect on the log-likelihood.

The leading-order contribution of a given subhalo with index $n$ to the model image can be obtained by Taylor-expanding Equation (1). Recasting this equation in terms of the photon count map lensed by the host halo, external shear, and all subhalos except the $n$th subhalo, $\bar{k}_{m}$, we obtain the lensed count map

[TABLE]

The second term is negative over some regions of the lensed emission, and positive in other regions. Instead, \Big{\langle}\lvert\vec{\alpha}_{n}\cdot\vec{\nabla}k_{m}\Big{\rvert}\rangle better encapsulates the level of perturbation that a given subhalo introduces to the image already lensed by all other mass components. Therefore, it is a relatively more accurate estimator of the statistical significance of a given subhalo compared to $M_{\rm{c},n}$. We refer to this quantity as the subhalo relevance. Subhalos with low relevance, although possibly real, cannot be constrained using the observed photometric data. An example can be seen in Figure 5, where the true subhalo at (1.8, 1.1) arcsec is relatively less constrained, despite having a mass above $10^{8}M_{\odot}$.

Figure 9 shows the correlation between the relevance of a subhalo and the log-likelihood increase when including the subhalo into the lens model while leaving all other parameters fixed. The correlation is not one-to-one, however, due to the fact that the host galaxy and the external shear can dominate the deflection field in some regions more than others, which causes subhalos with equal relevance to have different statistical significances. Therefore, in general, significant subhalos have high relevance, but highly relevant subhalos can have very low significance.

Subhalo mass fraction

The shape of the subhalo mass function is poorly constrained, given the small number of strongly lensed systems analyzed so far. Nevertheless, the fraction of mass tied in subhalos in the vicinity of the critical curve of the host halo can be accurately measured for strong lenses with high SNR. Therefore, the subhalo mass fraction provides an observational summary statistic of the subhalo mass function.

In order to determine the posterior fraction of mass locked in subhalos, we again first calculate it for each sample and then take posterior moments. We define the subhalo mass fraction for a given sample from the metamodel as the ratio between the mean surface mass density due to all subhalos and that due to the host halo inside an annuli centered at the host position with radius $\theta_{\rm{E,hst}}$ and thickness 0.2 arcsec,

[TABLE]

where “ann” represents the mentioned annulus around the host halo. The resulting posterior distribution of the subhalo mass fraction at the Einstein radius is shown in Figure 10 along with its evolution as a function of MCMC time. Note that we do not directly impose a prior probability density on the subhalo mass fraction, $f_{\rm{E,sub}}$. The figure shows that the posterior distribution of $f_{\rm{E,sub}}$ agrees with the true value, and that it is sensitive to the existence of subhalos below detection. The green line shows the true subhalo mass fraction.

Finally, we show the posterior deflection maps in Figure 11, which are obtained as follows. For each fair sample from the lens metamodel, we first calculate the deflection vector, i.e., deflection direction and magnitude, on all pixels. Then, for each pixel, we find the 50th percentile of the calculated deflection values, and refer to the resulting map as the posterior median deflection map. Figure 11 then illustrates the deflection field separately for the host halo (left), external shear (center), and subhalo population (right). The radius of the green and blue dashed lines, which are also mostly overlapping, denote the circles whose radii coincide with the Einstein radii of the true and model host halos, respectively. Because the host galaxy is elliptical, however, the circles do not correspond to the critical curve – they are only intended to guide the eye. The area of a circle is proportional to the host mass contained inside the Einstein radius. The true Einstein radius is 1.5 arcsec, yielding a mass inside the Einstein radius of 5$\times$1011 $M_{\odot}$.

The posterior total deflection field is found to be constrained at the percent level across most of the image. However, this translates to $\sim$ 100% uncertainty in the deflection field due to subhalos. The posterior deflection field due to subhalos is only constrained well near the multiple images, where the uncertainties can be as low as 20%.

4.3 Systematic bias in a one-subhalo lens model

Fixed-dimensional lens modeling requires iteratively testing lens models with a fixed number of subhalos (or a fixed amount of any other form of additional model complexity) in addition to the smooth lens model. In the simplest approach, a single subhalo can be added to the smooth lens model to fit an observed photon count map. In this section, we compare the previous results from transdimensional sampling to such a case, where the lens model contains a single subhalo. To make a robust comparison, we use the same simulated photon count map as shown in Figure 5. Then, instead of allowing a variable number of subhalos in the lens model, we turn off transdimensional proposals in PCAT and fix the number of subhalos to 1. Note that the latter scheme is still Bayesian and probes all within-model covariances, but neglects across-model covariances.

Figure 12 illustrates the posterior median convergence field obtained in the transdimensional (left) and one-subhalo (right) inferences. These posterior convergence maps are obtained similarly to the posterior deflection maps. In other words, the convergence maps are calculated for each sample from the posterior by numerically differentiating the deflection field, and then the median convergence in each pixel is presented as the posterior median convergence map.

In the left panel of Figure 12, transdimensionality allows subhalos to be born and die across the image, probing the goodness of fit of all combinations of subhalo multiplicities, positions and deflection profiles allowed by the prior. Such configurations mostly include mild likelihood improvements below what would be called a detection threshold in conventional cataloging. The resulting posterior median convergence map reveals low-significance diffuse features in the underlying subhalo convergence field, as well as constraining statistically significant true subhalos. In contrast, the posterior median convergence map obtained from the one-subhalo inference on the right contains a single overdensity, which is also inferred by the transdimensional inference. However, the transdimensional approach also (partially) deblends subhalos in the crowded region, to the extent allowed by the information available in the data and degeneracies in the likelihood, despite the fact that most sample subhalos presented in Figure 5 do not indicate associations with the true subhalos.

The failure of the one-subhalo inference to reveal the remaining features in the underlying true convergence field can be traced back to the fact that the absence of additional model subhalos causes the macro lens model and the single subhalo to be biased so as to minimize the residual deflection field. Without a way to probe transdimensional covariances, the sampler explores a maximum in the conditional likelihood, given a one-subhalo lens model. Note that the inference of lower-significance features in the convergence map does not imply that the additional subhalos are formally detected. In fact, the overall increase in the goodness of fit between the two approaches is only $\Delta\log P(D|M)\sim 10$. In order to claim a detection of these model subhalos at $5\sigma$, they would be required to individually improve the goodness of fit above $\Delta\log P(D|M)\sim 35$ given that they have five degrees of freedom.

The difference between the two inference schemes becomes more evident when the posterior predictions are compared to the underlying true convergence field due to subhalos. Figure 13 illustrates the residual between the true subhalo convergence map and the posterior median convergence maps obtained in the transdimensional (left) and one-subhalo (right) inferences, respectively. Note that the disagreement is significantly smaller in the case of a transdimensional approach, which correctly predicts the mean subhalo convergence field over most of the lens plane and has smaller disagreement close to the subhalos. However, note that even the transdimensional approach cannot fully constrain the true subhalo convergence field in the vicinity of most subhalos. This is expected because the simulated photon count map contains very little information on their properties.

We further show the posterior distribution of the positions, and deflection strength of the model subhalo in the one-subhalo inference in Figure 14. It shows that the joint distribution is highly skewed with heavy tails. Furthermore, the posterior median of the corresponding inferred mass of the subhalo, $4.5\times 10^{8}M_{\odot}$, is biased low compared to the mass of the closest true subhalo at (1.2, -0.5) arcsec, $8\times 10^{8}M_{\odot}$.

Another problem with fitting the lens image using a fixed number of subhalos is that the posterior distribution on the macro lens model can be significantly biased, depending on the orientation of the non-modeled subhalos. This is because the macro lens model compensates for the absence of subhalos below the detection threshold, resulting in a seemingly good fit that mismodels the underlying true lens metamodel. An example is shown in Figure 15, where the posterior distribution of the host halo ellipticity is plotted for the cases of nominal (transdimensional) and one-subhalo inferences, using a similar but independently drawn data set. One-subhalo inference is observed to generate multimodality in the posterior distribution of the ellipticity of the host.

4.4 Extending the metamodel to lower subhalo masses

The lower deflection strength limit of probabilistic cataloging can be made smaller at the expense of making the posterior prior-dominated (thus resulting in no information gain) and significantly increasing the computation time as required for MCMC convergence. Therefore, given the observed data, there is a reasonable lower limit to the deflection strength, e.g., 0.01 arcsec, that can balance information gain and computational complexity.

Nevertheless, for self-consistency, we illustrate an inference where there is no mismodeling at the low end of the subhalo mass function, unlike our nominal results in Section 4.2, where the minimum subhalo deflection strength allowed in the fitting metamodel was a factor of three higher than the true minimum subhalo deflection strength. Using the same data set as shown in Section 4.2, Figure 16 compares the posterior distribution of the number of subhalos in our nominal results (left) and in a run where $\alpha_{\rm{s,min}}=0.003$ arcsec (right). Given a lower $\alpha_{\rm{s,min}}$, the run shown in the right panel explores a higher number of subhalos whose masses are lower, on average. Note that the two posteriors are consistent with each other. Because we have a uniform prior on the expected number of subhalos, $\mu_{\rm{sub}}$, the fact that there is a preferred scale in the number of subhalos – and it is consistent with the true number of subhalos – implies that the posterior is informed by the underlying subhalos below the detection threshold.

5 No-signal (null) test

Any inference framework is expected to produce null results when subject to an input data set that does not contain the signal of interest. We therefore tested our inference framework on a mock data set that was not affected by any subhalo. For this purpose, we generated another mock image as a Poisson draw of an image forward-modeled by the macro lens model, background emission, and PSF (same as that used to generate the main data set presented in the paper), but without any subhalo. We then run PCAT on this image. Figure 17 shows the histogram of the number of subhalos in the posterior, which is consistent with the true value of 0. Furthermore, none of the subhalos in the posterior ensemble of catalogs is more significant than 4$\sigma$. The fact that there is a non-zero number of subhalos in the posterior is largely a consequence of the degeneracy between the mass of the main deflector and the total mass of subhalos inside the critical curve. The total mass in the subhalos, as shown in Figure 18, is also mostly zero, although the posterior does contain some mismodeled subhalo mass due to degeneracies in the lensing problem. The 84th percentile (vertical blue dashed-dotted line) extends to about 3$\times 10^{8}M_{\odot}$, which coincides with the maximum likelihood sample (vertical red line), indicating that these few samples overfit the image.

6 Sensitivity to the choice of $\alpha_{\rm{s},min}$

An important nuisance parameter of our model is the minimum deflection strength allowed by the prior. This is critical because the subhalo mass function is expected to be steep, causing the posterior total number of subhalos to sensitively depend on the minimum deflection strength allowed by the metamodel. The lower $\alpha_{\rm{s},min}$ is allowed to go, the more prior-dominated the posterior becomes. Conversely, taking $\alpha_{\rm{s},min}$ too high causes mismodeling because those subhalos with marginal significances that could fit subtle features in the data are thus removed from the metamodel.

In Figure 19 we illustrate the effect of varying this nuisance parameter away from the nominal value of 0.01 arcsec. In the left panel, $\alpha_{\rm{s},min}$ is reduced to 0.005 arcsec, whereas it is increased to 0.02 arcsec in the right panel. Although the posterior in the left panel agrees better with the ground truth at the low-mass end, this happens at the expense of oversplitting the most massive subhalos. Furthermore, the lower end of the posterior subhalo mass function becomes prior-dominated, and could miss a potential cutoff or flattening required by the data. The 0.02 arcsec run on the right performs comparably at the high masses, but fails to explore the subhalo mass function at the low-mass end.

7 Discussion

Our results using mock data show that probabilistic cataloging is an effective tool to probe lens models with a variable number of deflectors. Our approach significantly differs from conventional cataloging methods because we do not discard information or bias the posterior by imposing a detection threshold on model subhalos, i.e., requiring that model subhalos improve the goodness of fit above some threshold. By subsequent marginalization over the catalog space of subhalos, we obtain robust estimates of the underlying mass model and population characteristics of subhalos, without necessarily individually detecting them. This contrasts with the conventional cataloging, where a few formally detected subhalos and nondetections are used to constrain the subhalo mass function. However, because these detections are not fair draws from the underlying subhalo mass function, this approach cannot account for covariances in the subhalo catalog space, potentially yielding biased results with erroneously small uncertainties.

Another key ingredient of probabilistic cataloging is the exploration of across-model covariances, which is not accounted for in conventional cataloging. For example, given an observed lensed image with a single underlying subhalo, the resulting photon count features can be fitted by a single high-mass subhalo, or multiple lower-mass subhalos. Moreover, because of the inherent degeneracies of the lensing problem, an uncountable set of other configurations can be consistent with the observed data, depending on the positions of subhalos with respect to the lensed emission. Point estimates in the hypothesis space can therefore lead to biased results. Probabilistic cataloging offers a principled way of characterizing these covariances and propagating uncertainties.

In principle, probabilistic cataloging of subhalos could be used to sample arbitrarily low-mass subhalos. However, because the number of subhalos to be modeled steeply increases toward lower subhalo masses, convergence of chains with many low-mass subhalos becomes an issue. Furthermore, because current data sets can only allow the detection of only the few most massive subhalos, if any, the constraints on the subhalo mass function at the low-mass end are largely driven by the prior because the likelihood is essentially flat there.

Currently, the most extensive and homogeneous collection of optical images of galaxy-galaxy type strongly lensed systems is the SLACS collection (Bolton et al., 2005, 2008). These targets have been spectroscopically selected from the Sloan Digital Sky Survey (SDSS) data set conditional on encountering multiple strong emission lines at a redshift higher than that of the target galaxy. They were later followed up by the ACS detector on the HST. However, the limited depth of the HST only allows the most massive subhalos, if any, to be constrained (Vegetti et al., 2010, 2012). Next-generation, wide-field, high angular resolution telescopes such as WFIRST and EUCLID are expected to revolutionize studies of dark matter substructure. The number of high-SNR strong lens images is expected to go up to thousands (Treu, 2010; Oguri & Marshall, 2010). Furthermore, greater exposure depth should allow constraints on multiple subhalos per strong lens, making it possible to constrain the subhalo mass function to lower masses. These infrared (IR) telescopes will also carry spectrometers that, in synergy with the photometric redshifts provided by LSST, will allow simultaneous constrains on the redshifts of the model sources and deflectors. Follow-up of high-SNR systems with the JWST and interferometry at submillimeter wavelengths, e.g., ALMA, may allow constraints on the subhalo properties of systems at high redshift.

The current framework we present in this paper can be improved in several ways. First, multi-band images can be incorporated into the framework in a straightforward manner, by introducing color parameters for the source and host emissions. This would make available information contained in the different colors of sky background, source, and foreground galaxies. Furthermore, the source plane emission can be parameterized on an adaptively refined grid, instead of imposing that the source plane emission is from a galaxy with a perfect Sérsic profile. This would yield a more realistic model of the source plane emission, and a more principled method to marginalize over uncertainties in the emission from the source plane. This is especially relevant for late-type background galaxies, where emission has strong spatial features, such as outflows, bright spots, or dust obscured regions on the disc. We leave this study, along with the application of our method to real data, to future work.

Probabilistic cataloging can also be extended to use temporal information by sampling from the posterior distribution of a lens metamodel given exposures of the strongly lensed system taken at different times. Because lensing is also sensitive to the angular diameter distances to the lens and source planes, this would allow an independent measurement of the Hubble constant, $H_{0}$, as in Suyu et al. (2017), but with constraints that have been marginalized over dark subhalos that can potentially bias the macro lens parameters. Also, combining multiple exposures may improve sensitivity to substructure via time-delay lensing (Keeton & Moustakas, 2009; Cyr-Racine et al., 2016), although microlensing could potentially bias this inference (Tie & Kochanek, 2017).

Furthermore, probabilistic cataloging provides a principled way to gather information from independent strong lens systems, when constraining the probability distributions of subhalo properties. Given images of different strong lenses, the catalog space consistent with these images can be jointly sampled with a common hyperparameter set characterizing the population characteristics of the subhalos. This is a significant improvement with respect to traditional cataloging, where incorporating nondetections into the inference may become ambiguous because fixed-dimensional approaches do not consistently account for the fact that non-detections can contain real subhalos or that detections may contain none, two or, more.

In conclusion, we suggest a transdimensional, hierarchical, and Bayesian framework to infer the subhalo mass model and population characteristics in strong lens systems. Using gravitational lensing as a probe, the framework models the lens plane mass distribution and background light emission to draw fair samples from the posterior probability distribution over the lens metamodel, given some photometric data set. The resulting posterior can be marginalized to constrain the subhalo mass fraction and the mass distribution on the lens plane.

We thank Christopher W. Stubbs, Daniel Eisenstein, Ana Bonaca, Priyamvada Natarajan, Manoj Kaplinghat, and Stephen K.N. Portillo for useful discussions during the project. We also appreciate useful feedback from Simon Birrer, Simona Vegetti, and Brendon J. Brewer during the review process. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. F.-Y. C.-R. acknowledges the support of the National Aeronautical and Space Administration ATP grant NNX16AI12G at Harvard University.

Appendix A Deflection profile of subhalos

When projected onto the lens plane, the deflection field due to an untruncated, spherically symmetric NFW subhalo becomes azimuthally symmetric and can be described by a radial profile (Golse & Kneib, 2002)

[TABLE]

where $\bar{\Sigma}_{n}(\theta)$ is the mean surface mass density of the $n$th subhalo inside a subhalo-centric circle of radius $\theta$, and

[TABLE]

is the critical surface mass density on the lens plane, which has a value of 3.1$\times 10^{9}M_{\odot}$ kpc*-2* for our simulated data. For notational brevity, we use the rescaled angular distance, $\theta^{\prime}\equiv\theta/\theta_{{\rm s},n}$. The function $F(\theta^{\prime})$ is given by

[TABLE]

When the truncation is taken into account, the deflection profile becomes (Baltz et al., 2009)

[TABLE]

where

[TABLE]

In the limit of large $\tau_{n}=\theta_{{\rm c},n}/\theta_{{\rm s},n}$, this expression reduces to Equation (A1). Hence, the three parameters we use to describe the resulting deflection profile of the $n$th subhalo are the normalization $\alpha_{\rm{s},n}$, projected scale, and cutoff radii $\theta_{{\rm s},n}$, $\theta_{{\rm c},n}$.

Appendix B Markov chain convergence

When sampling from a target probability distribution by constructing Markov chains, one needs to ensure that all walkers reach the same stationary distribution. Toward this purpose, PCAT monitors the Gelman-Rubin (GR) test statistic, which estimates the factor by which the overall statistical uncertainty of a random variable could be reduced by running the chains longer. Because the chain is transdimensional, however, the statistic is calculated over the data space, i.e., by evaluating the statistic for the number of photons predicted by the metamodel in each pixel.

Figure 20 illustrates the GR test statistic evaluated over the image pixels for a typical PCAT output. At fixed number of samples, the across-chain variance generally increases with the typical number of subhalos. Therefore, given a choice of $\alpha_{\rm{s},min}$, and hence a typical number of subhalos, we ensure that the chain is run long enough to reduce the GR test statistic down to $\lesssim 1.1$. For a typical run with $20-40$ model subhalos, this requires $20\times 10^{6}$ samples to be taken before burn-in and thinning. This also implies that the reported constraints have a $\lesssim 10\%$ uncertainty that adds in quadrature to the statistical uncertainty that we derive by assuming that all chains have sampled from the same target distribution, i.e., the posterior distribution over the desired metamodel.