We present a strong lens analysis of SDSS J1004$$+$$4112, a unique quasar lens produced by a massive cluster of galaxies at $$z =$$ 0.68, using newly developed software for gravitational lensing. We find that our parametric mass model well reproduces all observations, including the positions of quasar images as well as those of multiply imaged galaxies with measured spectroscopic redshifts, time delays between quasar images, and the positions of faint central images. The predicted large total magnification of $$\mu \sim$$ 70 suggests that the lens system is indeed a useful site for studying the fine structure of a distant quasar and its host galaxy. The dark halo component is found to be unimodal, centered on the brightest cluster galaxy and the Chandra X-ray surface brightness profile. In addition, the orientation of the halo component is quite consistent with those of the brightest cluster galaxy and member galaxy distribution, implying that the lensing cluster is a relaxed system. The radial profile of the best-fit mass model is in good agreement with a mass profile inferred from the X-ray observation. While the inner radial slope of the dark halo component is consistent with being $$-$$1, a clear dependence of the predicted A–D time delay on the slope indicates that an additional time-delay measurement will improve constraints on the mass model.

## Introduction

SDSS J1004$$+$$4112 is a unique quasar lens system (Inada et al. 2003; Oguri et al. 2004; Inada et al. 2005, 2008). A quasar at $$z =$$ 1.734 is multiply imaged into five images, with their maximum image separation of 14$$^{\prime\prime}\!\!\!.$$7, produced by a massive cluster of galaxies at $$z =$$ 0.68. It is one of only two known examples of cluster-scale quasar lenses, the other being the triple lens SDSS J1029+2623 with the maximum image separation of 22$$^{\prime\prime}\!\!\!.$$5 (Inada et al. 2006; Oguri et al. 2008). In addition to the quasar images, SDSS J1004$$+$$4112 contains spectroscopically confirmed multiply imaged galaxies at $$z\sim 3$$ (Sharon et al. 2005). Both the quasar images and the lensing cluster have been detected in Chandra X-ray observations (Ota et al. 2006). Moreover, the time delays between some of the quasar images have also been measured from long-term optical monitoring observations (Fohlmeister et al. 2007, 2008).

Such a wealth of observational data available enable detailed investigations of the central mass distribution of the lensing cluster. Indeed, there have been several attempts to model the mass distribution of SDSS J1004$$+$$4112, using either a parametric or non-parametric method. Oguri et al. (2004) adopted a two-component (halo plus central galaxy) model to successfully reproduce the positions of four quasar images, but even the parities and temporal ordering of the quasar images could not be determined because of model degeneracies. Kawano and Oguri (2006) extended mass modeling along this line, and explored how time-delay measurements can distinguish different mass profiles. Fohlmeister et al. (2007) pointed out that it is important to include perturbations from member galaxies to reproduce the observed time delay between quasar images A and B. On the other hand, Williams and Saha (2004) and Saha et al. (2007) performed non-parametric mass modeling to show possible substructures in the lensing cluster. From non-parametric mass modeling, Saha et al. (2006) and Liesenborgs et al. (2009) concluded that the radial mass profile is consistent with that predicted in $$N$$-body simulations (Navarro et al. 1997). We summarize previous mass modeling in table 1.

Reference | Model* | Constraints (quasar)† | Constraint (galaxy)† |
---|---|---|---|

Inada et al. (2003) | SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Oguri et al. (2004) | NFW$$+$$SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Williams and Saha (2004) | non-parametric | pos (4) | $$\cdots$$ |

Sharon et al. (2005) | SPL$$+$$gals | pos (5) | pos (5) |

Kawano and Oguri (2006) | gNFW$$+$$SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Saha et al. (2006) | non-parametric | pos (5) | pos (8) |

Fohlmeister et al. (2007) | NFW$$+$$deV$$+$$gals$$+$$pert | pos$$+$$flux (5) | $$\cdots$$ |

Saha et al. (2007) | non-parametric | pos (5), $$\Delta t$$ (2) | pos (8) |

Inada et al. (2008) | NFW$$+$$SIE$$+$$gals$$+$$pert | pos$$+$$flux (5), $$\Delta t$$ (2) | $$\cdots$$ |

Liesenborgs et al. (2009) | non-parametric | pos (5), $$\Delta t$$ (3) | pos (7) |

