Probing dark matter substructure in the gravitational lens HE0435-1223 with the WFC3 grism

TL;DR

This study uses WFC3 grism observations of the quadruply lensed quasar HE0435-1223 to place constraints on dark matter substructure, demonstrating the potential of narrow-line flux ratios to test Cold Dark Matter models.

Contribution

First application of WFC3 grism narrow-line imaging to probe dark matter substructure in a lensed quasar, providing a new method to detect or rule out low-mass halos.

Findings

No significant dark matter subhalos with M_{600} > 10^{8} M_ ext{sun} detected near the images.

Results are consistent with ΛCDM predictions for a single lens system.

Method demonstrates high sensitivity, promising for future dark matter studies with ~20 lenses.

Abstract

Strong gravitational lensing provides a powerful test of Cold Dark Matter (CDM) as it enables the detection and mass measurement of low mass haloes even if they do not contain baryons. Compact lensed sources such as Active Galactic Nuclei (AGN) are particularly sensitive to perturbing subhalos, but their use as a test of CDM has been limited by the small number of systems which have significant radio emission which is extended enough avoid significant lensing by stars in the plane of the lens galaxy, and red enough to be minimally affected by differential dust extinction. Narrow-line emission is a promising alternative as it is also extended and, unlike radio, detectable in virtually all optically selected AGN lenses. We present first results from a WFC3 grism narrow-line survey of lensed quasars, for the quadruply lensed AGN HE0435-1223. Using a forward modelling pipeline which enables us to robustly account for spatial blending, we measure the [OIII] 5007 \AA~ flux ratios of the four images. We find that the [OIII] fluxes and positions are well fit by a simple smooth mass model for the main lens. Our data rule out a $M_{600}>10^{8} (10^{7.2}) M_\odot$ NFW perturber projected within $\sim$1\farcs0 (0\farcs1) arcseconds of each of the lensed images, where $M_{600}$ is the perturber mass within its central 600 pc. The non-detection is broadly consistent with the expectations of $\Lambda$CDM for a single system. The sensitivity achieved demonstrates that powerful limits on the nature of dark matter can be obtained with the analysis of $\sim20$ narrow-line lenses.

Figures

Tables

Table 1. Table 1: Gravitational lens model parameters for the main deflector and external shear inferred from image positions and [OIII] fluxes which are given in Tables A1 and A2 respectively.

Parameter	Value	Description
${\theta }_E$	1$\stackrel{\prime \prime }{.}$200$\pm 0.003$	Spherical Einstein radius
q	$0.91\pm 0.03$	b/a
PA	$-8\pm 5$	Degrees E of N
${\gamma }_{ext}$	$0.063\pm 0.007$	External shear amplitude
${\theta }_{\gamma }$	$-18\pm 2$	Direction of external shear (Deg. E of N)

Table 2. Table 2: F140W image positions in units of arcseconds with North up, East left coordinates (rotated relative to Figure 1).

Image	dRa	dDec	Uncertainty
A	2.476	0.608	0.008
B	0.997	1.157	0.008
C	0	0	0.008
D	1.530	-1.005	0.008
G	1.314	0.067	0.06

Table 3. Table 3: Subset of flux ratio measurements of HE0435-1223 from the literature, selected to represent variation across wavelength and with time, plotted in Figure 3. a Formal measurement uncertainties for Wisotzki et al. ( 2003 ) less than 0.1 0.1 0.1 % and dominated by unknown systematics.

Filter	Date	$f_A/f_C$	$f_B/f_C$	$f_D/f_C$	Reference
$u^{\prime }$	9/2007	2.36$\pm 0.07$	$0.93\pm 0.04$	$1.03\pm 0.05$	Blackburne et al. (2011)
CIV	9/2002	1.40a	1.05	0.770	Wisotzki et al. (2003)
$g^{\prime }$	9/2007	1.95$\pm 0.1$	$0.96\pm 0.08$	$0.9\pm 0.07$	Blackburne et al. (2011)
$V$	8-9/2009	$1.37\pm 0.07$	$0.97\pm 0.04$	$0.79\pm 0.04$	Ricci et al. (2011)
F555W	8/2003	$1.84\pm 0.1$	$1.08\pm 0.09$	$0.95\pm 0.06$	Kochanek et al. (2006)
$R$	10/2007	1.89$\pm 0.05$	$0.99\pm 0.04$	$0.88\pm 0.05$	Blackburne et al. (2011)
$R$	10-12/2009	$1.42\pm 0.08$	$0.93\pm 0.04$	$0.75\pm 0.04$	Courbin et al. (2011)
$i^{\prime }$	9/2007	$1.67\pm 0.09$	$0.95\pm 0.07$	$0.77\pm 0.06$	Blackburne et al. (2011)
F814W	8/2003	$1.69\pm 0.04$	$1.12\pm 0.04$	$0.82\pm 0.03$	Kochanek et al. (2006)
$z^{\prime }$	9/2007	$1.64\pm 0.05$	$0.96\pm 0.04$	$0.8\pm 0.04$	Blackburne et al. (2011)
$J$	9/2007	$1.57\pm 0.09$	$0.99\pm 0.08$	$0.69\pm 0.06$	Blackburne et al. (2011)
H$\beta$	8/2015	1.36 $\pm 0.04$	1.00 $\pm 0.04$	0.78$\pm 0.03$	this work
[OIII]	8/2015	0.97$\pm 0.07$	0.98$\pm 0.07$	0.66$\pm 0.05$	this work
F160W	8/2003	$1.57\pm 0.05$	$1.0\pm 0.03$	$0.79\pm 0.03$	Kochanek et al. (2006)
F160W	10/2012	$1.30\pm 0.05$	$0.92\pm 0.05$	$0.66\pm 0.05$	Wong et al. (2017)
$K_s$	9/2007	$1.39\pm 0.03$	$0.98\pm 0.03$	$0.77\pm 0.02$	Blackburne et al. (2011)
$K$	8/2008	$1.84\pm 0.1$	$1.37\pm 0.08$	$0.74\pm 0.06$	Fadely & Keeton (2011)
5 GHz	11/2012	$1\pm 0.1$	$0.8\pm 0.07$	$0.47\pm 0.07$	Jackson et al. (2015)

Full text

††pagerange: Probing dark matter substructure in the gravitational lens HE0435-1223 with the WFC3 grism–3††pubyear: 2016

1 Introduction

