DISSECTING THE GRAVITATIONAL LENS B1608+656. I. LENS POTENTIAL RECONSTRUCTION^*

The first step in solving Equation (3) for the potential perturbation is to obtain the source-intensity gradients and the intensity deficit, which appear in the correction equation. We follow Suyu et al. (2006) to obtain the source-intensity distribution on a grid of pixels given the current iteration's lens potential model. In this source reconstruction approach, the data (observed image) are described by the vector d_j, where j = 1, ..., N_d and N_d is the number of data pixels. The source intensity is described by the vector s_i, where i = 1, ..., N_s and N_s is the number of source-intensity pixels. The observed image is related to the source intensity via $d_j = {\mathsf f}_{ji}s_i + n_j$, where ${\mathsf f}_{ji}$ is the so-called blurred lensing operator (mapping matrix) that incorporates the lens potential (which governs the deflection of light rays) and the PSF (blurring),¹⁰ and n_j is the noise in the data characterized by the covariance matrix ${\bm{\mathsf C}}_{\mathrm{D}}$. In the inference of s_i, we impose a prior on s_i, which can be thought of as "regularizing" the parameters s_i to avoid overfitting to the noise in the data. Following Suyu et al. (2006), we use quadratic forms of the regularization (specifically, zeroth-order, gradient, and curvature forms of regularization). The Bayesian inference of the source-intensity distribution (s_i) given the observed image (d_j) is a linear inversion and is a solved problem. Having obtained the source intensity, we can calculate the intensity deficit and source-intensity gradients.

We pixelize the lens potential to allow for a flexible parametrization scheme. To solve Equation (3), we cast it into a matrix equation and invert the linear system. To write Equation (3) in a matrix form, we discretize the lens potential on a rectangular grid of N_p pixels (which is less than the number of data pixels N_d so that the potential and source-intensity pixels are not underconstrained) and denote the potential perturbation by δψ_i, where i = 1, ..., N_p. The intensity deficit on the image grid is $\delta I_j=d_j - {\mathsf f}_{ji}s_i$, where j = 1, ..., N_d (using the notation from source-intensity reconstruction, d, $\bm{\mathsf f}$, and s are the data vector, the blurred lensing operator, and the source-intensity vector, respectively). Equation (3) now becomes

$$\begin{equation} {\bm \delta I} = \bm{\mathsf t \delta \psi} + {\bm n}, \end{equation} \tag{ 4 }$$

where $\bm{\mathsf t}$ is a N_d × N_p matrix, which incorporates the PSF, the source-intensity gradient, and the gradient operator that acts on δψ (see the appendix for the explicit form of $\bm{\mathsf t}$), and n is the noise in the data. The above equation is equivalent to

$$\begin{equation} {\bm d} = \bm{\mathsf fs} + \bm{{\mathsf t} \delta \psi}+{\bm n}. \end{equation} \tag{ 5 }$$

We can infer the potential corrections δψ given the data d, source intensity s, and source-intensity gradients that are encoded in $\bm{\mathsf t}$. In the inference, we impose a prior on δψ. The posterior probability distribution is

$$\begin{equation} \overbrace{P(\bm{\delta \psi}|{\bm d}, \bm{\mathsf f}, {\bm s}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi })}^{\rm {posterior}} = \frac{ \overbrace{P({\bm d}|\bm{\delta \psi}, \bm{\mathsf t}, \bm{\mathsf f}, {\bm s})}^{\rm {likelihood}} \overbrace{P(\bm{\delta \psi}|\mu, \bm{\mathsf g}_{\delta \psi })}^{\rm {prior}}}{\underbrace{P({\bm d}|\bm{\mathsf f}, {\bm s}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi })}_{\rm {evidence}}}, \end{equation} \tag{ 6 }$$

where μ and $\bm{\mathsf g}_{\delta \psi }$ are the (fixed) strength and form of regularization for the potential correction inversion, and all irrelevant (in)dependences have been dropped. Modeling the noise as Gaussian, the likelihood is

$$\begin{equation} P({\bm d}|\bm{\delta \psi}, \bm{\mathsf t}, \bm{\mathsf f}, {\bm s}) = \frac{\exp (-E_{\rm {D}}({\bm d}|\bm{\delta \psi}, \bm{\mathsf t}, \bm{\mathsf f}, {\bm s}))}{Z_{\rm {D}}}, \end{equation} \tag{ 7 }$$

where

$$\begin{eqnarray} &&E_{\rm {D}}({\bm d}|\bm{\delta \psi}, \bm{\mathsf t}, \bm{\mathsf f}, {\bm s}) = \frac{1}{2}({\bm d}-\bm{\mathsf f}s - \bm{\mathsf t}\bm{\delta \psi})^{\rm {T}} \bm{\mathsf C}_{\mathrm{D}}^{-1} ({\bm d}-\bm{\mathsf f}s - \bm{\mathsf t}\bm{\delta \psi})\nonumber\\ \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&\hspace{-58pt} = \frac{1}{2}\chi ^2, \end{eqnarray} \tag{ 9 }$$