This work | gNFW$$+$$Jaffe$$+$$gals$$+$$pert | pos$$+$$flux (5), $$\Delta t$$ (3) | pos (27) |

Reference | Model* | Constraints (quasar)† | Constraint (galaxy)† |
---|---|---|---|

Inada et al. (2003) | SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Oguri et al. (2004) | NFW$$+$$SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Williams and Saha (2004) | non-parametric | pos (4) | $$\cdots$$ |

Sharon et al. (2005) | SPL$$+$$gals | pos (5) | pos (5) |

Kawano and Oguri (2006) | gNFW$$+$$SIE$$+$$pert | pos$$+$$flux (4) | $$\cdots$$ |

Saha et al. (2006) | non-parametric | pos (5) | pos (8) |

Fohlmeister et al. (2007) | NFW$$+$$deV$$+$$gals$$+$$pert | pos$$+$$flux (5) | $$\cdots$$ |

Saha et al. (2007) | non-parametric | pos (5), $$\Delta t$$ (2) | pos (8) |

Inada et al. (2008) | NFW$$+$$SIE$$+$$gals$$+$$pert | pos$$+$$flux (5), $$\Delta t$$ (2) | $$\cdots$$ |

Liesenborgs et al. (2009) | non-parametric | pos (5), $$\Delta t$$ (3) | pos (7) |

This work | gNFW$$+$$Jaffe$$+$$gals$$+$$pert | pos$$+$$flux (5), $$\Delta t$$ (3) | pos (27) |

Name of models: SIE $$=$$ singular isothermal ellipsoid, pert $$=$$ external perturbation (e.g., external shear), NFW $$=$$ Navarro, Frenk, & White (NFW) profile, SPL $$=$$ softened power-law model, gals $$=$$ perturbations from member galaxies, typically modeled by truncated isothermal profiles, gNFW $$=$$ generalized NFW profile, dev $$=$$ de Vaucouleurs profile, Jaffe $$=$$ pseudo-Jaffe profile.

“pos” indicates constraints from image positions, “flux” from fluxes of images, and $$\Delta t$$ from time delays. Numbers in parentheses show the number of images used as constraints.

In this paper, we revisit strong lens modeling of SDSS J1004$$+$$4112 while adopting a parametric mass model. We include many observational constrains currently available, including central images of lensed quasars and galaxies, flux rations, and the time delay between quasar images (see table 1). In particular, this paper represents first parametric mass modeling that includes *both* the quasar time delays and the positions of multiply imaged galaxies as constraints. We compare our best-fit mass model with the Chandra X-ray observation of this system (Ota et al. 2006).

This paper is organized as follows. We describe our mass model in section 2. We show our results in section 3, and give our conclusion in section 2. New software for gravitational lensing, which is used for mass modeling, is presented in Appendix. Throughout the paper we adopt matter density of $$\Omega_{\rm M} =$$ 0.26 and the cosmological constant $$\Omega_\Lambda =$$ 0.74, but regard the dimensionless Hubble constant, $$h$$, as a parameter. With this choice of cosmological parameters, a physical transverse distance of 1$$h^{-1}\ $$kpc at the redshift of the lensing cluster ($$z =$$ 0.68) corresponds to 0$$^{\prime\prime}\!\!\!.$$20. We denote the angular diameter distance from observer to the lens as $$D_{\rm l}$$, from the observer to the source as $$D_{\rm s}$$, and from the lens to the source as $$D_{\rm ls}$$.

## Mass Modeling

### A Parametric Model

We modeled the main halo of the lensing cluster as the generalized NFW profile (e.g., Jing & Suto 2000). Its radial density profile is given by

In this model, the inner slope is parametrized by $$\alpha$$ (0 $$\lt \alpha \lt$$ 2); the original NFW profile corresponds to $$\alpha =$$ 1. We adopted the following modified concentration parameter as a model parameter:

We then introduced an ellipticity in the projected mass distribution by replacing $$r$$ in $$\kappa(r)$$ by the following quantity:

In this paper we took the $$x$$-axis to West and the $$y$$-axis to North, and therefore the position angle, $$\theta_{\rm e}$$, is measured East of North.

We included the brightest cluster galaxy (BCG) G1 modeled by pseudo-Jaffe Ellipsoid (Keeton 2001a). The convergence takes the following form:

Perturbations from member galaxies were also included. We modeled individual galaxies by the pseudo-Jaffe ellipsoid without the core radius [equation (9) with $$s =$$ 0]. We fixed the positions, relative luminosities, ellipticities and position angles of member galaxies to the observed values, and assumed that the velocity dispersions and truncation radii scale with the luminosity:

The mass-to-light ratio became constant with this scaling. We included 14 member galaxies within $$\lesssim\ $$20$$^{\prime\prime}$$ from the center, which were selected from the *gri*-band Subaru Suprime-cam images (Oguri et al. 2004). We adopted $$r$$-band luminosities for the scaling.

To achieve a better fit, we also included several additional perturbations. We considered general perturbations whose lens potentials $$\phi$$ are described as (see, e.g., Evans & Witt 2003; Kawano et al. 2004; Congdon & Keeton 2005; Yoo et al. 2006)

In this paper we included four perturbation terms with $$m =$$ 2 (external shear: e.g., Keeton et al. 1997), 3, 4, and 5.

### Observational Constraints

We adapted the positions of five quasar images measured by Inada et al. (2005) using the Hubble Space Telescope Advanced Camera for Surveys (HST/ACS) F814W image. Considering possible effects of microlensing or small-scale structure, we adopted conservative positional errors of 0$$^{\prime\prime}\!\!\!.$$04, and also relative magnitudes of 0.3 (0.8) for image B–D (E). In addition, we included the measured time delays between images A and B (Fohlmeister et al. 2007) and between images A and C (Fohlmeister et al. 2008). When fitting the time delays, we allowed the Hubble constant to vary with a Gaussian prior of $$h =$$ 0.72$$\ \pm\ $$0.04.

We also included multiply imaged galaxies, identified by Sharon et al. (2005), as constraints. We revisited deep multi-band HST/ACS images (F435W, F555W, and F814W; GO-10509 and GO-10716), and identified several features associated with each lensed image. We used the positions of all these features for our mass modeling. We included central images as well (Liesenborgs et al. 2009). We assumed larger positional errors of 0$$^{\prime\prime}\!\!\!.$$4 than those of the quasar images, partly because the determination of the centroids of the extended galaxy images are much less accurate.

Figure 1 shows the HST/ACS image of SDSS J1004$$+$$4112, together with the positions of multiple images summarized in tables 2 and 3. A notable feature of this cluster strong lens system, which can easily be seen in the figure, is that multiple images are distributed in a very wide range of radius, ranging from the central images very near to the cluster center to lensed galaxy images at $$\sim$$30$$^{\prime\prime}$$ from the cluster center. This is apparently good for constraining the density profile of the lensing cluster.

Name | $$\Delta x$$ [$$^{\prime\prime}$$] | $$\Delta y $$ [$$^{\prime\prime}$$] | $$\Delta m $$ | $$\Delta t$$ [d] |
---|---|---|---|---|

A | 0.000 | 0.000 | $$\equiv\,$$0 | $$\equiv\,$$0 |

B | $$-$$1.317 | 3.532 | 0.35$$\,\pm\,$$0.3 | $$-$$40.6$$\,\pm\,$$1.8 |

C | 11.039 | $$-$$4.492 | 0.87$$\,\pm\,$$0.3 | $$-$$821.6$$\,\pm\,$$2.1 |

D | 8.399 | 9.707 | 1.50$$\,\pm\,$$0.3 | $$\cdots$$ |

E | 7.197 | 4.603 | 6.30$$\,\pm\,$$0.8 | $$\cdots$$ |

Name | $$\Delta x$$ [$$^{\prime\prime}$$] | $$\Delta y $$ [$$^{\prime\prime}$$] | $$\Delta m $$ | $$\Delta t$$ [d] |
---|---|---|---|---|

A | 0.000 | 0.000 | $$\equiv\,$$0 | $$\equiv\,$$0 |

B | $$-$$1.317 | 3.532 | 0.35$$\,\pm\,$$0.3 | $$-$$40.6$$\,\pm\,$$1.8 |

C | 11.039 | $$-$$4.492 | 0.87$$\,\pm\,$$0.3 | $$-$$821.6$$\,\pm\,$$2.1 |