A key prediction of Cold Dark Matter (CDM) is that the halo mass function should follow an unbroken power-law with $dN/dM\propto M^{-\alpha}$, where $\alpha=1.9\pm 0.1$ from cluster down to planet mass halos (Diemand et al., 2008; Springel et al., 2008). CDM models of the distribution and growth of structure match observations with remarkable success over an enormous range of distance and size scales (e.g. Planck Collaboration et al., 2014). On smaller scales, tests of CDM become more difficult owing to the uncertain physics of star formation in sub-Milky Way mass halos ($M_{vir}<10^{12}M_{\odot}$). A famous example of this is the ‘Missing Satellite Problem’, so called due to the fact that CDM simulations predict that thousands of subhalos should be gravitationally bound to the Milky Way, while only $\sim$ tens of luminous satellite galaxies have been observed (Moore et al., 1999; Klypin et al., 1999; Strigari et al., 2007; Weinberg et al., 2008; Drlica-Wagner et al., 2015). Extrapolations based on the depth and area of the Sloan Digital Sky Survey (SDSS) and the Dark Energy Survey (DES) forecast that the Milky Way may host as many as $\sim$100 ultra-faint ($L<10^{3}L_{\odot}$) satellite galaxies (see, e.g. Hargis et al., 2014; Drlica-Wagner et al., 2015). However, measurements of the halo mass function at these low luminosities becomes extremely difficult. In particular, the stars of low luminosity galaxies occupy only the inner $\sim 100$ pc of their dark matter halo. This is a small fraction of the total halo virial radius, which, in absence of the effects of tidal stripping extends out to several kpc even for a very low mass $10^{7}M_{\odot}$ halo. One consequence of this is that satellite galaxies inhabit halos with kinematically consistent masses over 5 orders of magnitude in luminosity (Strigari et al., 2008).

Simulations with varying implementations of semi-analytic and/or numerical baryonic feedback have demonstrated that it is possible to suppress star formation in dark matter subhalos sufficiently to matched the observed Milky Way satellite luminosity function with an underlying CDM halo mass function at low redshift (e.g., Thoul & Weinberg, 1996; Gnedin, 2000; Kaufmann et al., 2008; Macciò et al., 2010; Springel, 2010; Guo et al., 2011; Zolotov et al., 2012; Brooks et al., 2013; Starkenburg et al., 2013; Weinberg et al., 2008; Menci et al., 2014; Lu et al., 2014; Wetzel et al., 2016). However, models which can match the luminosity function of satellites at redshift zero around the Milky Way, are not always successful at reproducing the luminosity function around higher mass hosts, or at higher redshifts (e.g. Nierenberg et al., 2013, 2016).

Strong gravitational lensing provides a powerful test of CDM as it enables a measurement of the subhalo mass function without requiring stars or gas to detect subhalos. In a strong gravitational lens, a background source is multiply imaged with the image positions and magnifications depending on the first and second derivatives of the gravitational potential, respectively. The image magnifications are particularly sensitive to low-mass perturbations, with a lower sensitivity limit determined by the source size. If the source is of the scale of $\mu$as (such as a quasar accretion disk), then the image magnifications will be significantly affected by stars in the plane of the lens galaxy (a.k.a microlensing). In contrast, milliarcsecond scale sources are not significantly lensed by stars, but are sensitive to the presence of perturbing subhalos which have characteristic Einstein radii of $\sim$mas, corresponding to typical masses $10^{4}<M/M_{\odot}<10^{9}$ given typical lens configurations.

Traditionally, radio loud quasar sources have been used to detect subhalos, as they are extended enough (tens to hundreds of parsecs, e.g. Jackson et al. (2015)) to avoid microlensing, and also are not affected by differential dust extinction in the plane of the lens galaxy. Dalal & Kochanek (2002) used the lensed magnification of 6 radio loud quasar lenses and PG1115+08 (Weymann et al., 1980) to estimate the average fraction of mass in substructure relative to the mass in the smooth halo component ($f_{sub}$) around these lenses, finding this fraction to be broadly consistent with predictions from CDM although with large uncertainties. This was an important proof of method. Progress requires a larger sample of gravitational lenses which enables a measurement of not only the fraction of mass in substructure but also the slope of the subhalo mass function. Furthermore, Xu et al. (2015), Hsueh et al. (2016), and Gilman et al. (2017) have demonstrated that systematic uncertainties may occur in flux-ratio measurements if the deflector is not accurately modelled due to insufficiently deep optical imaging, as may be the case for studies which rely entirely on radio imaging.

There are several paths forward for increasing the sample of lenses which can be used to measure the subhalo mass function. Gravitational imaging for instance, can be used to detect subhalos perturbing the lensed positions of background galaxies. This method currently has a limiting mass sensitivity of $\sim 10^{8}M_{\odot}$ (Vegetti et al., 2012, 2014; Hezaveh et al., 2016; Birrer et al., 2017). The sensitivity limit is determined in part by imaging spatial resolution, which enables the measurement of the small astrometric perturbations to the strongly lensed backgrounds source caused by substructure. The next generation of telescopes and adaptive optics will lower this limit for background galaxy sources, and deep VLBI imaging of extended radio jets will potentially enable the detection of masses as low as $\sim 10^{6}M_{\odot}$ with gravitational imaging.

Observations of quasar lenses at longer wavelengths can provide a microlensing-free probe of substructure. Redward of $\sim 4\mu$m rest-frame, the quasar continuum emission is expected to be dominated by the dusty torus, which is sufficiently large to be unlensed by stars (e.g., Sluse et al., 2012). For a typical source redshift, this implies that mid-IR imaging at wavelengths greater than $10\mu$m can probe dark matter substructure. This method has been applied successfully to several systems with larger image separations (e.g., MacLeod et al., 2009; Minezaki et al., 2009; Chiba et al., 2005). JWST can provide the spatial resolution required to extend mid-IR flux ratio measurements to systems with smaller image separations. Jackson et al. (2015) demonstrated that deep radio observations of radio-quiet lensed quasars can also successfully be used as a microlensing free probe of substructure, albeit with larger flux uncertainties than radio-loud systems ($\sim 8-10\%$ compared with $\sim 3-5\%$ for radio loud systems). They estimate that this method can be applied to approximately half of optically selected quasar lenses.