and Z_D is the normalization for the probability. We express the prior in the following form:

$$\begin{equation} P(\bm{\delta \psi}|\mu, \bm{\mathsf g}_{\delta \psi }) = \frac{\exp (-\mu E_{\mathrm{\delta \psi }}(\bm{\delta \psi}|\bm{\mathsf g}_{\delta \psi }))}{Z_{\delta \psi }(\mu)}. \end{equation} \tag{ 10 }$$

We use quadratic forms of the regularizing function E_δψ. In particular, we use the curvature form of regularization (see for example, Appendix A of Suyu et al. 2006 for an explicit expression of the curvature form of regularization). We use this regularization instead of the zeroth-order or gradient forms because the lens potential should in general be smooth, being the integral of the surface mass density. Curvature regularization in the potential corrections effectively corresponds to zeroth-order regularization in the surface mass density corrections. This implies a prior preference toward zero surface mass density corrections, thus suppressing the addition of mass to the initial mass model unless the data require it.

Maximizing the posterior of parameters δψ, we obtain the most probable solution

$$\begin{equation} \bm{\delta \psi}_{\mathrm{MP}} = \bm{\mathsf A}^{-1} {\bm D}, \end{equation} \tag{ 11 }$$

where

$$\begin{eqnarray} \nonumber \bm{\mathsf A} & = & \bm{\mathsf B} + \mu \bm{\mathsf C}, \\ \nonumber \bm{\mathsf B} & \equiv & \nabla \nabla E_{\mathrm{D}}(\bm{\delta \psi}) = \bm{\mathsf t}^{\mathrm{T}} \bm{\mathsf C}_{\mathrm{D}}^{-1} \bm{\mathsf t}, \\ \nonumber \bm{\mathsf C} & \equiv & \nabla \nabla E_{\mathrm{\delta \psi }}(\bm{\delta \psi}), \\ \nonumber D & = & \bm{\mathsf t}^{\mathrm{T}} \bm{\mathsf C}_{\mathrm{D}}^{-1} ({\bm d}-\bm{\mathsf f}{\bm s}), \\ \nonumber \rm {and\ } \nabla & \equiv & \frac{\partial }{\partial \bm{\delta \psi}}. \end{eqnarray}$$

The matrices $\bm{\mathsf A}$, $\bm{\mathsf B}$, and $\bm{\mathsf C}$ have dimensions N_p × N_p and are, by definition, the Hessians of the exponential arguments in the posterior, the likelihood, and the prior probability distributions, respectively.

As discussed in detail in, for example, MacKay (1992) and Suyu et al. (2006), the evidence is irrelevant in the first level of inference where we maximize the posterior of parameters δψ to obtain the most probable parameters δψ_MP. However, the evidence is crucial for the second level of inference for model comparison, where a model incorporates the lens potential, PSF, and regularizations of both the source intensity and the potential correction. If we assert that models are equally probable a priori, then the evidence gives the relative probability of the model given the data. In other words, the ratio in the evidence values of two models tells us how much more probable the first model is relative to the second model if we assume that the two models are a priori equally probable. Since the evidence gives only the relative probability, the data set needs to be kept the same for model comparison.