D | 8.399 | 9.707 | 1.50$$\,\pm\,$$0.3 | $$\cdots$$ |

E | 7.197 | 4.603 | 6.30$$\,\pm\,$$0.8 | $$\cdots$$ |

The quasar redshift is $$z_{\rm s} =$$ 1.734. The positional error is assumed to 0$$^{\prime\prime}\!\!\!.$$04 for all the lensed quasar images.

Name | $$z_{\rm s}$$ | $$\Delta x$$ [$$^{\prime\prime}$$] | $$\Delta y $$ [$$^{\prime\prime}$$] |
---|---|---|---|

A1.1 | 3.33 | 3.93 | $$-$$2.78 |

A1.2 | 1.33 | 19.37 | |

A1.3 | 19.23 | 14.67 | |

A1.4 | 18.83 | 15.87 | |

A1.5 | 6.83 | 3.22 | |

A2.1 | 3.33 | 4.13 | $$-$$2.68 |

A2.2 | 1.93 | 19.87 | |

A2.3 | 19.43 | 14.02 | |

A2.4 | 18.33 | 15.72 | |

A2.5 | 6.83 | 3.12 | |

A3.1 | 3.33 | 4.33 | $$-$$1.98 |

A3.2 | 2.73 | 20.37 | |

A3.3 | 19.95 | 13.04 | |

A3.4 | 18.03 | 15.87 | |

A3.5 | 6.83 | 3.02 | |

B1.1 | 2.73 | 8.88 | $$-$$2.16 |

B1.2 | $$-$$5.45 | 15.84 | |

B1.3 | 8.33 | 2.57 | |

B2.1 | 2.73 | 8.45 | $$-$$2.26 |

B2.2 | $$-$$5.07 | 16.04 | |

B2.3 | 8.33 | 2.57 | |

C1.1 | 3.28 | 10.25 | $$-$$3.06 |

C1.2 | $$-$$7.55 | 15.39 | |

C1.3 | 8.49 | 2.72 | |

C2.1 | 3.28 | 9.95 | $$-$$3.36 |

C2.2 | $$-$$7.30 | 15.44 | |

C2.3 | 8.49 | 2.72 |

Name | $$z_{\rm s}$$ | $$\Delta x$$ [$$^{\prime\prime}$$] | $$\Delta y $$ [$$^{\prime\prime}$$] |
---|---|---|---|

A1.1 | 3.33 | 3.93 | $$-$$2.78 |

A1.2 | 1.33 | 19.37 | |

A1.3 | 19.23 | 14.67 | |

A1.4 | 18.83 | 15.87 | |

A1.5 | 6.83 | 3.22 | |

A2.1 | 3.33 | 4.13 | $$-$$2.68 |

A2.2 | 1.93 | 19.87 | |

A2.3 | 19.43 | 14.02 | |

A2.4 | 18.33 | 15.72 | |

A2.5 | 6.83 | 3.12 | |

A3.1 | 3.33 | 4.33 | $$-$$1.98 |

A3.2 | 2.73 | 20.37 | |

A3.3 | 19.95 | 13.04 | |

A3.4 | 18.03 | 15.87 | |

A3.5 | 6.83 | 3.02 | |

B1.1 | 2.73 | 8.88 | $$-$$2.16 |

B1.2 | $$-$$5.45 | 15.84 | |

B1.3 | 8.33 | 2.57 | |

B2.1 | 2.73 | 8.45 | $$-$$2.26 |

B2.2 | $$-$$5.07 | 16.04 | |

B2.3 | 8.33 | 2.57 | |

C1.1 | 3.28 | 10.25 | $$-$$3.06 |

C1.2 | $$-$$7.55 | 15.39 | |

C1.3 | 8.49 | 2.72 | |

C2.1 | 3.28 | 9.95 | $$-$$3.36 |

C2.2 | $$-$$7.30 | 15.44 | |

C2.3 | 8.49 | 2.72 |

The positional error is assumed to 0$$^{\prime\prime}\!\!\!.$$4 for all the lensed galaxy images. The redshifts are measured spectroscopically.