Strongly-lensed quasar narrow-line emission provides an alternate probe of substructure, with comparable precision to radio-loud lensing studies. This method, originally proposed by Moustakas & Metcalf (2003), is extremely promising as it enables the measurement of substructure with current observational facilities in virtually all of the tens of optically selected quadruple quasar lenses predicted to be found in DES (Ostrovski et al., 2017; Agnello et al., 2015, http://strides.astro.ucla.edu/) and other wide field imaging surveys, and the hundreds forecast to be found in LSST (Oguri & Marshall, 2010). Sugai et al. (2007) demonstrated this method for the gravitational lens RXS 1131-1231 using seeing-limited observations with the Subaru integral field spectrograph. This lens has an unusually large separation of $>1$ arcsecond between each of the images.

Higher resolution imaging is necessary for the majority of quad lens systems, which often have at least one pair of images separated only by a few tenths of an arcsecond. Adaptive optics can provide the necessary spatial resolution to isolate the lensed images. Nierenberg et al. (2014, hereafter N14) used the integral field spectrograph, OSIRIS (Larkin et al., 2006), at Keck with adaptive optics to measure spatially resolved narrow-line flux ratios in the gravitational lens B1422+231 (Patnaik et al., 1992). Adaptive optics is an effective tool, however it can only be applied to those systems in which the narrow-line emission of interest falls in a suitable wavelength range for adaptive optics correction, for instance either $H$ or $K$ band in the case of Keck OSIRIS. Furthermore, adaptive-optics requires the presence of a nearby, bright, tip-tilt star, although often the lensed quasar itself is bright enough for this purpose. Space-based spatially resolved spectroscopy provides an alternative for those systems which fall outside of these wavelength bands, or are at a declination unsuitable to Keck OSIRIS spectroscopy.

In this work we demonstrate that the Hubble Space Telescope infrared grism on the Wide Field Camera 3 (WFC3) can be used to measure spatially-resolved narrow-line image fluxes, with comparable sensitivity to substructure to ground-based results from Keck. We present an analysis of WFC3 grism observations of HE0435-1223 (hereafter HE0435, Wisotzki et al., 2002) to demonstrate our reduction mechanism. The deflector is an early-type galaxy at redshift 0.4546 (Morgan et al., 2005) and the source is a quasar at redshift 1.693 (Sluse et al., 2012). This system has been extensively studied since its discovery; it has been monitored over a decade, and its relatively long time delay and significant intrinsic variability have made it a powerful probe of the group density profile and the Hubble constant (Wong et al., 2017; Bonvin et al., 2017; Courbin et al., 2011; Kochanek et al., 2006). Thanks to this attention, there are numerous multi-band and spectroscopic measurements available for comparison and significant effort has gone into measuring the properties of the environment of the deflector (Momcheva et al., 2006, 2015; Wong et al., 2011; Sluse et al., 2016; Wilson et al., 2016, e.g.), which is important for comparing detected subhalo properties with predictions from CDM.

In Section 2 we describe the observations and initial data reduction for this system. In Section 3 we describe the spectral extraction pipeline developed to measure narrow-line fluxes. In Section 4 we report measured quasar spectral features, and integrated emission fluxes, and compare the measured fluxes to results from broad band and radio studies. In Section 5 we perform a simple gravitational lens inference to test for the presence of substructure. In Section 6, we test for the effects of resolved narrow-line emission on our results. In Section 7 we discuss the constraint from this system. In Section 8 we provide a brief summary of the main conclusions.

We assume a flat $\Lambda$CDM cosmology with $h=0.7$ and $\Omega_{\rm m}=0.3$. All magnitudes are given in the AB system (Oke, 1974).

2 Observations and Initial Reduction

We observed the gravitationally lensed quasar HE0435 as a part of HST-GO-13732 (P. I. Nierenberg), a grism survey of narrow-line emission in six quasar lenses. The target was observed on August 30, 2015 for 2062s with the G141 Grism, and for 400s with F140W direct imaging. The observation was taken at a dispersion angle of 147 degrees East of North so that the dispersed quasar images would be maximally separated from each other along the direction perpendicular to the dispersion axis. In order to recover sub-pixel information, we split the observations into a four point dither pattern with half-integer sub-pixel offsets following the procedure of Brammer et al. (2012) (see also Schmidt et al., 2014; Momcheva et al., 2016). For each dither position, we took a 100s direct exposure with F140W, immediately followed by a 515s G141 exposure. The F140W direct images at each dither position were used to obtain accurate wavelength solutions for each G141 exposure (Brammer et al., 2012).

Raw F140W and G141 exposures were individually processed with AstroDrizzle (Gonzaga & et al., 2012) in order to reject cosmic rays, remove geometric distortion and perform flat-field subtraction (Koekemoer et al., 2011; Brammer et al., 2012). The F140W exposures were drizzled onto a 0$\aas@@fstack{\prime\prime}$06 pixel scale (approximately half the native pixel size). The upper left panel of Figure 1 shows the drizzled F140W image of the gravitational lens, and nearby spiral galaxy G1.

The G141 exposures were interlaced and combined onto a 0$\aas@@fstack{\prime\prime}$06 pixel scale, corresponding to an observed wavelength resolution of $\sim 22$ Å  per pixel following Brammer et al. (2012), Schmidt et al. (2014) and Momcheva et al. (2016). Unlike drizzling, interlacing does not introduce correlated pixel flux errors. Figure 1 shows the final interlaced grism data with arrows indicating the 5007 and 4959 Å [OIII] doublet which is partially blended given the grism resolution.

3 Spectral Extraction

Our goal is to measure the narrow-line emission flux in each lensed quasar image, taking into account blending between distinct spectral components after they are dispersed by the grism. The lower left panel of Figure 1 highlights the blending by showing how the light from the ring, quasar images and galaxies combine in a model grism image. In order to rigorously account for the overlapping spectra in the grism image, we employ a forward modelling approach. We discuss this method in detail in the following subsections; in brief, we generate a model direct image and use the 3D-HST grism simulation code (Brammer et al., 2012) to iteratively map proposed component spectra into a simulated 2D grism image, which is then compared with the original interlaced grism image to compute a $\chi^{2}$ goodness of fit. In Section 3.1 we discuss how the model direct image is generated, in 3.2 we discuss the 1D models we use for the spectral components, and in 3.3 we describe the statistical inference.

3.1 Direct image model

The grism image is effectively a convolution of an object spectrum with its direct image. Thus given a model for the direct image, it is possible to generate a predicted grism spectrum. We model the direct image as having seven distinct spatial plus spectral components: Four separate quasars, the main deflector, the quasar host galaxy which is lensed into a ring, and the nearby galaxy G1. Here we discuss how the direct image model is generated for each of the components. These direct model images are then combined with model 1D spectra, as discussed in Section 3.2, in order to generate model grism images.

We model the four quasar images as point sources, using a nearby star to model the point spread function (PSF). We optimize the point source positions and fluxes using galfit (Peng et al., 2002, 2010). A possible concern for using a drizzled star image as the model for the point source is that the true PSF is not accurately captured at the exact location of the lensed images. Furthermore, the FWHM of the grism PSF varies slightly with wavelength, which the 3D-HST pipeline does not account for in the forward modelling process. We have checked that the exact PSF model does not impact our inference on the [OIII] flux ratios by running the entire modelling process described in the next two subsections with a total of five different PSF models: 1) A median combination of stars in the F140W FOV, 2) and 3) The median star blurred by 10 and 15%, 4) and 5) two different nearby stars. While the choice of PSF model affected the overall fit to the 2D grism image, we found that it had no impact on our inference of the relative [OIII] fluxes. The quasar image positions are listed in Table A1, with uncertainties given by the variation in best fit galfit positions for the different PSF models.

