Image Flux Ratios of Gravitationally Lensed HS 0810+2554 with High-resolution Infrared Imaging

We report near simultaneous imaging using LMIRCam on the LBTI of the quadruply imaged lensed quasar HS 0810+2554 at wavelengths of 2.16, 3.7, and 4.78 μm with a full width at half maximum spatial resolution of 0.″13, 0.″12, and 0.″15 respectively, comparable to Hubble Space Telescope optical imaging. In the z = 1.5 rest frame of the quasar, the observed wavelengths correspond to 0.86, 1.48, and 1.91 μm respectively. The two brightest images in the quad, A and B, are clearly resolved from each other with a separation of 0.″187. The flux ratio of these two images (A/B) trends from 1.79 to 1.23 at wavelengths from 2.16 to 4.78 μm. The trend in flux ratio is consistent with the 2.16 μm flux originating from a small sized accretion disk in the quasar that experiences only microlensing. The excess flux above the contribution from the accretion disk at the two longer wavelengths originates from a larger sized region that experiences no microlensing. A simple model employing multiplicative factors for image B due to stellar microlensing (m) and substructure millilensing (M) is presented. The result is tightly constrained to the product m × M = 1.79. Given the observational errors, the 60% probability contour for this product stretches from m = 2.6, M = 0.69 to m = 1.79, M = 1.0, where the later is consistent with microlensing only.