We also added several Gaussian priors to the mass model. Based on the measurement by Inada et al. (2008), we assumed the velocity dispersion of the central galaxy, G1, to be $$\sigma =$$ 352$$\ \pm\ $$13 km s$$^{-1}$$. ^{1} The position of G1 was fixed to the observed position, (7$$^{\prime\prime}\!\!\!.$$114, 4$$^{\prime\prime}\!\!\!.$$409) in our coordinate system, whose origin was at the quasar image A. From the observed position and shape, we assumed the ellipticity and the position angle to 0$$^\circ\!\!\!.$$3$$\ \pm\ $$0$$^\circ\!\!\!.$$05 and 152$$^\circ\ \pm\ $$5$$^\circ$$, respectively. Furthermore, we added a weak prior to the truncation radius, $$a =$$ 8$$^\prime\prime{\ \pm\ }$$4$$^{\prime\prime}$$, based on the observed correlation between the velocity dispersion and the truncation radius (Natarajan et al. 2009).

### Model Optimization

We used software named *glafic* for all calculations of the lens properties and model optimizations (see Appendix 1). We employed a standard $$\chi ^2$$ minimization to find the best-fit mass model. Specifically, we adopted a downhill simplex method to find a minimum. To speed up the calculations, we estimated $$\chi ^2$$ in the source plane, which is found to be sufficiently accurate for our purpose. See Appendix 2 for detailed discussions about the source plane $$\chi ^2$$ minimization.

## Result

### Best-Fit NFW Model

First, we fixed the inner slope of the dark halo component to $$\alpha =$$ 1 (i.e., the NFW profile) and derived the best-fit mass model. Figure 2 shows the best-fit critical curves. It can be seen that the best-fit model reproduces the observed multiple image families quite well. In addition, it successfully reproduces the observed time delays between quasar images. The best-fit model has $$\chi ^2 =$$ 31 for 39 degrees of freedom, suggesting that our choice of errors was reasonable. The contribution from each source is reasonably similar, $$\chi ^2 =$$ 4.7 from the quasar, 21 from the galaxy A, 0.9 from galaxy B, and $$2.6$$ from galaxy C. The best-fit model predicts magnifications of the five quasar images to $$\mu_{\rm A} =$$ 29.7, $$\mu_{\rm B} =$$ 19.6, $$\mu_{\rm C} =$$ 11.6, $$\mu_{\rm D} =$$ 5.8, and $$\mu_{\rm E} =$$ 0.16. The total magnification for all the quasar images is $$\mu_{\rm tot} =$$ 67. The model also predicts the time delay between quasar images A and D to be $$\Delta t_{\rm AD} = \Delta t_{\rm D} - \Delta t_{\rm A} =$$ 1218 d, and that between quasar images A and E to be $$\Delta t_{\rm AE} = \Delta t_{\rm E} - \Delta t_{\rm A} =$$ 1674 d. The AD time delay is slightly smaller than the lower limit reported by Fohlmeister et al. (2008), $$\Delta t_{\rm AD} \gt $$ 1250 d.

We find that the best-fit centroid of the dark halo (NFW) component is (6.92$$^{+0.20}_{-0.32}$$, 4.25$$^{+0.31}_{-0.24}$$) at the 95% confidence limit, which is quite consistent with the observed position of G1. The result is in marked contrast with Oguri et al. (2004), in which significant offsets between the center of the dark halo component and that of G1 have been reported based on the modeling of quasar images A–D. Such a large offset is no longer allowed because of additional constraints from multiply imaged galaxies, time delays, and central images. The best-fit center of the dark halo component is also consistent with the X-ray center reported by Ota et al. (2006). Furthermore, the best-fit position angle of $$\theta_{\rm e} =$$ 152$$^\circ\!\!\!.$$9 for the dark halo component is quite consistent with the position angle of G1, and also that of the member galaxy distributions studied in Oguri et al. (2004). The concentricity and alignment between the dark matter, BCG, and X-ray implies that the lensing cluster is a highly relaxed system (see also Liesenborgs et al. 2009). The best-fit ellipticity of the halo component is $$e =$$0.26.

The best-fit parameters for the perturbations terms are ($$\epsilon, \theta_\epsilon$$) $$=$$ (0.040, 51.8) for $$m =$$ 2, (0.019, 114) for $$m =$$ 3, (0.013, 47) for $$m =$$ 4, and (0.010, 16.5) for $$m =$$ 5. Thus the perturbations are rather small, but they are still important for an accurate reproduction of lensed images, particularly for those of the quasar (see also Oguri et al. 2004).