In order to disentangle the lens galaxy light from the prominent ring and bright quasar images, we start with the empirical model of the deflector by Wong et al. (2017) derived from very deep (9337 s), F160W imaging. This model is generated from a superposition of Chameleon profiles (Dutton et al., 2011), and is based on a simultaneous fitting of the quasar point images, a model for the lensed background quasar host galaxy, and the lens. To generate a model for the F140W light profile, we start by fitting two Sérsic profiles to the empirical F160W model. Next, we hold all of the parameters for the Sérsic profiles obtained in the previous step fixed except for the total flux, and fit the F140W direct image as a combination of the lens and quasar images. We subtract the best fit galaxy and QSO models from the direct F140W image and are left with a residual image composed primarily of the ring. In the third step, we subtract the residual ring image from the original F140W image, and re-fit the ring-subtracted image with the galaxy and QSO models, this time allowing all of the Sérsic parameters for the galaxy model to vary freely. The final model for the galaxy is taken from the third step. We have tested several different iterations of this process, including allowing for less flexibility in the galaxy model to verify that the inferred narrow-line fluxes and image positions are not sensitive to the exact galaxy model, although the overall fit to the 2D grism image varies.

The ring model is generated by subtracting the best fit lens galaxy and quasar models from the direct image. We mask a small region near the centre of each QSO image where the PSF subtraction is noisy. We have confirmed that the inferred [OIII] fluxes are not sensitive to the exact size of the masked region. Finally for G1, we simply use a small cutout of the direct image, which is possible because it is isolated in the direct image.

3.2 1D Spectral Models

In this section we describe the analytic models we use for the 1D spectra. Although our grism data extend over a significantly larger wavelength region, we confine our model and comparison with the data to a small wavelength region around [OIII]. Figure 2 shows the fitted regions for each QSO image. The extent of the region is approximately rest frame 4500-5200 Å, but it varies between the lensed QSO images and is selected to achieve two goals. First, to provide sufficient spectral coverage to obtain a good constraint on the broad Fe, H$\beta$ and continuum features which overlap with the [OIII] emission (Figure 2). Second, the modeled region is extended where there is a possibility of an emission feature blending with the [OIII] emission of a neighbouring 2D spectrum. An example of the latter case is found in the redward extension of the fitted region for QSO spectrum C, in order to include the possible broad Fe flux contribution from image C to the image B narrow line emission (Figure 1).

Although gravitational lensing is not wavelength-dependent, the spectra of the four QSO images may not necessarily be related by a simple multiplicative magnification factor owing to the intrinsic variability of the QSO (which can lead to different image spectra owing to the varying image arrival times), and to the differential effect of stellar microlensing as a function of intrinsic source size. Thus we construct a model for the four QSO spectra which enables variations due to these effects. The QSO spectrum in the wavelength range of interest is composed of broad Fe and H$\beta$ emission, continuum emission and narrow [OIII] and H$\beta$ emission.

We model the broad H$\beta$ emission as a Gaussian with an intrinsic redshift offset from the [OIII] emission left as a free parameter. Although other studies have found that the H$\beta$ line profile can have a complicated structure, the low resolution of the grism data does not warrant a higher order model; Figure 2 demonstrates that a Gaussian is sufficient to fit the line profile in this case. The broad Fe emission is modelled using IZw1 templates which we interpolate in velocity space following Bennert et al. (2011) and Woo et al. (2006). The broad Fe emission velocity is independent of the broad emission width in our model. The model allows for a redshift offset between the broad Fe emission and the [OIII] emission.

We model the continuum emission as a straight line rather than a power-law, given the small fractional size of the wavelength region of interest. The slope and amplitude of the line are left as free parameters. The broad-line and continuum emission both arise from regions which can be affected by microlensing, which we discuss in more detail in Section 4. To account for this possibility we allow the width and amplitude of the H$\beta$ emission and the slope and amplitude of the continuum emission to vary independently between the model spectra. We did not find significant evidence for variations in the broad Fe velocities between images, and so kept them fixed for our final analysis, however we allowed the Fe amplitudes to vary independently from the H$\beta$ amplitudes.

Unlike the broad and continuum emission, narrow [OIII] and H$\beta$ emission come from a sufficiently extended source (greater than tens of parsecs) to not be affected by either stellar microlensing or intrinsic variability (Moustakas & Metcalf, 2003; Müller-Sánchez et al., 2011; Bennert et al., 2006b, a). Owing to this, we assume that both the line widths and the relative amplitudes of [OIII] and narrow H$\beta$ should be constant between the lensed images. We model the [OIII] doublet and H$\beta$ narrow-lines as Gaussians, and assume that they have the same redshift, which is valid given the spectral resolution of the grism. The ratio of the [OIII] doublet 4959 and 5007 amplitudes is fixed to the quantum-mechanically predicted value of 1/3.

The 1D models for the deflector, ring and G1 spectra are modelled as straight lines over the short wavelength region of interest, with amplitudes and slopes as free parameters. We do not find evidence requiring the inclusion of emission or absorption features in any of these spectra relative to the measurement uncertainties and given the brightness of the QSO spectra (see e.g. Figure 2).

We assume that the image fluxes are not affected by differential dust extinction. In the rest frame of the lens, the [OIII] emission lines lie at roughly $\sim 9300$ Å. At this wavelength, total dust extinction in lens galaxies, and early-type galaxies in general, is typically of order only a few hudredths of a magnitude (e.g Falco et al., 1999; Ferrari et al., 1999), which is well within our overall flux measurement uncertainty. This assumption is further supported by the similarity of the broad-band optical colours of the images (Wisotzki et al., 2003). The images also have mutually consistent CIV (lens rest frame $\sim 2790$ Å) and H$\beta$ (lens rest frame $\sim 9300$ Å) broad-line flux ratios.

3.3 Inference of QSO spectral parameters

We infer the probability distribution of the parameters of the 1D spectral models using a Bayesian forward modelling approach with the emcee Markov Chain Monte Carlo software package (Foreman-Mackey et al., 2013). For each step, the MCMC algorithm proposes parameters for the 1D spectra of all seven distinct spectral components (four QSO images, the main galaxy, the lens ring and G1). We then simulate dispersed images of each separate component and add them to generate a full model 2D grism image. Finally, the $\chi^{2}$ of the fit is computed relative to the original 2D interlaced image. Figure 1 illustrates how the model 2D direct image components are dispersed into the model 2D grism image for each MCMC step.