The posterior ($P(\bm{\delta \psi}|{\bm d}, \bm{\mathsf f}, {\bm s}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi })$) and the evidence ($P({\bm d}|\bm{\mathsf f}, {\bm s}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi })$) in Equation (6) are conditional on the source-intensity distribution. Ideally, we would have an expression of the posterior for both s and δψ: $P({\bm s}, \bm{\delta \psi}|{\bm d}, \bm{\mathsf f}, \lambda, \bm{\mathsf g}_{\rm S}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi })$, where λ and $\bm{\mathsf g}_{\rm S}$ are, respectively, the strength and form of regularization for s. We would also obtain the evidence by marginalizing both the source-intensity and the potential correction values, $P({\bm d}|\bm{\mathsf f}, \lambda, \bm{\mathsf g}_{\rm S}, \bm{\mathsf t}, \mu, \bm{\mathsf g}_{\delta \psi }) = \int \rm {d}{\bm s}\, \rm {d}\bm{\delta \psi}\,P({\bm d}|{\bm s}, \delta \psi, \bm{\mathsf f}, \bm{\mathsf t})\,P\,({\bm s}, \bm{\delta \psi}| \lambda, \bm{\mathsf g}_{\rm S}, \mu, \bm{\mathsf g}_{\delta \psi })$. However, due to the iterative nature of the method (i.e., s and δψ are not inferred simultaneously), we do not have such expressions for the posterior and the evidence. Pragmatically, we use the evidence from the source reconstruction (given the corrections δψ), $P({\bm d}| \bm{\mathsf f}, \bm{\delta \psi}, \lambda, \bm{\mathsf g}_{\rm S})$, for comparing the potential models, PSF, and regularizations. Specifically, after iterating through the source-intensity reconstructions and lens potential corrections, we use the final corrected lens potential for one last source-intensity reconstruction and use the evidence from this final source reconstruction for comparing models. This approximation is valid provided that the probability distributions of δψ and the regularization constant are sharply peaked at the most probable values. Suyu et al. (2006) showed that the delta function approximation for the regularization constant is acceptable; simulations of the iterative potential reconstruction method suggest that the probability of δψ after the final iteration is sharply peaked. Therefore, the probability of a given potential model, PSF, and form of regularization is $P(\bm{\mathsf f}, \bm{\mathsf g}_{\rm S}|{\bm d}) \propto \int \rm {d}\lambda \, \rm {d}\bm{\delta \psi}\, P({\bm d}| \bm{\mathsf f},\bm{\delta \psi}, \lambda, \bm{\mathsf g}_{\rm S})\,P\,(\bm{\mathsf f}, \bm{\mathsf g}_{\rm S}) \sim P({\bm d}| \bm{\mathsf f},\hat{\bm{\delta \psi}}, \hat{\lambda }, \bm{\mathsf g}_{\rm S}) P(\bm{\mathsf f}, \bm{\mathsf g}_{\rm S})$, where $\hat{\bm{\delta \psi}}$ and $\hat{\lambda }$ are the most probable solutions. Assuming that all models are equally probable a priori (i.e., $P(\bm{\mathsf f}, \bm{\mathsf g}_{\rm S})$ is constant), the evidence from the source reconstruction serves as a reasonable proxy to use for model comparison.

There is an uncertainty associated with the evidence values due to finite source-intensity resolution as a result of the source pixelization. The source reconstruction region is initially chosen such that the mapped source region on the image plane encloses the Einstein ring. This ensures that the source region contains the entire source-intensity distribution. Throughout the iterative pixelated potential reconstruction, the source region and pixelization are kept the same. In the final source reconstruction for evidence computation, the evidence value depends on the pixelated source region because the goodness of fit on the image plane generally changes, especially in areas of significant intensity gradients, as one shifts the source region. To estimate the uncertainty in the evidence values, we perform the last source reconstruction for various source regions that are shifted by a fraction of a pixel from the optimized one in the potential reconstruction. The range of the resulting evidence values for the various source regions then allow us to quantify the uncertainty in the evidence. In addition to the uncertainty due to source pixelization, the evidence also depends on the amount of regularization on δψ, which is discussed in Section 2.2.2.

Solving for the potential perturbations is very similar to solving for the source-intensity distribution in Suyu et al. (2006) except for the following technical details.

• 1.  
  In each iteration, the perturbative potential correction is obtained only in an annular region instead of over the entire lens potential grid due to the need for the source-intensity gradient (see Equation (3)) to be measurable. Since the extended source intensity is only non-negligible near the Einstein ring, we only have information about the source-intensity gradients in this region. In practice, the annular region is the mapping of the finite source reconstruction grid that encloses the extended source with a minimal number of source pixels (for computational efficiency). The annulus of potential corrections obtained at each iteration is extrapolated for the next iteration by minimizing the curvature in the potential corrections. This allows the shape of the annular region to change as needed when the lens potential gets corrected. In addition, the forms of the regularization matrix, as discussed in Appendix A of Suyu et al. (2006), are modified accordingly to take into account the nonrectangular reconstruction region (described in more detail in the third point below).
• 2.  
  Since Equation (5) is a perturbative equation in δψ, the inversion needs to be over-regularized to enforce a small correction in each iteration. Empirically, we set the regularization constant, μ, at roughly the peak of the function μE_δψ (within a factor of 10), which corresponds to the value before which the prior dominates. The resulting evidence value from the final source-intensity reconstruction weakly depends on the value of μ, and we include this dependence in the uncertainty of the evidence value.
• 3.  
  The potential corrections are generally nonzero at the edge of the annular reconstruction region. This calls for slightly different structures of regularization compared to those written in Appendix A in Suyu et al. (2006) for source-intensity reconstruction (since the source grids are chosen to enclose the entire extended source such that edge pixels have nearly zero intensities). The regularizations are still based on derivatives of δψ; however, no patching with lower derivatives should be used for the edge pixels because the zeroth-order regularization at the top/right edge will incorrectly enforce the δψ values to zero in those areas. The absence of the lower derivative patches leads to a singular regularization matrix,¹¹ which is problematic for evaluating the Bayesian evidence for lens potential correction. However, since we do not use the evidence values to compare the forms of regularization for the potential corrections (because we use only the curvature form) nor to compare the lens potential and PSF model, the revised structure of regularization is acceptable. We have found this structure of regularization for potential corrections to work for various types of sources (with varying sizes, shapes, number of components, etc.).