Introduction
High-resolution N-body simulations of the cold dark matter (CDM) structure formation of the universe predicts that there should be hundreds of dark matter subhalos with mass ∼10 4 -10 9 M e within a massive halo (∼10 12 M e ). The fact that only a few dozen dwarf galaxies have been observed within our own Milky Way (McConnachie 2012) has led to the so-called "missing satellite" problem (Klypin et al. 1999;Moore et al. 1999;Sawala et al. 2014). The inclusion of baryonic physics (Brooks & Zolotov 2014;Dutton et al. 2016;Bullock & Boylan-Kolchin 2017;Buck et al. 2018;Garrison-Kimmel et al. 2019) in models of the formation of galaxies has made progress on resolving this issue, but many aspects of dark matter substructure are still uncertain (Despali & Vegetti 2017;Gomer & Williams 2018;Despali et al. 2019;Hsueh et al. 2019).
While these subhalos would remain undetected due to their lack of visible matter, they should still be detectable through strong gravitational lensing (e.g., Kochanek & Dalal 2004). Gravitational lensing provides a unique method for detecting the presence or absence of substructure in a galaxy outside of the Local Group. Any perturbation in the lens potential, be it stars or dark matter substructure, will create perturbations in the observed flux from some of the lensed images, if it is sufficiently compact and close to the images. It has been shown by Mao & Schneider (1998), Metcalf & Madau (2001), and Chiba (2002) that the presence of substructure in the lens mass model can reproduce the observed flux ratios in multiply imaged lensed quasars, something that smooth lens models cannot do. Hsueh et al. (2019) use a sample of seven quadruply imaged lensing systems and assuming a CDM cosmology and contribution from low mass halos along the line of sight, they infer an average total mass fraction in substructure that is in rough agreement with the predictions from CDM hydrodynamical simulations. Their result is significantly different when compared to previous studies that did not include line-of-sight halos (Xu et al. 2015). Their data set is dominated by images with good flux ratios at radio wavelengths, but with less contribution to the analysis from midinfrared (MIR) observations which probe smaller scales. Several lensed systems, including HS 0810+2554, have images that are probably resolved at radio wavelengths and were dropped from their initial sample of 14. Nierenberg et al. (2019) report observations of 8 quads, separating out emission from the quasar narrow line emission region (which is not subject to microlensing) and conclude that simple, smooth halo mass distributions are inadequate, indicating the need for more complex mass distributions. Although HS 0810 was in their sample, the small angular distance between images A and B made the analysis by Nierenberg et al. (2019) difficult, and narrow line flux ratios were not reported. Dobler & Keeton (2006) point out that the angular size of the source and the size of the perturber must be on the same order of magnitude for the source to be appreciably magnified. MIR observations have the advantage that in the rest frame of the lensed quasar, they span a range in wavelengths that are sensitive to the emission of both the very small sized accretion disk in the quasar which is subject to the effects of both stellar microlensing and substructure millilensing, and the much larger dusty torus which is subject only to millilensing (Sluse et al. 2013). Therefore, if we compare the flux ratios in a lensed quasar at different wavelengths, corresponding to different source sizes, we can quantify the size and location of substructure in the lens galaxy. This approach was taken by Fadely & Keeton (2011, 2012 using K and L′ bands for six lenses with image separations >1″ and by Chiba et al. (2005) at 11 μm for two lenses with image separations ∼1.0 and 0.5″. The closer the two brightest images in a quad are together, the less influence larger scale variations in the halo potential (e.g., ellipticity) and the effects of time delay will have (Lemon et al. 2017). Thus, quad images with very close separation are preferred in the search for evidence of halo substructure.
In this paper we present observations of the quad lensed system HS 0810+2554 (Reimers et al. 2002;Chartas et al. 2014;Hartley et al. 2019;Jackson et al. 2015;Nierenberg et al. 2019) at wavelengths from 2 to 5μm. The background quasar is at a redshift of z=1.5 but the lensing galaxy is at an unknown redshift. The two brightest images are separated by only 0 19 and at a redshift of z=1.5, the observations span a range in rest-frame wavelength of 0.8-2.0μm. Over this wavelength range the flux arises purely from the accretion disk at the shorter wavelength, and from a combination of accretion disk and dusty torus at the two longer wavelengths (e.g., Kobayashi et al. 1993;Suganuma et al. 2006). . The focal plane scale was 0 0108 per pixel. Images were made in the Ks (2.16 μm), L′ (3.7 μm), and M (4.78 μm) filters and the natural seeing was 0 8. We also made a test observation along with a standard star on UT 2017 February 17 in the L′ filter (only) with natural seeing of 1 2. Filters and integration times are listed in Table 1. The AO secondary was operating in natural guide star mode for both epochs and had to lock onto the optical images of the two adjacent bright images in the quad together. This probably explains the slightly elongated images we obtained for the quad compared to the standard star.

Observations
For the 2017 observations in the L′ filter, several dither positions on the array were used, but most often with a separation of 4 5. In 2019 we imaged HS 0810 at two dither positions on the array, separated by 8″. Images at each dither position in the 2019 observations were taken in sets of 10, 100, and 200 integrations at Ks, L′, and M respectively. Slightly elongated images, probably due to the AO locking onto the blended image, were evident. The axis of elongation rotates between wavelengths due to the rotation of the sky on the focal plane over the length of time the observations were made. The images occupied a small portion of the array, minimizing errors in flat-fielding. Simple aperture photometry was used to determine the brightness of the well isolated sources C and D in the image, whereas an elliptical point-spread function (PSF) was used to separate the slightly blended images of A and B, as explained in Section 3. Formal statistical errors were computed using the measured rms fluctuations in the background, which was determined to be flat within the vicinity of sources C and D. For blended sources A and B, a small systematic error of 0.01 mag. was added to account for any potential residual uncertainties in the background level (see Section 3) due to the more complicated PSF fitting.
Due to time and weather constraints in 2019, we were unable to observe standard stars on the same night as HS 0810 to fix the fluxes at each wavelength on a photometric scale. We have examined standard star observations on previous nights at Ks, L′, and M and found that the relative response of the detector between filters, adjusted for airmass and stellar colors, is consistent to within a few percent. We used this information to build the relative fluxes between filters for our observations of HS 0810. In this way we can place the fluxes of the images in HS 0810 on a spectral energy distribution (SED; with a flux normalization) with confidence. This does not affect the relative brightness between images at each wavelength which are, within the photometric and systematic errors, very well determined.

Analysis
Grayscale images of a star (HD 74721) and HS 0810 from the 2017 test observations are shown in Figure 1. The stellar image has had dithers subtracted, but has not been processed   for bad pixel and pattern removal nor has the image been flat fielded. The image of HS 0810 has been fully processed. The individual lensed images are labeled according to Reimers et al. (2002). A radial profile of the stellar image is shown in Figure 2, where the measured full width at half maximum (FWHM) is 0 10. Grayscale images of the quad from the 2019 observations are shown in Figure 3. The lensing galaxy itself is barely visible in the Ks band, and finding its centroid is impossible using our data. Consequently we cannot measure the positions of the lensed images with respect to the lensing galaxy. The positions relative to image A are, within an error of ±0 003, identical to the values listed in the CfA-Arizona Space Telescope Lens Survey (CASTLES) website (Muñoz et al. 1998). The brightness for images C and D were computed using standard aperture photometry on the images. For images A and B, their individual PSF overlap. For these two images, the PSF was modeled with a simple two-dimensional Voigt profile with a slight ellipticity. The Gaussian width, dampening wings parameter, and orientation (long axis) were derived from image C at each wavelength. The measured values for the Gaussian width (1σ) were 0 055, 0 051, and 0 064 at Ks, L′, and M respectively. The dampening wings made a small, but detectable contribution to the overall PSF. The ellipticity was 15% at Ks and L′, and 5% at M.
A synthetic image of the blend was then formed using two model Voigt profiles with adjustable spacing and peak intensity. This model was subtracted from the image and the spacing and peak intensity of each profile adjusted until the residual intensity was less than 1% of the total intensity of the A+B combination. A cut through the centers of images A and B in the L′ filter for the 2019 observations is shown in Figure 4 along with the model fit. A slight tilt to the baseline background can be seen in the image cut that is not modeled. Since our PSF modeling is more complicated than the simple aperature photometry used for images C and D, a systematic error of 0.01 mag was added to the error budget to account for uncertainties in the method. The rest of the error budget was determined by comparing separate sets of observations that were coadded to make the final images and computing the statistical variation between these subsets. The results are given in Table 2, with the ratios of images B, C, and D to A expressed as magnitude differences, e.g., for images B and A, Δm=−2.5log(B/A). Figure 5 shows images A and B with scaled contours overlaying the grayscale image at each wavelength. In all three filters image A, which is the minimum in the lensing arrival time surface, is brighter than image B, which is the saddle point. This is consistent with these images being affected by substructure-stars or subhalos, where the saddle is expected to be demagnified with respect to the neighboring minimum (Schechter & Wambsganss 2002;Saha & Williams 2011). Clearly visible is the trend toward more equal brightness with increasing wavelength. Our goal is to determine if the ratio of B/A as a function of wavelength can be used to determine if there is substructure in the dark matter halo of the lensing galaxy or if all of the B/A ratios are consistent with microlensing only. Because the two images are very close on the sky in comparison to the scale of the quad system, their light paths are traveling through nearly identical regions in the gravitational potential of the lensing galaxy. If there were no microlensing and the dark matter halo associated with the lensing galaxy was smooth azimuthally in the region of images A and B, images A and B should have identical fluxes.
We estimated the time delay between images A and B using a free-form lens modeling software PixeLens (Williams & Saha 2004), which is dependent on the (unknown) redshift of the lensing galaxy. Mosquera & Kochanek (2011) made an estimate of z=0.89 for the lensing system based on basic photometric properties of the galaxy. For z=0.3 (z = 0.6, z = 0.89), the time delay is ∼0.05 day (∼0.15 day, ∼0.33 day) respectively. These delays are much shorter than typical intrinsic quasar variability (MacLeod et al. 2012), and comparable to the ∼3 hr of telescope time we spent making the observations. For all practical purposes, our observations of the A and B quasar images can be considered simultaneous in the rest frame.
Our flux ratios, expressed as magnitude differences, are also listed in Table 3, which includes ratios at other wavelengths from Chartas et al. (2016), Jackson et al. (2015), and Reimers et al. (2002), along the epoch of observation. We have been unable to find information on when the CASTLES observations cited in Chartas et al. (2016) were made, so no epoch is given in the Table. Note that there is clear variation in the short wavelength (X-ray-Ks band) B/A ratio over time, ranging from Δm=0.2 in the X-rays in 2002, to 0.635 at Ks in 2019. According to Mosquera & Kochanek (2011), a typical timescale for a microlensing event in the HS 0810 system due to crossing the Einstein radius of a single star is ∼24 yr, a timescale comparable to the elapsed time between the earliest and most recent observations. For a typical projected density of stars in the lensing galaxy, microlensing will likely result from the combined effect of many stars along the line of sight, rather than a single star. Modeling by Wambsganss (1990) shows that microlensing is a nonlinear process, and both a longer term trend in magnification in addition to variations on timescales shorter than the crossing time, combine to create the light curve of an image. This is similar to what we observe for the B/A ratio, which is always in the same sense (A > B), but with significant variations over the last 20 yr.
To establish the range in mass and size of halo substructure our observations are sensitive to, we list the relevant source plane sizes in Table 4. The values for the stellar Einstein radius R E , the accretion disk radius R AD , and the radius of the broad line region R BLR , were taken from Mosquera & Kochanek (2011). We estimated the size of the outer radius of the dusty torus from Table 1 in Sluse et al. (2013) assuming HS 0810 has a luminosity of L∼5×10 45 ergs s −1 and that the torus is represented by the "extended" model. The luminosity estimate was made using the black hole mass given in Mosquera & Kochanek (2011) and the mass-luminosity relation in Woo & Urry (2002). We chose the extended torus model since the infrared flux from HS 0810 is greater than for the typical quasar (Figure 6), suggesting a more extensive emitting region.
With these size estimates we can roughly estimate the mass range for any halo substructure that can significantly influence the flux ratios we measure by using a Sérsic radial profile to compute an Einstein radius. Using a suite of numerical simulations and fitting radial profiles with an Einasto function, Navarro et al. (2004) derive a shape parameter of ∼0.17 that fits the simulated halos well. Using Figure 8 in Dhar & Williams (2010), this approximately corresponds to a Sérsic profile index of m S ∼6. We use a Sérsic profile because the Einstein radius can be analytically calculated for a given index and half-mass-radius R 1/2 . For a mass of 10 4 M e and assuming the lensing galaxy is at z=0.89, the Einstein radius ranges from 1.1 pc for R 1/2 =0 (point mass) to 0.5 pc for R 1/2 =2 pc projected onto the source plane. These sizes are comparable to the size of the infrared torus and represents a rough lower limit to the mass we can detect. The distance between images A and B in the lensing plane is ∼1.45 kpc, again assuming the lensing galaxy is at z=0.89. This corresponds to an Einstein radius for a point mass of 2.1×10 10 M e , or 1×10 10 M e for R 1/2 =2 kpc. This represents a rough upper limit to the mass that can influence our observations. The mean spectral energy distribution for quasars has been computed by Padovani et al. (2017) and Krawczyk et al. (2015). They find that the combination of a small, subparsec sized accretion disk with a simple power-law SED well matches the rest-frame optical flux. At longer wavelengths, starting at about 1.2μm, flux from a much larger (few parsecs) dusty torus begins to contribute and eventually dominates the SED. Krawczyk et al. (2015) found a mean power-law index for the accretion disk of F ν ∝ ν −0.4 . The dust in the torus is at a temperature of ∼1500 K (e.g., Kobayashi et al. 1993;Suganuma et al. 2006), and the flux from the dusty torus drops very rapidly shortward of 1.2μm, which is on the Wien side of the emission spectrum for the dust. Our 3.7 and 2.16μm filters correspond to wavelengths of 1.48 and 0.86μm in the rest frame of the lensed quasar. At these rest wavelengths, the flux from a T=1500 K blackbody in νF ν at the shorter wavelength is a factor of 12 lower than at the longer wavelength. Given that the dust emission in our 3.7μm filter is only about 40% (see below) of the total, the contribution to our 2.16μm filter from hot dust emission is only ∼3%, and can be ignored.  We plot the SED for the entire (within a single photometric aperture) HS 0810 quad in Figure 6, along with the mean SED for quasars taken from Krawczyk et al. (2015). The flux values plotted for HS 0810 were taken from the WISE, 2MASS (Skrutskie et al. 2006), and UKDISS (Hambly et al. 2008) catalogs, as well as a quasar catalog by Souchay et al. (2015). As noted above, the K filter is likely dominated by the small accretion disk, and consequently could show considerable variation in flux over time due to stellar microlensing. For this reason we plot the observed J and K band fluxes as a vertical bar, indicating the range of fluxes between these data sets. For both the 2MASS and UKIDSS surveys, the J, H, and K observations were made at the same time, so the J-K colors should represent the quasar accretion disk continuum at the epoch of observation. The 2MASS and UKIDSS H filter observations are strongly contaminated by Hα and do not represent the continuum. The J band observations suggest a somewhat steeper (bluer) power-law index of −0.3, but this filter could be contaminated by [O III], although the effect should be small.
We only have relative flux levels between our infrared filters and we have no simultaneous red and optical photometry. This makes it difficult to empirically determine the power-law slope for the accretion disk in HS 0810 at the epoch of our observations. Based on the J-K colors from 2MASS and UKIDSS, we will fix the power-law slope at −0.3 as shown in Figure 7. We have vertically scaled our relative flux measurements for the A+B pair to match the trend between the WISE W1 and W2 filters. Flux from images A and B by far dominate our images and would be expected to dominate the overall SED, and we will assume the shape of the overall SED represents the shape of the SED for the combination of images A and B. With this normalization, our Ks (2.16 μm) point fits within the lower range of previous K band observations, although this is not a given.

Lensing Model
In Figure 7 we plot the SED for HS 0810 where the dashed blue line shows a simple accretion disk power law with n µ n -F 0.3 that passes through our Ks flux. Working with this figure, we developed a simple model for the flux from images A and B by making the following assumptions.
(1) The Ks flux is due entirely to a small sized (fraction of a parsec) accretion disk that can suffer microlensing in the lensing galaxy.
(2) The longer wavelength flux consists of a contribution from the accretion disk (extrapolated from the Ks flux) and a contribution from a larger dusty torus that does not experience any microlensing. Expressed mathematically, we have: I λ is the total flux from A and B, AD indicates flux from the accretion disk alone, and IR represents the flux from the torus alone, as indicated in Figure 7. The factor m is the ratio of I A /I B due to microlensing. If lensing due to substructure in the dark matter halo was present in addition to microlensing, an additional magnification factor M must be added. This factor, which we will call millilensing, effects both the emission from the accretion disk and the infrared emission from the torus. Observations of HS 0810 at radio wavelengths (Jackson et al. 2015) and in optical narrow emission lines (Nierenberg et al. 2019), neither of which should be subject to microlensing, are consistent with no millilensing. For this reason we initially set this value to M=1.0. Given the observed SED and our assumptions listed above, we can solve for the factor m that provides a best fit to the data at all wavelengths and compare to the observations, where: The results of our model fit for the case with only microlensing, no millilensing (M = 1), are shown in Figure 8, where we have expressed the flux ratios as a magnitude difference for comparison with the data (Section 3). The data Figure 5. Contours for images A and B arbitrarily scaled to produce 11 contours from sky to peak brightness. Note the clear trend toward more equal brightness with increasing wavelength.  Table 2 along with the observed differences. Our simple model fits the data well and is consistent with A/B ratios at radio wavelengths (Jackson et al. 2015;Hartley et al. 2019) of 1.02±0.06 and model results using narrow line emission (Nierenberg et al. 2019) that imply a flux ratio of 1.0±0.07. If substructure were present at a level that could significantly influence the observed flux ratios, then the measured ratios for B/A at the two longer wavelengths would lie above or below the simple model that is compatible with the flux ratio at Ks. In this case, the product m×M determines the flux ratio at Ks and the inclusion of the millilensing factor changes the flux ratios at L′ and M from what is observed. In Figure 9 we plot the probability contours for our model in the m, M plane. The red point is the case for microlensing only, no millilensing. The contours take on the shape of a very narrow band and there is a significant degeneracy along the line m×M=1.79. The extent of this degeneracy can be reduced only by the inclusion of longer wavelengths where the contribution of the accretion disk becomes insignificant. The small errors on our flux ratios make the contours very tight perpendicular to the m×M=1.79 line. Changing the powerlaw slope we used for the accretion disk moves the contours along the line m×M=1.79 where steepening the power-law index (bluer) will move the contours up to the left, and vice verse.
The contribution for HS 0810 from the accretion disk in our simple model continues to be a factor even at longer wavelengths (see Sluse et al. 2013). Our model predicts a flux ratio of B/A=0.9/1 will not occur until wavelengths longer than ∼10μm. The FWHM of the LBT at this wavelength is 0 3, but with a high signal-to-noise observation it is still possible to extract an accurate flux ratio given the very stable telescope PSF at this wavelength and the known separation of the two point sources. The source brightness at 10μm is ∼30 mJy, which is sufficient given LBTIs high sensitivity at infrared wavelengths. Although the upcoming JWST will be able to image much deeper at thermal MIR wavelengths than ground-based telescopes such as the LBT, with a smaller 6 m primary it will be difficult to resolve adjacent images in quad systems such at HS 0810. Davidson et al. (2019) report detection of all four images at a wavelength of 2.1 mm in the continuum and in CO line emission using ALMA. At this long wavelength the emitting region is likely considerably larger than the hot, inner dusty torus we observe in the MIR, and the  Krawczyk et al. (2015). Red open circles: combined flux in our filters scaled vertically to match the trend between the WISE W1 and W2 filters. The vertical bars indicate the range in fluxes between 2MASS and UKIDSS J and K band observations. The H band (1.65 μm) flux is not plotted due to strong contamination by Hα. Figure 7. Similar to Figure 6 except that a power-law fit (dashed blue line) to the short wavelength portion is used to extrapolate the flux contribution from the physically small accretion disk to longer wavelengths. This power-law fit is slightly steeper than the mean power-law index shown in Figure 6 (see the text). This yields the "AD" terms in Equation (1). The vertical dashed lines indicate the location of our three filter bandpasses, and the difference between the accretion disk SED and the observed flux gives the "IR" terms in Equation (1). images may be insensitive to millilensing on scales smaller than a few parsecs, but capable of detecting substructure on larger scales. MIR imaging on the next generation of very large, ground-based telescopes equipped with an infrared AO system presents a future opportunity to further constrain the presence of millilensing due to substructure in dark matter halos. Alternatively to long wavelength imaging, observations of emission from the quasar narrow line region, also unaffected by microlensing, show promise (e.g., Nierenberg et al. 2019), but also require high angular resolution.

Conclusions
We have imaged the lensed quasar system HS 0810+2554 at wavelengths spanning 2-5 μm with LMIRCam on the LBTI using the adaptive optics system. Flux ratios were computed for all four images in each filter. The two brightest images are 0 187 apart, and have very high signal-to-noise at all wavelengths. At the shorter wavelength of 2.16 μm (rest wavelength 0.86 μm) the flux is entirely consistent with flux from only the subparsec accretion disk of the quasar. At the two longer wavelengths, the flux is a combination of flux from the accretion disk and flux from a larger (few parsec) dusty torus. A simple model employing multiplicative factors for image B due to stellar microlensing (m) and substructure millilensing (M) was presented. The result is tightly constrained to the product m×M=1.79. The 60% probability contour for this product stretches from m=2.6, M=0.69 to m=1.79, M=1.0, where the latter is consistent with microlensing only. Note that we have not ruled out the presence of millilensing, but have shown that the microlensing only case is within the 60% probability contour.