4 Spectral Forward Modelling Results

Figure 2 shows the 1D model, data and residual ‘traces’ for the four lensed QSO images. These traces are obtained by integrating the flux along the y axis in the 2D image, weighted by the relative flux of the direct F140 model PSF along that axis. Jumps in flux are due to small misalignments between the dispersion axis and the detector axis. This comparison shows that the input model provides an excellent fit to the observed spectra.

From the spectral modelling we obtain flux ratios between the broad H$\beta$ fluxes and the [OIII] fluxes from the image pairs A/C, B/C, D/C. Given that the intrinsic quasar luminosity is not known, gravitational lensing analyses rely on ratios of image fluxes rather than their absolute values. In Figure 3, we compare these flux ratios with measurements from other studies across a range of filters and for fixed filters at multiple dates. These measurements are chosen to represent how the flux ratios vary with wavelength and time, and are only a small subset of the many measurements of this system obtained for time variability studies (e.g. Bonvin et al., 2017; Courbin et al., 2011; Kochanek et al., 2006). Table A2 contains references and observing dates for all flux ratios plotted in Figure 3.

The narrow [OIII] flux ratios are strikingly different from optical to near-IR flux ratios which are subject to contamination by microlensing and intrinsic QSO time variability.

HE0435 has been monitored for 15 years (Bonvin et al., 2017; Courbin et al., 2011; Kochanek et al., 2006), and during that time has shown highly variable broad band flux ratios due to stellar microlensing and intrinsic variability. The intrinsic variability particularly affects images B and D which have time delays of over a week relative to images A and C. Figure 3 highlights several repeat measurements of the system which show significant variability.

Based on simulations of QSO accretion disks and dusty tori, blueward of rest-frame $\sim$4$\mu$m, (observed $\sim$10 $\mu$m) the accretion disk makes a dominant contribution to the QSO emission (Sluse et al., 2013). From chromatic microlensing studies of this system, the quasar continuum emission has a half light radius of $\sim 10^{16.3\pm 0.3}$ cm (or 0.003 pc) at a rest-frame wavelength of 8000 Å. Sluse et al. (2012) estimate the MgII broad line region size to be $\sim 10^{18\pm 1}$ cm (or 0.03 pc). These sizes correspond to $0.5-5\mu$as at the lens redshift, and are thus affected by stellar microlensing. Bonvin et al. (2017) have inferred the approximate amplitude of the observed $R$ band (rest frame $\sim 2500$Å) stellar microlensing as a function of time for each of the images since 2003. Figure 4 shows an estimate from Bonvin et al. (2017) for the microlensing effect on the $R$ band flux ratios as a function of time assuming a ‘true’ flux ratio value indicated by the straight lines. The amplitude of microlensing depends on the source size thus bluer filters are more strongly affected by microlensing while redder continuum measurements are less affected. Blackburne et al. (2014) have performed a detailed study of differential microlensing as a function of wavelength for this system and their data are included in Figure 3, and in Table A2.

The broad line emission flux ratios of both H$\beta$ from this study and CIII] and CIV from Wisotzki et al. (2003) are closer to the narrow-line emission flux ratios, which is consistent with microlensing being a function of emission region size.

We can test for the effects of microlensing on our data by comparing the relative amplitudes of emission features in our inferred spectra. In Figure 5, we plot the marginalised models for the lensed image spectra from our analysis, normalised to the peak of the [OIII] flux at 5007 Å  in order to highlight how the emission features vary between the lensed images. Image A shows significant morphological differences, with continuum and broad H$\beta$ fluxes which are much higher relative to the [OIII] flux than the other three images. This indicates that there is significant source-size dependent lensing. This finding is consistent with the inferred $R$ band microlensing of image A as observed by Bonvin et al. (2017) and shown in Figure 4.

The narrow-line [OIII] flux ratios are consistent with 5 GHz radio measurements from Jackson et al. (2015), with A/C, B/C and D/C differing at 0.25, 1.8 and 2.2 $\sigma$ respectively. This is expected given that both sources are expected to be extended enough to avoid all microlensing contamination. Although the results do not differ significantly, we note that Jackson et al. (2015) found that their radio emission was somewhat resolved, with an intrinsic source size of $\sigma\sim 288$ pc, assuming the source had a Gaussian flux distribution. This affects the flux ratios predicted from gravitational lensing relative to a point source for a fixed deflector mass model. We discuss this further in Section 6, where we also place limits on the size of the narrow-emission region in our data and we examine the effects of a resolved narrow emission line region on our results.

5 Gravitational lens modelling

The lensed image positions and [OIII] flux ratios are sensitive to the mass distribution of the deflector. As discussed in the Introduction, the image fluxes are particularly sensitive to small scale perturbations caused by dark matter subhalos. In this section we perform a gravitational lensing analysis of the system using the image positions and [OIII] fluxes. We do not include the time delays as a constraint, given that they are minimally sensitive to perturbers near the lensed images, relative to the image fluxes and positions, and given the time delay measurement uncertainties for this system (Keeton, 2009; Keeton & Moustakas, 2009). Furthermore, unlike image positions and fluxes, which are most sensitive to the local mass distribution and can thus be well matched without including G1 explicitly (e.g. Sluse et al., 2012), the time delays are sensitive to the larger scale environment and using them as a constraint would require including the complex multi-plane lensing effect of G1 which is at a higher redshift than the main deflector (Wong et al., 2017; Bonvin et al., 2017; Jackson et al., 2015; Sluse et al., 2016). The macromodel parameters including the external shear term are left free to vary, and thus absorb the large-scale contributions from G1.

In Subsection 5.1 we discuss the optimum smooth mass model fit to the data. In Subsection 5.2 we place limits on the presence of substructure near the lensed images. In Subsection 5.3 we discuss the effects of a finite narrow-line emission region.

5.1 Smooth Model

We start with a simple, smooth mass distribution model of the system. We model the main deflector as a singular isothermal ($\rho(r)\propto r^{-2}$) ellipsoid (SIE) which has been shown to provide an excellent match to the combined stellar and dark matter mass distributions of elliptical galaxies well beyond the Einstein radius (e.g. Rusin et al., 2003; Rusin & Kochanek, 2005; Gavazzi et al., 2007; Treu, 2010; Gilman et al., 2017). The SIE has five free parameters; the centroid, the ellipticity and position angle, and the Einstein radius. The smooth model also includes parameters describing the magnitude and direction of external shear which can be generated by the group environment of the galaxy.