In the source reconstruction steps of this iterative scheme, we discover by using simulated data that over-regularizing the source reconstruction in early iterations helps the process to converge. This is because initial guess potentials that are significantly perturbed from the true potential often lead to highly discontinuous source distributions when optimally regularized (corresponding to maximal Bayesian evidence, which balances the goodness of fit and the prior), and over-regularization would give a more regularized source-intensity gradient for the potential correction. Unfortunately, we do not have an objective way of setting this over-regularization factor for the source reconstruction. Currently, at each source reconstruction iteration, we set the over-regularization factor such that the magnitude of the intensity deficit is at approximately the same level as that from the optimally-regularized case but with a smoother source-intensity distribution for numerical derivatives. This scheme ensures that we do not over-regularize when we are close to the true potential. Based on simulated test runs, the recovery of the true potential depends on the amount of over-regularization. When the initial guess is far from the true potential, over-regularization in the early iterations is crucial for convergence. We find that it is better to over-regularize in excess than in deficit. Too much over-regularization simply leads to more iterations to converge, whereas too little over-regularization may not converge at all.

For each iteration of source-intensity reconstruction, there is also a mask on the source plane to exclude source pixels that either (1) are not mapped by that iteration's lens potential on the data grid or (2) have no neighboring pixels for the computation of numerical derivatives. We generalize the regularizing function for this nonrectangular reconstruction region to have the right-most and top-most pixels (pixels adjacent to the edge or adjacent to the masked source pixels) patched with lower derivatives as we did for the edge pixels in Appendix A of Suyu et al. (2006). This patching ensures that the regularization matrix is nonsingular for the evaluation of the Bayesian evidence.

Based on simulated test runs, we find that a practical stopping criterion for the iterative procedure is to terminate when the relative potential corrections between all image pairs are (δψ₁ − δψ₂)/(ψ₁ − ψ₂) < 0.1%, where 1 and 2 label the images in any pair. After this criterion is reached, further iterations give a negligible contribution to the predicted Fermat potential differences between the images.

The restriction to using only isophotal data implies that the potential correction we obtain at each iteration may be affected by the "mass-sheet degeneracy" (Falco et al. 1985). However, the addition of mass sheets is suppressed by the curvature form of the regularization for the potential correction and also by the large amount of over-regularization. We refer to Kochanek et al. (2006) and Paper II for a detailed description of the mass-sheet degeneracy; here we review a few key points that are relevant for the potential corrections. In essence, an arbitrary symmetric paraboloid, gradient sheet, and constant can be added to the potential without changing the predicted lensed image:

$$\begin{equation} \psi _{\nu }(\vec{\theta }) = \frac{1-\nu }{2} \vert \vec{\theta } \vert ^2 + \vec{a}\cdot \vec{\theta } + c + \nu \psi (\vec{\theta }), \end{equation} \tag{ 12 }$$

where $\psi (\vec{\theta })$ is the original potential, $\psi _{\nu }(\vec{\theta })$ is the transformed potential, and ν, $\vec{a}$, and c are constants. The constants $\vec{a}$ and c have no physical effects on the lens systems as they merely change the origin on the source plane (which is unknowable) and change the zero point of the potential (which is not observable). The parameter 1 − ν refers to the amount of mass sheet, which can be seen in the corresponding convergence transformation: $\kappa _{\nu }(\vec{\theta })=(1-\nu)+\nu \kappa (\vec{\theta })$. To make sure we remain "close" to the initial potential model, we set ν = 1 and fix three points in the corrected potential after each iteration to the corresponding values of the initial potential. Setting ν = 1 ensures that the size of the extended source intensity remains approximately the same, and the three fixed points allow us to solve for $\vec{a}$ and c in Equation (12) to remove irrelevant gradient sheets and constants in the reconstructed potential. We choose the three points to be three of the four (top, left, right, and bottom) locations of the annular reconstruction region that are midway in thickness between the annular edges. The three points are usually chosen to be at places with lower surface brightness in the ring. This technique of "fixing" the mass-sheet degeneracy is demonstrated in Section 2.2.4 using simulated data.