One of the most important quantities to characterize a strong lensing system is the Einstein radius, $$r_{\rm Ein}$$. We compute the Einstein radii, $$r_{\rm Ein}$$, for our best-fit mass model using the following relation:

We find that $$r_{\rm Ein} =$$ 8$$^{\prime\prime}\!\!\!.$$14 for the quasar redshift $$z_{\rm s} =$$ 1.734, and 13$$^{\prime\prime}\!\!\!.$$38 for the redshift of the lensed galaxy A, $$z_{\rm s} =$$ 3.33. If we compute $$r_{\rm Ein}$$ only from the dark-halo component, excluding any contributions from galaxies, we obtain $$r_{\rm Ein} =$$ 4$$^{\prime\prime}\!\!\!.$$84 for $$z_{\rm s} =$$ 1.734 and 10$$^{\prime\prime}\!\!\!.$$31 for $$z_{\rm s} =$$ 3.33, which are quite different from those computed from the total mass distribution. This suggests that the effect of the BCG G1 on the lens system is quite significant.

### Comparison with X-Ray

Next, we compare the best-fit radial mass profile derived from strong lensing with that inferred from the Chandra X-ray observation (Ota et al. 2006). In brief, from the Chandra observation the extended X-ray emission from the lensing cluster was detected out to $$\sim$$1$$^\prime\!\!.$$5 from the cluster center, with a temperature of $$\sim$$6.4 keV. Assuming isothermal profiles, Ota et al. (2006) constrained the projected mass profile and argued that the mass within 100 kpc agrees well with the mass expected from strong lensing. Here, we compare our result of new mass modeling with the X-ray result.

Figure 3 compares the projected two-dimensional mass profiles from gravitational lensing and X-ray measurements. We confirmed that the profiles agree quite well with each other, including the radial slopes of the profiles. This agreement suggests that the effect of the halo triaxiality, which affects the apparent two-dimensional lensing masses particularly near the center of the cluster (see, e.g., Oguri et al. 2005; Gavazzi 2005), is not significant.

However, it should be noted that the best-fit halo mass of $$M_{\rm vir} =$$ 1.0 $$\times$$ 10$$^{15}\ h^{-1}\ M_\odot$$ and the concentration of $$c_{-2} =$$ 2.8 are quite different from those inferred from X-ray, $$M_{\rm vir} \sim$$ 4.3 $$\times$$ 10$$^{14}\ h^{-1}M_\odot$$ and $$c_{-2} \sim$$ 6.1. One reason for this is the strong degeneracy between $$M_{\rm vir}$$ and $$c_{-2}$$ inherent to strong-lens mass modeling. Basically, strong lenses constrain the central core mass of the cluster, which is a strong function of both $$M_{\rm vir}$$ and $$c_{-2}$$. The determinations of $$M_{\rm vir}$$ and $$c_{-2}$$ solely from strong lensing relies on extrapolation of the subtle change of the radial slope out to much larger radii. The robust determination of these parameters from lensing therefore requires addition constraints from weak gravitational lensing.

### Generalized NFW Profile

We now allow the inner slope of the dark-halo component, $$\alpha$$, to vary to see how well the current strong lens data, including the time delays, which are quite helpful to break any degeneracy in the mass models (e.g., Kawano & Oguri 2006), can constrain the inner-density profile. Specifically, for each fixed value of $$\alpha$$ we optimize the other parameters. Figure 4 shows the $$\chi ^2$$ difference as a function of $$\alpha$$. We find that our mass modeling is quite consistent with the NFW profile, i.e., $$\alpha =$$ 1. We constrain the range of the slope to 0.76 $$\lt \alpha \lt$$ 1.41 at the 95% confidence limit. A profile as steep as $$\alpha =$$ 1.5 is clearly rejected. There is a clear correlation between $$\alpha$$ and the velocity dispersion of galaxy G1, such that the best-fit velocity dispersion decreases with increasing $$\alpha$$, which approximately conserves the central core mass of the total matter distribution.