We optimize these model parameters relative to the observed lensed image positions and [OIII] fluxes using gravlens (Keeton, 2001b, a), and find an overall best fit $\chi^{2}$ of 1.7 for one degree of freedom. Table 1 lists the mean and one sigma uncertainties for the lens model parameters. Figure 3 shows the best-fit model fluxes as straight lines. The $\chi^{2}$ of 1.7 for one degree of freedom indicates that the data has a roughly $\sim 20\%$ chance of being drawn from the best fit model, thus we do not find significant evidence indicating that a more complex model with a significant subhalo contribution is necessary to fit the data. Our inferred lens model parameters are comparable to other values in the literature which are inferred from microlensing-free data. In particular, our inferred Einstein radius and external shear parameters are consistent with the single lens model parameters from Sluse et al. (2012). As expected, our inferred spherical equivalent Einstein radius of 1$\aas@@fstack{\prime\prime}$200$\pm$0$\aas@@fstack{\prime\prime}$003 for the single lens model is somewhat higher than other values in the literature which explicitly include G1 and find values for the Einstein radius of the main deflector ranging from 1$\aas@@fstack{\prime\prime}$07$\pm$0$\aas@@fstack{\prime\prime}$02 (Jackson et al., 2015) to 1$\aas@@fstack{\prime\prime}$182$\pm$ 0$\aas@@fstack{\prime\prime}$002 arcseconds (Wong et al., 2017). We also infer a somewhat higher value for the external shear than models which include G1 in addition to an external shear; $\gamma_{ext}=$0.063$\pm$0.007 compared with 0.039$\pm 0.004$ and 0.030$\pm 0.003$ for Jackson et al. (2015) and Wong et al. (2017) respectively.

This result differs from the result of Fadely & Keeton (2012) who found that the $L$ band image fluxes and positions could not be fit with a smooth model although they did not report their best fit lens model parameters. They note that observed $L$ band corresponds to the rest frame 14000 Åat the redshift of the quasar, which should have a significant flux contribution from the dusty torus in addition to the accretion disk continuum emission (Sluse et al., 2013; Wittkowski et al., 2004; Hönig et al., 2008). In order to test for possible microlensing of the continuum component of the $L$ emission, they analysed two years worth of monitoring data from Kochanek et al. (2006). Unluckily over this time scale there was not significant evidence for microlensing induced variability in image A. Longer baseline data over 15 years (Bonvin et al., 2017) reveals significant microlensing of image A, as shown in Figure 4, which likely affected $L$ band flux ratios.

In the following subsection we place limits on the presence of perturbing subhalos near the lensed images.

5.2 Limits on the presence of substructure

Each of the image positions and fluxes provides a local constraint on the presence of small scale structure. In this subsection we test the limits on the presence of a single perturbing subhalo given our [OIII] flux and position measurements. We test the measurement sensitivity to two different perturber masses of $M_{600}\sim 10^{8}M_{\odot}$ and $10^{7}M_{\odot}$ where $M_{600}$ is the integrated mass within 600 pc of the centre of the perturber. These masses are chosen to be above and below the limit where the ‘Missing Satellite Problem’ is observed in the Milky Way (e.g. Strigari et al., 2008).

As demonstrated in Nierenberg et al. (2014), the perturber mass distribution can significantly affect the predicted lensing signal with fixed $M_{600}$. This is due to the fact that shallower mass profiles must have a higher overall normalisation than steeper mass profiles in order to achieve the same interior $M_{600}$ integrated mass. This in turn causes the shallower mass profile to have a longer range impact on the observed image fluxes and magnifications. Here we demonstrate the lensing effect for two mass profiles, a singular isothermal sphere (SIS) which has an $m(r)\propto r^{-2}$, and a Navarro, Frenk and White (NFW Navarro et al., 1996) halo, which has a shallow interior profile of $m(r)\propto r^{-1}$, which transitions to a steeper value $m(r)\propto r^{-3}$ outside of a scale radius. We obtain the scale radius from the mass-concentration relation predicted by Macciò et al. (2008) assuming a WMAP5 cosmology (Dunkley et al., 2009).

For each perturber mass, and for each mass profile, we iteratively place the perturber at a fixed position, re-optimize the smooth model parameters, and compare the new best fit $\chi^{2}$ with the original $\chi^{2}$ in the absence of a perturbation. We choose the grid spacing qualitatively to ensure that relevant angular dependences are captured in the $\chi^{2}$ distribution as a function of position. We find that variations are well captured by a spacing of 0$\aas@@fstack{\prime\prime}$1 for the $M_{600}\sim 10^{8}M_{\odot}$ mass perturber and 0$\aas@@fstack{\prime\prime}$01 for the $\sim 10^{7}M_{\odot}$ respectively.

Figure 6 shows the projected two and three sigma exclusion regions ($p<5$% and $0.3$%) for a singular isothermal sphere perturber with Einstein radius of 0$\aas@@fstack{\prime\prime}$01 and 0$\aas@@fstack{\prime\prime}$001 respectively. Assuming the perturber is in the plane of the lens galaxy, these Einstein radii correspond to integrated masses of $\sim 10^{8.2}$ and $10^{7.2}M_{\odot}$ respectively within 600 pc of the centre of the perturber, making them comparable to the Milky Way satellites Fornax and Sagittarius (Strigari et al., 2008), albeit with steeper mass profiles. Based on the average minimum radius at which including a perturber results in a model probability which is lower than $0.3\%$, the average exclusion radius is $\sim$0$\aas@@fstack{\prime\prime}$4 (0$\aas@@fstack{\prime\prime}$1), 0$\aas@@fstack{\prime\prime}$3 (0$\aas@@fstack{\prime\prime}$08), 0$\aas@@fstack{\prime\prime}$4 (0$\aas@@fstack{\prime\prime}$09) and 0$\aas@@fstack{\prime\prime}$3 (0$\aas@@fstack{\prime\prime}$06) for images A, B, C and D for the 0$\aas@@fstack{\prime\prime}$01 (0$\aas@@fstack{\prime\prime}$001) Einstein radius perturber.