To summarize, the steps for the iterative and perturbative potential reconstruction scheme via matrices are as follows. (1) Reconstruct the source-intensity distribution given the initial (or corrected) lens potential based on Suyu et al. (2006). (2) Compute the intensity deficit and the source-intensity gradient. (3) Solve Equation (5) for the potential corrections δψ in the annulus of reconstruction. (4) Update the current potential using Equation (2): ${\bm \psi}_{\rm {next \ iteration}} = {\bm \psi}_{\rm {current \ iteration}}-\bm{\delta \psi}$. (5) Transform the corrected potential ${\bm \psi}_{\rm {next \ iteration}}$ via Equation (12) so that ν = 1 and the transformed corrected potential has the same values as the initial potential at the three fixed points. (6) Extrapolate the transformed corrected potential for the next iteration. (7) Interpolate the transformed corrected potential onto the resolution of the data grid for the next iteration's source reconstruction. (8) Repeat the process using the extrapolated and finely gridded reconstructed potential, and stop the process when the relative potential correction between any pair of images is <0.1%.

The ACS data were reduced using the multidrizzle package (Koekemoer et al. 2002) in an early version of the HAGGLeS image-processing pipeline (P. J. Marshall et al. 2009, in preparation), producing drizzled images with a 003 pixel scale. The drizzled ACS images are shown in Figure 3. The corresponding output weight images from multidrizzle give the values for the inverse variance of each pixel. We approximate the noise covariance matrix as diagonal and use the variance pixel values for the diagonal entries, even though drizzling will correlate the noise between adjacent pixels. It is assumed that the effect of drizzling can be modeled as having a diagonal covariance matrix with the diagonal elements rescaled (Casertano et al. 2000). In practice, we do not need to do the rescaling because the ranking of the models using the relative log evidence values from the source reconstruction is insensitive to rescaling of the covariance matrix.

A pixelated representation of a continuous intensity distribution generally introduces error in the interpolated intensity values between pixels, especially for intensity distributions with sharp features. This error should be incorporated into the likelihood function. Therefore, for modeling the source-intensity distribution on a grid (in Sections 5.6 and 6), we also include the error due to pixelization on the image and source planes (which we call "regridding error") in the image covariance matrix. We express the regridding error on the image plane in terms of the data (instead of on the source plane and transforming it to the image plane) in order to obtain a noise map that is independent of the pixelated lens modeling. The regridding error associated with pixel i is

$$\begin{equation} \big(\sigma _{\rm grid}^2\big)_i = \frac{1}{12} \mu _i \frac{\Delta \beta ^2}{\Delta \theta ^2} \sum _{{ j \in\ {\rm pixels}\ {\rm adjacent\ to\ }\ i}}^{N_{\rm adj}} \frac{(d_j-d_i)^2}{N_{\rm adj}}, \end{equation} \tag{ 15 }$$

where μ_i is the lensing magnification at pixel i, Δβ is the source pixel size, Δθ is the image pixel size, N_adj is the number of pixels adjacent to pixel i, and d_i (d_j) is the image intensity at pixel i (j). The summation divided by Δθ² in the above equation is a conservative estimate on the error due to pixelization on the image plane. Since sharper features in the image have larger gradients (hence, larger values for the summations), the regridding error is higher in these areas by construction. The remaining quantities in the equation, μ_iΔβ²/12, account for the uncertainty in the predicted image (the source image mapped to the image plane) due to the pixelization of the source-intensity distribution. The factor 1/12 is the second moment of a uniform distribution between −0.5 and 0.5. When one constructs the predicted image by mapping each image pixel to the source plane and reading off the source-intensity value, the mapped source position (of an image pixel) is generally not centered on a source pixel, but have on average a $(1/\sqrt{12})$-pixel shift from the center of the source pixel. Therefore, $\Delta \beta /\sqrt{12}$ is the effective size of the source pixel, which is then magnified by (on average) $\sqrt{\mu _i}$ due to lensing. In the pixelated potential reconstruction, we approximate the magnification at each image pixel (which requires the second derivative of the potential) by the value computed from the initial potential because (1) the approximation enforces the regridding error to be independent of the pixelated potential modeling and (2) the corrected potential values are obtained on an annular region only a few pixels thick. Having obtained an estimate for the regridding error, we add it in quadrature to the variance from the weight image to obtain the entries of the approximated diagonal covariance matrix.

The inclusion of the regridding error is important for source-intensity reconstructions with sharp intensity features (such as the presence of a bright core); it has the effect of stabilizing the evidence values with respect to choices in the source pixelization. Without including the regridding error, a pixelated description of, for example, a source-intensity distribution with a bright core would be highly sensitive to the centering of the core on the source pixels. A small mismatch could create large image residuals near the cores that would veto an otherwise good lensing model, which has the rest of the extended features well described. Such an undesirable effect is mostly removed by the inclusion of the regridding error. For B1608+656, the ratio of the regridding error to the error from the multidrizzle weight image is around ∼30 near the image centroids and ∼1 in other parts in the Einstein ring.