As discussed in subsection 3.2, the Chandra X-ray observation suggests that the lensing cluster may have a larger value of the concentration parameter than the best-fit NFW model predicts. We include this effect by adding a prior of $$c_{-2} =$$ 6$$\ \pm\ $$1.5 to the mass model to see how the constraint on $$\alpha$$ is modified. The result shown in figure 4 indicates that the constraint is basically shifted to a lower $$\alpha$$, i.e., a shallower inner-density slope. The resulting range is 0.62 $$\lt \alpha \lt$$ 1.14 at the 95% confidence limit.

In figure 4, we also show the total magnification factor for the five quasar images, $$\mu_{\rm tot}$$, and the time delay between quasar images A and D, $$\Delta t_{\rm AD}$$, predicted by the best-fit model for each fixed $$\alpha$$. We find that $$\mu_{\rm tot}$$ is decreasing and $$\Delta t_{\rm AD}$$ is increasing with increasing $$\alpha$$, which are consistent with the well-known dependence of the magnification and time delay on the radial density slope (e.g., Wambsganss & Paczyński 1994; Oguri & Kawano 2003). Thus, the reported lower limit of $$\Delta t_{\rm AD} \gt $$ 1250 d (Fohlmeister et al. 2008) prefers a steeper inner slope (larger $$\alpha$$), which appears to be opposed to the effect of the prior from X-ray discussed above. In either case, the strong dependence of $$\Delta t_{\rm AD}$$ on $$\alpha$$ implies that additional A–D time delay measurements will greatly help to constrain the mass model further.

## Summary

We have revisited the parametric mass modeling of SDSS J1004$$+$$4112, a unique quasar–cluster lens system, using newly developed mass modeling software. We included several new constraints, including the positions of spectroscopically confirmed multiply imaged galaxies, time delays between some quasar images, and faint central images. Our model, comprising a halo component modeled by the generalized NFW profile, member galaxies including the brightest cluster galaxy G1, and perturbation terms, very successfully reproduced all observations including time delays.

Unlike earlier claims based on parametric mass modeling, we have found that the center and orientation of the dark-halo component is in good agreement with those of member galaxies and the Chandra X-ray observation, implying that the cluster is highly relaxed. The radial profile from strong lensing is also in excellent agreement with the mass profile inferred from the X-ray observation. Our mass modeling prefers the dark-halo component with the inner slope close to $$\alpha =$$ 1, being consistent with the so-called NFW density profile. Additional measurements of the time delay between quasar image A and D will be useful to constrain the mass model further. The predicted total magnification of $$\mu_{\rm tot} =$$ 67 for the NFW profile is quite large compared with those for typical galaxy-scale lenses, because of the shallower density profiles for clusters. This makes the lens system a good site for studying the fine structure of the quasar through microlensing (Richards et al. 2004; Green 2006; Lamer et al. 2006; Pooley et al. 2007), or for studying the host galaxies of distant quasars (Ross et al. 2009). We note that our result based on parametric mass modeling is broadly consistent with recent non-parametric mass modeling by e.g., Liesenborgs et al. (2009).

Our result suggests that the core of the lensing cluster at $$z =$$ 0.68 is highly evolved. Recently, Limousin et al. (2010) showed that the cluster MACS J1423.8$$+$$2404 at $$z =$$ 0.54 is similarly highly relaxed based on comparisons of the mass, light, and gas distributions. Therefore our result may point to the fact that relaxed clusters are quite common already at $$z \sim$$ 0.6.

It is clear that the lensing cluster SDSS J1004$$+$$4112 is currently one of the best-studied high-redshift clusters whose inner density profile is very tightly constrained by strong lensing. Additional constraints on this cluster by weak lensing, the Sunyaev-Zel’dovich effect, and spectroscopic identifications of member galaxies will be important to advance our understanding of high-redshift clusters.

We would like to thank the referee, Marceau Limousin, for useful comments and suggestion.

## Lens Modeling Software *glafic*

*glafic*

We have developed a comprehensive software package, called *glafic*,^{2} that can be used for a wide variety of analysis for gravitational lensing. Its features include efficient computations of lensed images for both point and extended sources, the handling of multiple sources, a wide range of lens potentials, and the implementation of a noble technique for mass modeling.