Figure 7 shows the projected two and three sigma exclusion regions for an NFW perturber with scale radius of 1$\aas@@fstack{\prime\prime}$0 and 0$\aas@@fstack{\prime\prime}$1 corresponding to integrated masses of $\sim 10^{8}$ and $10^{7.2}M\odot$ within the central 600 pc of the perturber. Again, assuming a best fit model probability lower than $0.3\%$ be excluded, the average radius of exclusion is A: 1$\aas@@fstack{\prime\prime}$2 (0$\aas@@fstack{\prime\prime}$1), B: 0$\aas@@fstack{\prime\prime}$3 (0$\aas@@fstack{\prime\prime}$08), C: 1$\aas@@fstack{\prime\prime}$1 (0$\aas@@fstack{\prime\prime}$09), D: 0$\aas@@fstack{\prime\prime}$8 (0$\aas@@fstack{\prime\prime}$06) for images A, B, C and D for the 1$\aas@@fstack{\prime\prime}$0 (0$\aas@@fstack{\prime\prime}$1) perturber. As expected, the NFW perturber must be further from the lensed images than an SIS perturber with the same mass to avoid significantly perturbing their fluxes and positions relative to the best fitting smooth mass model.

6 Finite Source Effects

Until this point we have assumed that the narrow-line emission in HE0435 is unresolved with HST resolution, even after being strongly lensed and thus magnified. In this section, we explore the effect of a resolved emission region on our results. Jackson et al. (2015) found that the 5 GHz radio emission fluxes for HE0435 were best fit with a smooth deflector mass distribution and a resolved source with with an intrinsic size of $\sigma\sim 34$mas (288 pc) for a Gaussian profile. This model significantly improved the fit to their data relative to an unresolved emission model. Given this, it is important that we test whether the narrow-line emission is resolved in our data, and what the effects would be of incorrectly assuming the narrow-line emission was unresolved.

The narrow-line region is extended and in some cases has been observed to extend out several hundreds of parsecs (e.g., Bennert et al., 2006a). The narrow-line flux is not uniform, and can be dominated by the central tens of parsecs even for very luminous quasars (Müller-Sánchez et al., 2011). Sluse et al. (2007) demonstrated that the [OIII] emission in the gravitational lens RXS J1131$-$1231 is resolved in their data, and thus inferred a minimum size of $\sim 150$ pc. Nierenberg et al. (2014) found that the narrow [OIII] emission for B1422+231 (Patnaik et al., 1992) was marginally resolved with a dispersion of $\sim 10$ mas (50 pc).

We perform two tests here. First, we simulate narrow-line emission for three different source sizes and then redo the forward modelling inference using these simulated extended narrow-line images in place of the original point source model for the direct image to test whether this leads to an improved fit to the observed spectra. Secondly, we simulate gravitational lenses with three different narrow-line source sizes, and model them under the (in this case, incorrect) assumption that the narrow-emission is unresolved to test for the effect on the inferred flux ratios.

To test the model fit with a resolved narrow-line emission, we repeat the forward modelling inference described in Section 3 with an additional model component. As before, the QSO broad and continuum emission regions are modelled as being emitted by a point source. In order to account for a resolved narrow-line source, we generate a new model direct image for the narrow-line emission only. This extended emission model is generated assuming that the intrinsic source is a Gaussian with width $\sigma_{\rm{NL}}$. We use the best fit gravitational lens model inferred from the continuum image positions and [OIII] flux ratios in order to generate lensed extended emission models for the direct image. We generate models for $\sigma_{\rm{NL}}$ of 50 pc, 100 pc, and 288 pc. The latter source size is the best fit size for the 5 GHz radio emission for this system found by Jackson et al. (2015).

We then re-estimate the best fit 1D spectral parameters following the steps in Section 3 with two adjustments. First we use the new simulated extended source model as the direct image model for the narrow-line emission. The other QSO spectral components are modelled as being emitted by point sources as before. Secondly, because the simulated direct image for the [OIII] fluxes by definition fixes the relative image fluxes, the spectral model has only a single parameter for the overall normalization of the narrow-line flux. Thus there are three fewer model parameters than the point source model in which the [OIII] fluxes vary independently.

The left panels of Figure 8 show the best fit simulated grism images for the narrow-line components only of the three extended models and the point source model for comparison. Note that the narrow H-$\beta$ emission is also modelled as being extended as we assume it is emitted from the same region as the [OIII] emission. While the 50 and 100 pc models differ only marginally from the point source model, the 288 pc source size is clearly extended with HST resolution. The [OIII] doublet for images B and D in particular are nearly completely blended.

In the right panels of Figure 8 we compare the best fit model 1D trace to the data (analogous to Figure 2) for each of the source sizes. The $\chi^{2}$ comparison between the data and the simulated grism image grows progressively worse as the source size increases, with best fit $\chi^{2}$ values of 13228, 13368 and 13990 for 3909 DOF for the 50, 100 and 288 pc narrow-line regions, compared with 13112 for 3906 DOF for the point source model222The point source model has three extra degrees of freedom as the narrow-line fluxes are allowed to vary independently. Increasing the narrow-emission source size results in a decreasing best fit peak narrow-line flux, despite the fact that the intrinsic narrow-line emission width is a free model parameter. This is due to the significantly extended emission which cannot be well fit in 2D.

We can also test what the effect would be on the inferred flux ratios, if we use a point source model for the narrow-line flux when the flux is actually resolved. For this we considered only the 50 pc and 100 pc models as the 288 pc model provides a markedly worse fit to the data. We simulated a mock spectrum using the resolved narrow-line models for the narrow-line component, and then inferred the image fluxes under the assumption used elsewhere in this work, that the narrow-emission is unresolved in our data. The resulting flux ratios were: A/C $=0.97\pm 0.08$, B/C $=0.99\pm$0.07 and D/C $=0.65\pm 0.05$ for the 50 pc source, and A/C $=1.0\pm 0.1$, B/C $=1.0\pm$0.1 and D/C $=0.65\pm 0.08$ for the 100 pc source. Both results are consistent with the input fluxes used to generate the resolved narrow-line source mock lenses, indicating that our results would not be biased by a marginally resolved narrow-line region.

7 Discussion

We have demonstrated that the WFC3 IR Grism provides sufficient spatial and spectral resolution to precisely measure lensed narrow-line emission with similar precision to continuum emission studies while avoiding the effects of microlensing, variability and differential dust extinction. We measure the [OIII] flux ratios to be significantly offset from optical/near-IR continuum measurements of this system, particularly for image A. This is consistent with the results from long term monitoring for this system, which have indicated that image A is significantly affected by microlensing, which systematically affects smaller emission regions. The [OIII] flux ratios are consistent with the 5 GHz radio measurement from Jackson et al. (2015), which are also extended and thus not affected by microlensing.

In order to fit the measured F140W image positions and [OIII] fluxes, we rule out the presence of a perturbing NFW subhalo projected within roughly $\sim$1$\aas@@fstack{\prime\prime}$0 (0$\aas@@fstack{\prime\prime}$1) for a $M_{600}\sim 10^{8}$ $(10^{7})M_{\odot}$ perturber near the lensed images. At the redshift of the lens these angular sizes correspond to approximately $\sim 6$ and 0.6 kpc respectively. It is informative to compare these limits with an approximate prediction of the number of subhalos at these radii based on CDM models.