The NICMOS F160W image was taken from Koopmans et al. (2003). Drizzled images on rectangular grids for different instruments are generally not on the same resolution and not aligned. This is the case for the NICMOS and ACS images. We use SWarp ¹³ to align the combined NICMOS image to the ACS images. The final SWarped NICMOS F160W image with 003 pixel scales is shown in Figure 4.

Zoom In Zoom Out Reset image size

The ACS PSF is both spatially and temporally varying (e.g., Rhodes et al. 2007). One source of temporal variation is the "breathing" of the telescope while it orbits, which causes the focal length (and, hence, the PSF) of the telescope to change. Instead of adopting a universal PSF, we take the approach of modeling several PSFs using different means, and quantitatively comparing them using the Bayesian analysis described in Section 2.2.1. This has the advantage of using the data (the observed image) to rank the models. For each of the two drizzled ACS images, we create five models for the PSF either based on the TinyTim package (Krist & Hook 1997) or from the unsaturated stars in the field: (1) drizzled PSF ("PSF-drz") from a set of TinyTim simulations (following Rhodes et al. 2007), (2) single (nondrizzled) TinyTim PSF ("PSF-f3") with a telescope focus value of −3, (3) the closest star ("PSF-C") located at ∼9'' in the northeast direction from B1608+656 in the drizzled ACS field with a Vega magnitude of 21.3 in F814W, (4) bright star #1 ("PSF-B1") that is located at ∼19 southwest of B1608+656 in the drizzled ACS field with a Vega magnitude of 18.7 in F814W, and (5) bright star #2 ("PSF-B2") that is located at ∼16 south of B1608+656 in the drizzled ACS field with a Vega magnitude of 19.1.

The TinyTim frame(s) were drizzled and resampled to pixel sizes of 003 to match the resolution of the ACS images. We keep in mind that the TinyTim PSFs (PSF-drz and PSF-f3) may be insufficient due to the time-varying nature of the PSF and the aging of the detector since the TinyTim code was written. We expect the closest star to B1608+656 (PSF-C) to be a good approximation to the PSF because the spatial variation of the PSF across ∼9'' should be negligible and any temporal variations are the same as in the lens system. However, this closest star is not bright enough to see the secondary maxima in the PSF, so we additionally include two of the brightest stars in the drizzled field mentioned above. For each of the stars in F606W and F814W, we make a small cutout around the star (25 × 25 pixels for PSF-C, 51 × 51 pixels for PSF-B1, and 41 × 41 pixels for PSF-B2) and center it on a 200 × 200 grid, which is the size of the drizzled science image cutouts of B1608+656 that are used for the image processing.

The NICMOS PSF is thought to be more stable, and thus we assume a TinyTim model for it. The output TinyTim PSF is in the CCD frame of NICMOS with pixel size 0043. As with the F160W science image, the PSF was SWarped to be aligned with the ACS images with 003 pixels. Since there is only one PSF model for NICMOS, PSF specifications throughout the rest of this paper refer to the ACS PSFs.