In this code, we adopt the adaptive mesh refinement algorithm described in Keeton (2001b) for deriving lensed point-source images, although the code uses rectangular coordinates rather than polar coordinates. Critical curves are computed using a marching squares technique (Jullo et al. 2007). In simulations of extended sources, one can convolve point-spread functions and include a number of galaxies read from a catalog file, which should also be useful for weak-lensing work. For more details, readers are referred to a manual available at the *glafic*’s website.^{2}

## Source-plane χ² Minimization

Strong-lens modeling with the standard $$\chi ^2$$ minimization is sometimes time-consuming, especially when many lens potential components and images are involved. One way to overcome this problem is to evaluate $$\chi ^2$$ in the source plane instead of the image plane. Although the source plane $$\chi ^2$$ involves approximations (given that observational measurements are always made in the image plane) and is therefore less accurate than the image plane $$\chi ^2$$, it allows one to estimate $$\chi ^2$$ without solving the nonlinear lens equation. This technique has been adopted by several authors (e.g., Kayser 1990; Kochanek 1991; Koopmans et al. 1998; Keeton 2001b; Smith et al. 2005; Bradač et al. 2005; Halkola et al. 2006; Jullo et al. 2007; Sand et al. 2008), although the implementations were quite different for different papers. Here, we describe our implementation and argue about the accuracy of the technique.

For a given source position $$\boldsymbol{u}, \chi ^2$$ is estimated as

Here, $$\phi(\boldsymbol{x})$$ is the so-called lens potential. From the lens equation, the magnification factor, $$\mu_i$$, is computed as

Finally, the time delay $$\Delta t_i$$ is

Computing $$\chi ^2$$ in equation (A1) usually requires the derivation of $$\boldsymbol{x}_i$$ for a given $$\boldsymbol{u}$$, which is the most time-consuming part because the lens equation is *not* one-to-one mapping, and thus extensive solution finding in the image plane is needed. In the source plane $$\chi ^2$$ technique, we do not solve the lens equation, but just approximate $$\chi ^2$$ as follows. Assuming $$\boldsymbol{x}_i$$ and $$\boldsymbol{x}_{i,{\rm obs}}$$ are close with each other, we can write

Similarly, the magnification and time delay are approximated as

Note that $$m_0$$ and $$\Delta t_0$$ are nuisance parameters whose best-fit values can easily be derived as

Also note that one can adopt image fluxes as constraints rather than magnitudes. The modification for this is quite straightforward (see also Keeton 2001b).

We find that the source plane $$\chi ^2$$, if properly evaluated like above, is sufficiently accurate to be used in most cases of strong-lens mass modeling. As a specific example, we consider a mass model consisting of an SIE and external shear. The model parameters are $$\sigma =$$ 320 km s$$^{-1}, e =$$ 0.3, $$\theta_{\rm e} =$$ 20$$^\circ, \gamma =$$ 0.1, $$\theta_\gamma =$$ 90$$^\circ, z_{\rm l} =$$ 0.5, $$z_{\rm s} =$$ 2.0, $$\boldsymbol{u} =$$ (0.04, $$-$$0.02). We then add errors to the image positions, magnitudes, time delays, and fit the “observed” four-image system using the same mass model. We adopted a positional error of 0$$^{\prime\prime}\!\!\!.$$01, a magnitude error of 0.1, and a (relative) time delay error of 0.5 d. We ran a Markov Chain Monte Carlo (MCMC) for this model to explore $$\chi ^2$$ around the best-fit model, and at each step we computed both $$\chi^2_{\rm img}$$ and $$\chi^2_{\rm src}$$ to see their difference. The result shown in figure 5 indicates that $$\chi^2_{\rm src}$$ agrees well with $$\chi^2_{\rm img}$$ within a few percent level, which are sufficiently accurate to derive the best-fit mass model and errors on the best-fit parameters. We note that $$\chi^2_{\rm src}$$ is even more accurate when observational constraints are tighter.

We note that Halkola et al. (2006) also studied the accuracy of the source plane $$\chi ^2$$ based on modeling of the strong-lens cluster A1689. They found a clear correlation between $$\chi^2_{\rm img}$$ and $$\chi^2_{\rm src}$$, but with different values. This is because they did not take account of the full magnification tensor [equation (A7)], but simply adopted a magnification factor to approximate $$\chi ^2$$, as has often been done in the literature. Thus, our result highlights the quantitative importance of the proper approximation taking the full lens mapping into account.