We perform a very basic estimate assuming that all possible perturbers are within the virial radius of the HE0435 group and neglecting line-of-sight associations. From Sluse et al. (2016), the HE0435 group has a virial mass of $\log[M_{200}/M_{\odot}]=13.7\pm 0.4$, where $M_{200}$ is the mass within the region of the halo which has a mean density 200 times the critical density.

We estimated the number of subhalos based on the results of Han et al. (2016), scaled to the virial mass of HE0435, assuming that approximately half of subhalos are destroyed through tidal interactions and merging following the recommendation of Han et al. (2016), and that an additional $\sim 30$% are destroyed by tidal interactions with a central baryonic potential (e.g. Garrison-Kimmel et al., 2017)

We assume that the subhalo mass profiles are NFW with the mass-concentration relation given by Macciò et al. (2008) scaled to the redshift of the group following the relation by Prada et al. (2012) and neglecting scatter. We fix $M_{600}$ for subhalos at infall, where $M_{600}$ is the mass within the interior 600 pc of the subhalo. We make a simplifying assumption that $M_{600}$ is not affected by tidal stripping after infall. This yields an estimated number of $\sim 250$ surviving subhalos in the group with masses greater than $M_{600}>10^{8}M_{\odot}$ and $\sim 19000$ subhalos with $M_{600}>10^{7}M_{\odot}$ within $R_{200}=540$kpc.

We examined two different spatial distributions for the subhalos. First, that the subhalo number density follows the mass distribution of the host halo everywhere, as predicted by pure CDM simulations in the absence of tidal stripping. Second, that the subhalo number density follows the mass distribution of the host halo except within the three-dimensional scale radius, where we assume all subhalos are destroyed. This mimics an extreme version of the impact of the disk seen by Garrison-Kimmel et al. (2017). Both spatial distributions are normalized to have the same number of subhalos. These two subhalo spatial distributions are chosen to bracket limits of the possible effect of the baryonic potential on our predicted detection rate.

In the first case, where the subhalo number density simply follows the NFW profile of the lens halo, we expect approximately $\sim 0.8$ total subhalos with $M_{600}>10^{8}M_{\odot}$ to be found within each $\sim 6$ kpc projected excluded region. The mass function and size of the sensitivity region scales so that we expect to find $\sim 1$ $M_{600}>10^{7}M_{\odot}$ subhalo within each of the four smaller $\sim 0.6$ kpc projected regions around the images where we are sensitive to these lower masses.

In the second case in which all subhalos are removed within the NFW scale radius of the host, we expect to detect approximately 0.08 (0.1) subhalos with $M_{600}>10^{8}M_{\odot}$ ($10^{7}M_{\odot}$) per exclusion region.

We also explored the effect of a WDM mass function, using the result from Schneider et al. (2013) to estimate the shape of the subhalo mass function at infall given a 3 keV/c2 thermal relic dark matter particle, which is consistent with Ly$\alpha$ forest measurements (e.g. Viel et al., 2009). In this case, almost no subhalos with masses $M_{600}<10^{8}M_{\odot}$ survive. Unlike the central baryonic potential, which is predicted to have a largely mass-independent effect on the subhalo population, WDM selectively destroys low mass subhalos. Promisingly, subhalos with $M_{600}>10^{8}M_{\odot}$ are all expected to contain a significant number of stars (Strigari et al., 2008) in order to match comparisons between the Milky Way satellite luminosity function and CDM predictions (e.g. Strigari et al., 2008). This mean that CDM can be tested by comparing the rate of detections of substructure in narrow-line lenses with the rate predicted by luminous satellite studies, or measured by gravitational imaging studies which are sensitive to $M_{600}>10^{8}M_{\odot}$ subhalos, after accounting for variations in host halo mass and sensitivity. This test can be performed in a sample of $\sim 20$ lenses in which we would expect to detect $\sim 10$ $M_{600}>10^{7}M_{\odot}$ subhalos in the case of a CDM subhalo mass function, even in the case in which all subhalos are destroyed within the NFW scale radius.

We emphasise however, that a true comparison requires a CDM model which incorporates effects such as tidal stripping, as well as a marginalization over possible halo orientations, masses and formation histories (e.g. Jiang & van den Bosch, 2016), and a range of baryonic physics implementations (e.g. Chua et al., 2016; Despali & Vegetti, 2016). We leave such an analysis to a future paper, in which we jointly infer the properties of the subhalo mass function given our sample of narrow-line gravitational lenses measured with OSIRIS and the WFC3 grism.

8 Summary

We present a forward modelling method which uses the 3D-HST pipeline (Brammer et al., 2012) to measure spectra in the presence of significant spatial blending for G141 grism data. We apply this method to infer the lensed narrow, broad and continuum fluxes of the images in HE0435-1223. 2. 2.

The narrow [OIII] flux ratios for HE0435 are consistent with radio measurements from Jackson et al. (2015), and are significantly different from other emission measures which are subject to contamination by microlensing and intrinsic QSO time variability. 3. 3.

We find that the [OIII] fluxes and image positions are well modelled with a simple gravitational lens model consisting of a singular isothermal ellipsoid for the main galaxy in the presence of external shear. 4. 4.

Our data strongly disfavours a perturber with mass greater than $M_{600}=10^{8.2}(10^{7.2})M_{\odot}$ within $\sim 1$ (0.1) arcsecond of the lensed images, where $M_{600}$ is the projected perturber mass within its central 600 pc (best fit model probability $<0.3$%). 5. 5.

This demonstration that WFC3 grism measurement of narrow-line lensed quasars can be used to detect low-mass $M_{600}\sim 10^{7}M_{\odot}$ subhalos is extremely promising for future constraints of dark matter given the large number of quadruply imaged quasar lenses to be discovered in optical surveys such as DES and LSST, and with the follow-up which will be enabled by JWST.

Acknowledgments

Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program #13732. Support for program #13732 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. T.T. thanks the Packard Foundation for generous support through a Packard Research Fellowship. A.M.N thanks the Center of Cosmology and AstroParticle Physics for support via a CCAPP postdoctoral fellowship. C.S.K. is supported by NSF grant AST-1515876.

We thank K. Wong for providing model constraints for the deflector and inferred lensed image fluxes based on F160W data. We thank V. Bonvin and F. Courbin for providing microlensing light curve constraints.

Appendix A Table of measurements of flux ratios for HE0435-1223

Here we provide the F140W lensed quasar image position as well as a selection of the many measurements of system chosen to illustrate how the flux ratios vary with wavelength and time, as shown in Figures 3 and 4.