The dust correction method outlined above requires the intrinsic colors to be determined from the color maps. To construct the color maps, we need to unify the different resolutions of the images in different bands (due to the wavelength dependence of the PSF). We do so by deconvolving the F606W, F814W, and F160W images using their corresponding PSFs, and reconvolving the images with the F814W PSF for each set of the five ACS PSFs and the single NICMOS PSF described in Section 5.3. Reconvolved images are preferred to deconvolved images, because the latter show small-scale features (of a few pixels' size) that are artificial due to the amplification of the noise during the deconvolution process. We select the F814W PSF for the reconvolution because F814W will be used for the lens potential modeling, due to its high S/N compared with F160W and its less severe dust extinction compared with F606W. In working with the reconvolved images, we assume that the dust varies on a scale larger than the F814W PSF, which is true for the regions near the Einstein ring. For the deconvolution, we use IDL's max_entropy iterative routine that is based on the algorithm by Hollis et al. (1992). We were unable to deconvolve the ACS F814W image using PSF-f3. This suggests that PSF-f3 is a bad model, which we have expected due to temporal variations in the PSF. PSF-f3 is a single-epoch PSF whereas the F814W image was drizzled from multiple exposures. We, therefore, discard this PSF model.

For each set of PSF models (PSF-drz, PSF-C, PSF-B1, and PSF-B2 for ACS, and TinyTim PSF for NICMOS), we construct the color maps F606W−F814W, F606W−F160W, and F814W−F160W from the reconvolved F606W, F814W, and F160W images, respectively. Figure 5 shows the three color maps derived for PSF-B1. Regions with bluer color slightly west of G1 are shown in all three color maps. Since the centroid of this blue region is offset from the centroid of G1, we believe that this blue region arises from differential reddening and not from intrinsic color variations within G1, which is an elliptical galaxy (Surpi & Blandford 2003). Since elliptical galaxies typically contain little dust, Koopmans & Fassnacht (1999) and Surpi & Blandford (2003) suggested that the dust comes from G2, likely a dusty late-type galaxy, through dynamical interaction. This may explain why the spectrum of G1 shows signatures of a young stellar population plus a poststarburst population (Dressler & Gunn 1983; Myers et al. 1995; Surpi & Blandford 2003; Koopmans et al. 2003): gas from G2 may have been transferred to G1, where the tidal interactions may have triggered star formation.

Zoom In Zoom Out Reset image size

The color maps also show regions of bluer color around images C and D, and we again believe that these are mostly differential reddening due to the misalignment of the image positions and the centroids of these blue regions, especially in F606W−F160W and F814W−F160W. Furthermore, we find more dust at the crossing point of the isophotal separatrix (the figure-eight-shaped intensity contour) of the image pair A–C. This is encouraging, as lensing models indeed predict the crossing point to be closer to image A (see discussion in Section 5.4.2). However, these bluer regions near images C and D may also arise from the lensed source being intrinsically bluer than the surrounding emission. The F814W−F606W color for these blue regions is consistent with typical star-forming galaxies (e.g., Coleman et al. 1980). In the F606W−F814W color map, there is a faint ridge of redder color connecting images A and C. This may be due to the asymmetry in the stellar PSF model (with the star position not exactly centered within a pixel), which would cause the F606W and F814W isophotes to shift relative to each other after the deconvolution and reconvolution. For the color maps from the other PSF models, we find that the color maps from PSF-C and PSF-B2 look similar to that from PSF-B1 with varying amounts of noise due to varying brightnesses of the stellar PSFs. PSF-drz gave color maps that differ from those from the stellar PSFs (PSF-C, PSF-B1, and PSF-B2) because PSF-drz, especially in the F606W band, did not exhibit a single brightness peak but a string of equal brightness pixels at the center due to frame alignment difficulties during the drizzling process. This caused the brightest pixels in the Einstein ring to shift by ∼1 pixel after the deconvolution and the reconvolution process in F606W, and created artificial sharp highlights tracing the edge of the ring in the F606W−F814W color map. As will be seen in Section 5.6, this leads to PSF-drz and its resulting dust map giving a lower goodness of fit in the lens inversion, and hence being ranked lower compared with other models.

In each of the color maps, we define three color regions for the three color components: one within the Einstein ring for the lens galaxies (we assume G1 and G2 to have the same colors), one for the Einstein ring of the lensed extended source, and one for the lensed AGN (core of the extended source). Following Koopmans et al. (2003), we determine the bluest color within each region, assume that this part of the region was not absorbed by dust, and adopt this color as the intrinsic color. This assumes that each of the three components has a constant intrinsic color. This would allow us to obtain the differential reddening for each of the components across the lensed image; absolute reddening is not needed because a uniform dust screen does not affect lens modeling. Table 3 lists the intrinsic colors for each of the three pairs of color maps. The intrinsic colors of F606W−F814W are not identical to the difference between F606W−F160W and F814W−F160W, but agree within the uncertainties (0.02–0.1).

With the intrinsic colors determined for each PSF model, we obtain two dust maps (A_V maps) using (1) only the ACS F606W and F814W images and (2) the ACS F606W and F814W images together with the NICMOS F160W image. In this way, we can assess whether the inclusion of the lower S/N NICMOS image (with the much broader PSF) improves the dust correction.

The left-hand panel of Figure 6 is the resulting A_V dust map derived using PSF-B1 and using images in all three bands. The dust map shows the east-west dust lane through the system (absorbing light from C, G2, G1, and D) that is visible in the original drizzled ACS F606W and F814W images. There is little extinction near images A and B, but there are faint rings surrounding the images that are mostly due to imperfect F160W deconvolution. We note that the low S/N exterior to the Einstein ring results in the dust map being noisy in this area. We make sure that these noisy areas are not included in the Bayesian evidence computations in Sections 5.6 and 6. The right-hand panel of Figure 6 is the resulting dust-corrected F814W image that exhibits two signs of proper dust correction: the correctly shifted crossing point of the isophotal separatrix of the image pair A–C, as shown more clearly in Figure 7, and the smoother lens galaxy profiles. As a result of recovering the absorbed light, the dust-corrected image has higher intensity values than the uncorrected image. Therefore, we create a weight map for the dust-corrected image by scaling the multidrizzle weight image in order to keep the S/N of each pixel the same (before and after dust correction). This "dust-corrected weight image" will be used in the next section for determining the lens galaxy light.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dust maps obtained from the other PSF models with or without the inclusion of the NICMOS image show similar features except for the following two dust maps.

• 1.  
  The ACS-only (no NICMOS) dust map from PSF-B2 showed a faint ridge of dust connecting images A and C. As explained, this may be due to the asymmetrical/bad PSF model. Since the dust map otherwise exhibits the correct features, we keep this dust map for the next analysis step.
• 2.  
  The ACS-only dust map from PSF-drz showed prominent artificial lensing arc features due to the ∼1 pixel offset in the image positions/arcs in the deconvolved and reconvolved F606W and F814W images, respectively. Therefore, we discard this dust map of the ACS-only images for PSF-drz, but keep the dust map derived from using all three bands (that includes NICMOS).

After discarding the ACS PSF-f3 and the ACS-only dust map from PSF-drz, we have a total of seven dust maps (and resulting dust-corrected F814W images). All of these are reasonable dust corrections to use since they are derived using representative PSFs and intrinsic colors. We will compare these dust maps and PSF models in Section 5.6.

Table 5 summarizes the results of model comparison. The "Reg. Type" column denotes the preferred type of regularization for the source reconstruction based on the highest Bayesian evidence value (Suyu et al. 2006). It can be one of the three types that we use: zeroth-order, gradient, and curvature. The last column lists the log evidence values from the inversions. Assuming the different models to be equally probable a priori, we use these evidence values for model comparison. The log evidence values range from 1.1 × 10⁴ to 1.6 × 10⁴ with uncertainties of ∼0.03 × 10⁴ due to the finite source resolution.

The list shows that the three-band dust models have higher evidence values than the two-band dust models. This is attributed to the two-band dust models showing image residuals from the aforementioned artificial surpluses and deficits of light in the dust-corrected image. The inclusion of the NICMOS F160W image to the ACS images (F606W and F814W) for the dust correction is, therefore, crucial due to (1) the proximity in the wavelengths of the ACS images and (2) the reduction in the error associated with image misalignments and simplistic intrinsic color models.

The three-band dust models also have higher evidence values than the no-dust model. This further validates the three-band dust correction, as already indicated by Figure 7. The evidence value of the no-dust model is in midst of the values for the two-band dust models, suggesting that the systematic effects in the two-band dust maps are comparable to the corrections that the dust maps are meant to achieve, thus leading to little improvement in the lens modeling.

The difference between the evidence values in Models 1 and 2 (where the models only differ in the Sérsic light profiles) is, in general, smaller than the difference between one of these two models and another PSF/dust model. Therefore, the source reconstruction (part of lens modeling) seems to be less sensitive to the galaxy light profiles than the PSF/dust models. This is in agreement with our finding that the dust-corrected image depends more on the PSF and the intrinsic color models than on the form of the lens galaxy light and the dust associated with the lens galaxies. Models 1 and 2 with PSF-drz have log evidence values on the low side of the collection of models with three-band dust maps, which was expected with PSF-drz not having a single brightness central peak due to misalignments in the drizzling process. The other models with three-band dust maps (Models 3, 5, 7, and 9) have effectively the same evidence values within the uncertainties. The models with two-band dust maps (Models 4, 6, 8, and 10) lead to a range of evidence values with the PSF-C dust map being preferred to the PSF-B1 and PSF-B2 dust maps. The two-band dust maps suggest that the shape of the primary maximum in the PSF is more important in the modeling than the inclusion of secondary maxima since PSF-C, which we expect to have a more accurate shape for the primary PSF maximum than PSF-B1 and PSF-B2, does not have the secondary maxima whereas PSF-B1 and PSF-B2 do. The asymmetry in the PSF due to the star not being centered on a single pixel may also explain the less-preferred PSF-B1 and PSF-B2. The distinction between the various stellar PSFs vanishes with the three-band dust maps, possibly due to the higher amount of noise in the three-band dust map with the inclusion of the lower S/N NICMOS image. In this case, the effects of the PSF variations across the field are suppressed.

All models preferred either the zeroth-order or gradient form of regularization, but never the curvature form; however, we mention that the difference in the log evidence values between the different regularization schemes (≲3 × 10²) are on the order of the uncertainties due to source pixelization, and the resulting reconstructions for different types of regularizations are almost identical. This is because differences in evidence values between models are currently dominated by changes in goodness of fit rather than subtle differences between the prior forms. Only when the image residual is reduced will the prior (regularization) begin to play a greater role in avoiding the reconstruction to fit to noise in the data by keeping the source model simple.

This section has illustrated a method of creating sensible PSF, dust, and lens galaxy light models for the gravitational lens B1608+656. We have obtained a representative sample of models, and have compared these models quantitatively. This collection of PSF, dust, and lens galaxy light models leads to image residuals that cannot be beaten down further unless we improve the SPLE1+D (isotropic) simply-parameterized lens potential model by Koopmans et al. (2003) to take into account the two interacting galaxy lenses. The pixelated potential reconstruction of B1608+656 is the subject of Section 6.