A Mass Model for the Lensing Cluster SDSS J1004+4112: Constraints from the Third Time Delay

The lens system SDSS J1004+4112 was the first example of a quasar lensed by a cluster (Inada et al. 2003). The source and lens redshifts are z_s = 1.73 and z_l = 0.68, respectively. This leads to an image separation of 15'' and large magnified images of the host galaxy. It is the second largest quasar lens after SDSS J1029+2623 (Inada et al. 2006). Since its discovery, many additional observational constraints have been obtained: more galaxy cluster members (Oguri et al. 2004), the central quasar image (Inada et al. 2005, 2008), other background lensed sources (Sharon et al. 2005; Liesenborgs et al. 2009; Oguri 2010), the velocity dispersion of the brightest galaxy cluster (Inada et al. 2008), and the time delays between quasar images A, B, and C (Fohlmeister et al. 2007, 2008).

As new observational constraints were measured, increasingly sophisticated mass models of the lens cluster were made using both parametric (Inada et al. 2003; Oguri et al. 2004; Fohlmeister et al. 2007; Oguri 2010) and nonparametric (Williams & Saha 2004; Saha et al. 2006; Liesenborgs et al. 2009; Mohammed et al. 2015) approaches. Recently, the time delay of the fourth quasar image (Muñoz et al. 2022) has been measured at 2457 days between images D and C, and it is over one year longer than predicted (Oguri 2010). Indeed, the ability of the models to predict the unmeasured time delays has been remarkably poor.

Apart from the models of the mass distribution, this system has been widely studied not only in the optical but also from X-rays to radio (Lamer et al. 2006; Ota et al. 2006; Ross et al. 2009; Jackson 2011; Chen et al. 2012; McKean et al. 2021). However, there are still issues that remain unsolved like the origin of the enhanced asymmetric wings in the broad emission lines of quasar image A (Richards et al. 2004; Gómez-Álvarez et al. 2006; Lamer et al. 2006; Motta et al. 2012; Fian et al. 2018; Popović 2020; Fian et al. 2021).

Here we revisit the mass model making use of the newly measured time delay. The new delay should help to constrain the inner mass profile of the cluster halo and to look for deviations from a Navarro–Frenk–White (NFW) profile (Kawano & Oguri 2006; Oguri 2010). The paper is organized as follows. In Section 2 we present the available observational constraints. In Section 3 we describe the modeling method and the components of the mass model. Section 4 presents the results and they are discussed in Section 5. Throughout this paper we assume a flat ΛCDM cosmology with Ω_M = 0.26, Ω_Λ = 0.74 and H₀ = 72 km s⁻¹ Mpc⁻¹.

In addition to the lensed quasar, SDSS J1004+4112 has seven lensed background galaxies at three different redshifts. Their redshifts, image positions, and magnitudes (with respect to quasar image A) are from Oguri (2010) and listed in Table 1. For the quasar time delays we use the results of Muñoz et al. (2022) with the uncertainties set to five times the formal errors (see Table 1) since time delay uncertainties are almost always underestimated (see, e.g., Tie & Kochanek 2018). We also adopted the positions, ellipticities, position angles, and luminosity ratios (with respect to the central galaxy) of 14 cluster galaxy members (see Table 2) from Oguri (2010). For the brightest cluster galaxy (BCG) we adopt a position of (7114,4409), an ellipticity of e = 0.30 ± 0.05, a major axis PA of θ_e = 152 ± 5° (Oguri 2010), and a central velocity dispersion of 352 ± 13 km s⁻¹ (Inada et al. 2008).

We modeled the mass distribution of the lens using both lensmodel (Keeton 2001) and glafic (Oguri 2010). Both packages employ parametrized mass profiles and optimize the model using the downhill simplex method (see, e.g., Press et al. 1992). The final results are computed with glafic in order to ease the comparison with Oguri (2010). The different mass components and any constraints on their properties are described in the following subsections.

3.1. Dark Matter Halo

We modeled the dark matter (DM) mass distribution as a generalized Navarro–Frenk–White (gNFW) profile (Navarro et al. 1997; Jing & Suto 2000) with a 3D density profile

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{\rho }_{s}}{{\left(r/{r}_{s}\right)}^{\alpha }{\left(1+{{rr}}_{s}\right)}^{3-\alpha }}.\end{eqnarray} \tag{ 1 }$$

With the addition of its 2D ellipticity, the projected profile is completely described by seven parameters: M_vir (virial mass), x and y (position), e (ellipticity), θ_e (position angle), c₋₂ (concentration parameter), and α (inner slope). The relationships between (ρ_r , r_s ) and (M_vir, c₋₂) can be found in glafic's manual. ⁸ We do not impose any restrictions on the parameters, however the inner slope α is set to unity until the whole model is optimized.

3.2. BCG

The BCG is parametrized by a pseudo-Jaffe (pJaffe) profile

$$\begin{eqnarray}&&\kappa =2\pi {\left(\displaystyle \frac{\sigma }{c}\right)}^{2}\displaystyle \frac{{D}_{\mathrm{ls}}}{{D}_{\mathrm{os}}}\left[\displaystyle \frac{1}{R}-\displaystyle \frac{1}{\sqrt{{r}_{\mathrm{trun}}^{2}+{R}^{2}}}\right].\end{eqnarray} \tag{ 2 }$$

This model is defined by six parameters when ellipticity is added: σ (velocity dispersion), x and y (position), e (ellipticity), θ_e (position angle), and r_trun (truncation radius, where beyond this scale, the convergence falls as R⁻³). We fix the position to the observed value and include Gaussian priors on e and θ_e . The observed velocity dispersion does not directly correspond to the model dispersion and the Appendix details how we use the measured dispersion as a constraint to obtain a model velocity dispersion prior of 325 ± 20 km s⁻¹. The truncation radius is loosely constrained using the same prior as Oguri (2010) (r_trun = 8'' ± 4'') because of the observed correlation between the velocity dispersion and truncation radius (Natarajan et al. 2009).

3.3. Other Cluster Galaxy Members

The rest of the cluster galaxies are modeled by scaled pseudo-Jaffe ellipsoids (mass profile "gals" in the glafic software). Their σ and r_trun are scaled relative to the luminosity of the BCG as

$$\begin{eqnarray}&&\sigma ={\sigma }_{\ast }{\left(\displaystyle \frac{L}{{L}_{\ast }}\right)}^{1/4}\,{\rm{a}}{\rm{n}}{\rm{d}}\,{r}_{{\rm{t}}{\rm{r}}{\rm{u}}{\rm{n}}}={r}_{\ast }{\left(\displaystyle \frac{L}{{L}_{\ast }}\right)}^{1/2},\end{eqnarray} \tag{ 3 }$$

where σ_* and r_* are model parameters. The ratio L/L_*, ellipticity, and position angle of each galaxy member are constrained by the observed values given in Table 2. Since the reference galaxy is the BCG, the scale parameters σ_* and r_* are initially set equal to the velocity dispersion and truncation radius of the BCG, but, during the optimization, both are allowed to vary.

3.4. External Perturbations

We also include multipole perturbations to mimic deviations of the DM halo from a perfect ellipsoid. Like Oguri (2010) we used multipoles of order m = 2, 3, 4, and 5 (mpoles model) with the potential

$$\begin{eqnarray}&&\phi =-\displaystyle \frac{\epsilon }{m}{r}^{2}\cos m(\theta -{\theta }_{\epsilon }-\pi /2).\end{eqnarray} \tag{ 4 }$$

For m = 2, the perturbation is an external shear (named "pert" in glafic) so we denote the m = 2 parameters as γ and θ_γ instead. We center these potentials on the position of the DM halo.

3.5. Background Sources

The optimization was performed as follows. First we optimized all the parameters fitting the images on the source plane. Once we found the minimum χ², we ran the optimization fitting the images on the image plane. Next we computed the parameter uncertainties fitting the images on the source plane in order to maintain a sensible computing time. Finally, we ran the optimization on the image plane using as initial parameters the model that produced the minimum χ² in the uncertainties estimation step.

The central values and the 1σ and 2σ (corresponding to Δχ² = 1 and Δχ² = 4, respectively) uncertainties are given in Table 3. We define the uncertainties using the model most distant from the χ² minimum that is within the Δχ² limit. The χ² contours and profiles for the mass M_vir, concentration c₋₂, and inner slope α of the gNFW model are displayed in Figure 1. The parameters M_vir and c₋₂ are strongly correlated because they are defined on scales much larger than the Einstein ring, an issue we discuss in more detail below. The parameters do not have smooth χ² profiles because the minimizer has difficulties in finding the true minimum.

Zoom In Zoom Out Reset image size

Table 3. Central Values and Errors at 1σ and 2σ for the Mass Model Parameters

Model	Parameters		1σ	2σ	Best Fit
gNFW	Mvir (1015 h−1 M☉)	1.0	${}_{-0.2}^{+0.3}$	${}_{-0.4}^{+0.7}$	1.0
	x ('')	7.04	${}_{-0.07}^{+0.07}$	${}_{-0.14}^{+0.13}$	7.03
	y ('')	4.95	${}_{-0.10}^{+0.09}$	${}_{-0.19}^{+0.2}$	4.88
	e	0.36	${}_{-0.03}^{+0.02}$	${}_{-0.05}^{+0.04}$	0.36
	θe (deg)	158.6	${}_{-1.8}^{+1.5}$	${}_{-3}^{+3}$	158.1
	c−2	2.3	${}_{-0.8}^{+0.8}$	${}_{-1.7}^{+1.7}$	2.2
	α	1.18	${}_{-0.03}^{+0.02}$	${}_{-0.18}^{+0.11}$	1.20
pJaffe	σ (km s−1)	296	${}_{-4}^{+5}$	${}_{-12}^{+10}$	290
	e	0.40	${}_{-0.04}^{+0.06}$	${}_{-0.10}^{+0.12}$	0.37
	θe (deg)	152	${}_{-3}^{+4}$	${}_{-7}^{+9}$	154
	rtrun ('')	10	${}_{-2}^{+2}$	${}_{-4}^{+6}$	10
gals	σ* (km s−1)	351	${}_{-10}^{+12}$	${}_{-20}^{+20}$	337
	r* ('')	9	${}_{-2}^{+3}$	${}_{-4}^{+6}$	9
pert	γ (10−2)	5.8	${}_{-1.2}^{+1.0}$	${}_{-2}^{+2}$	5.7
	θγ (deg)	66	${}_{-5}^{+3}$	${}_{-9}^{+7}$	65
mpole3	(10−2)	1.6	${}_{-0.3}^{+0.3}$	${}_{-0.6}^{+0.6}$	1.4
	θ (deg)	−9	${}_{-3}^{+2}$	${}_{-7}^{+6}$	−7
mpole4	(10−2)	1.1	${}_{-0.2}^{+0.2}$	${}_{-0.4}^{+0.5}$	1.1
	θ (deg)	40	${}_{-4}^{+3}$	${}_{-7}^{+8}$	41
mpole5	(10−2)	1.39	${}_{-0.2}^{+0.13}$	${}_{-0.3}^{+0.3}$	1.25
	θ (deg)	16.2	${}_{-1.8}^{+2}$	${}_{-4}^{+6}$	17.6

Download table as:  ASCIITypeset image

The optimization on the image plane starting from the parameters which give the minimum χ² in the error computations is considered the best-fit model and it is shown in the last column of Table 3. The central values of the χ² profiles do not exactly correspond to the best-fit parameters (see columns 3 and 6 of Table 3) because the two results come from different optimization procedures. Nonetheless, all the parameters of the best model fit lie between the reported 1σ confidence intervals except for σ and σ_* which lie between their 2σ intervals.

The predicted image positions and critical curves of the four background sources for this model are shown in Figure 2. This best-fit model has a χ² value of 59.91 in the image plane which corresponds to ${\chi }_{\nu }^{2}={\chi }^{2}/\nu =1.50$. The observational constraints on the quasar contribute 11.76 to the total χ² and the A, B, and C groups of galaxies contribute 30.95, 2.92, and 8.16, respectively. Additionally, the χ² associated with the six parameter priors is 6.13. The model fits all three delays well, with a χ² = 1.33 for the delay constraints. The model predicts a time delay between the leading image C and the fifth image E of Δt_EC = 2853.90 days or 7.81 yr, although this delay will be virtually impossible to measure due to the faintness of image E. Additionally, Table 4 gives the convergence κ and shear γ at the quasar positions along with the total magnification.

Zoom In Zoom Out Reset image size

The goodness of the fit for the A group of galaxies is substantially worse than for the rest of the background galaxies because the model predicts seven images instead of the five observationally reported. This occurs because the images A1.2, A2.2, and A3.2 lie very close to a critical curve. A very small change in the source position can change the number of the images and that is precisely what occurred in our model where one image has unfolded into three (they are marked as "unfolded" in Figure 2). The A group images are really composed of extended arcs so the problem would be solved by increasing their position uncertainties, although the model would be more realistic if we modeled the galaxies as extended sources. Despite that, the reconstructed image positions properly reproduce all the main structures of the system (see Figure 3).

Zoom In Zoom Out Reset image size

One of the main results from our model is the value of the inner slope of the generalized NFW profile, $\alpha \,={1.18}_{-0.03(-0.18)}^{+0.02(+0.11)}$ at the 68% (95%) confidence level. We must bear in mind that this value is the asymptotic slope of the DM for r → 0. However at this regime the mass of the BCG becomes important and it is difficult to distinguish both components because the lensing effect is sensitive to the total mass of the system. A different parameterization of the BCG would alter the inner slope value of the generalized NFW profile. If we wanted to determine the mass distribution in the innermost regions of the cluster more precisely, we would probably need a velocity dispersion profile in order to separate the two components better.

Newman et al. (2013a, 2013b) made such models for seven massive galaxy clusters and found very shallow inner slopes of $0.50\pm 0.10{{\rm{(random)}}}_{-0.13}^{+0.14}{\rm{(systematic)}}$. On the other hand, Smith et al. (2017) inferred a DM inner slope for the cluster Abell 1201 of α = 1.01 ± 0.12 using the same procedure but with other assumptions. Schaller et al. (2015) and He et al. (2020) tested this method with EAGLE hydrodynamical simulations and found that the results were dependent on the assumptions about the mass-to-light ratio and the intrinsic degeneracies of the profiles. The simulations suggest that the DM halo distribution is compatible with the original NFW profile and our value. While our model likely has degeneracies between the inner slope of the BCG (which we hold fixed) and the DM halo, our result is compatible with these simulations of massive galaxy clusters.

Our inferred virial mass and concentration parameter for the generalized NFW cluster halo are in good agreement with those obtained with the lens model of Oguri (2010). However they differ from the values reported by Ota et al. (2006) based on X-ray observations, ${M}_{\mathrm{vir}}={4.2}_{-1.5}^{+2.6}\times {10}^{14}\ {h}^{-1}{M}_{\odot }$ and ${c}_{-2}={6.1}_{-1.2}^{+1.5}$. These differences can be explained by the degeneracy between M_vir, c₋₂ and α in strong lensing models pointed out by Oguri (2010) and shown for our model in Figure 1. The source of this degeneracy comes from the fact that the lensing effect is mainly sensitive to the enclosed mass projected inside the region delimited by the lensed images but the parameters M_vir and c₋₂ are defined on the much larger scale of the virial radius and so are extrapolations of our mass model. In this system, the outermost images are ∼125 kpc from the center of the cluster whereas the parameters M_vir and c₋₂ are extrapolations to a radius r_vir ≈ 2 Mpc. It is better to test the compatibility of the two results by comparing the mass enclosed inside a smaller radius to which both methods are sensitive. Ota et al. (2006) reported the projected mass inside 70 h⁻¹kpc to be ${3.5}_{-0.8}^{+1.3}\times {10}^{13}\ {h}^{-1}{M}_{\odot }$, in good agreement with our enclosed mass of 3.7 × 10¹³ h⁻¹ M_☉.

We find an offset between the BCG and the DM halo of ${3.8}_{-0.7(-1.3)}^{+0.6(+1.4)}$ kpc at the 68% (95%) confidence level. Previous parametric models have found similar offsets for this system (Oguri et al. 2004; Fohlmeister et al. 2007) although Oguri (2010) found an offset compatible with zero. From fitting 10,000 galaxy clusters with strong lensing measurements, Zitrin et al. (2012) found a mean separation of ${18}_{-12}^{+37}$ kpc where the uncertainty is the observed scatter. Kluge et al. (2020) estimated the separation between the BCG and the intracluster light tracing the DM distribution for local clusters, and found a mean offset of 36 kpc with a sample scatter of 33 kpc. Our inferred offset is compatible with both results.

We also find that the DM halo and the BCG are virtually aligned and with similar ellipticities. On the other hand, the velocity dispersion scale parameter for the cluster galaxy members differs from the BCG parameterization even though they were started with values scaled to the BCG. This difference suggests that BCGs are structurally different from the other cluster galaxy members because they have more DM content than the rest. Another feature in the BCG velocity dispersion is that the model fits the pseudo-Jaffe velocity dispersion to a smaller value than the prior. That could be explained by the fact that the observed velocity dispersion includes the effect of the DM halo whereas the model velocity dispersion parameter represents only the mass of the BCG.

In conclusion, we have made use of the recently measured time delay of image D and previous observational data to model the lensing cluster SDSS J1004+4112. The model parameters we obtained are broadly consistent with previous estimates but they are now better constrained, in particular the inner slope of the DM halo of the lensing cluster.

This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program GO-9744. J.A.M. and E.M. are supported by the Spanish Ministerio de Ciencia e Innovación with the grants PID2020-118687GB-C32 and PID2020-118687GB-C31. J.A.M. is also supported by the Generalitat Valenciana with the project of excellence Prometeo/2020/085. C.S.K. is supported by NSF grants AST-1908570 and AST-1814440.

The glafic package parametrizes the pseudo-Jaffe profile with a velocity dispersion which corresponds to the central velocity dispersion of a singular isothermal sphere (${\sigma }_{\mathrm{mod}}$). On the other hand, the observed velocity dispersion is the averaged velocity dispersion (σ_obs) along the line of sight inside a finite size slit. The observed velocity dispersion is produced by a combination of the BCG mass and the DM halo, however we assume that the DM halo is negligible for a central velocity dispersion and we only consider the dispersion velocity due to the BCG modeled as a pseudo-Jaffe ellipsoid. In order to obtain their relationship we make use of the expressions of the densities, masses, and velocities dispersion of Elíasdóttir et al. (2007). Additionally, we transform the expressions into dimensionless quantities plus a prefactor in order to simplify the expressions and we end with the following definitions:

• 1.  
  Surface density:

  $$\begin{eqnarray}&&\begin{array}{c}{\rm{\Sigma }}(R)=\displaystyle \frac{{\sigma }_{{\rm{mod}}}^{2}}{2{{Gr}}_{{\rm{trun}}}}\tilde{{\rm{\Sigma }}}(x),\\ {\rm{where}}\,\,\tilde{{\rm{\Sigma }}}(x)=\displaystyle \frac{1}{\sqrt{{f}^{2}+{x}^{2}}}-\displaystyle \frac{1}{\sqrt{1+{x}^{2}}}.\end{array}\end{eqnarray} \tag{ A1 }$$
• 2.  
  Volume density:

  $$\begin{eqnarray}&&\begin{array}{c}\rho (r)=\displaystyle \frac{\sqrt{q}\left(1-{f}^{2}\right){\sigma }_{{\rm{mod}}}^{2}}{2\pi {{Gr}}_{{\rm{trun}}}^{2}}\tilde{\rho }(z),\\ {\rm{where}}\,\,\tilde{\rho }(z)=\displaystyle \frac{1}{\left({f}^{2}+{z}^{2}\right)\left(1+{z}^{2}\right)}.\end{array}\end{eqnarray} \tag{ A2 }$$
• 3.  
  Enclosed projected mass:

  $$\begin{eqnarray}&&\begin{array}{c}{M}_{2{\rm{D}}}(R)=\displaystyle \frac{\pi {r}_{{\rm{trun}}}{\sigma }_{{\rm{mod}}}^{2}}{{Gq}}{\tilde{M}}_{2{\rm{D}}}(x),\\ {\rm{where}}\,\,{\tilde{M}}_{2{\rm{D}}}(x)=\sqrt{{f}^{2}+{x}^{2}}-\sqrt{1+{x}^{2}}+1-f.\end{array}\end{eqnarray} \tag{ A3 }$$
• 4.  
  Enclosed mass:

  $$\begin{eqnarray}\begin{array}{rcl}{M}_{3{\rm{D}}}(r) & = & \displaystyle \frac{2{r}_{{\rm{trun}}}{\sigma }_{{\rm{mod}}}^{2}}{{Gq}}{\tilde{M}}_{3{\rm{D}}}(z),\,\,{\rm{where}}\\ {\tilde{M}}_{3{\rm{D}}}(z) & = & \arctan \left(z\right)-f\arctan \left(z/f\right).\end{array}\end{eqnarray} \tag{ A4 }$$
• 5.  
  Projected velocity dispersion:

  $$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{P}^{2}(R) & = & \displaystyle \frac{4\left(1-{f}^{2}\right){\sigma }_{{\rm{mod}}}^{2}}{\pi \sqrt{q}}\tilde{{\sigma }_{P}^{2}}(x),\,\,{\rm{where}}\\ \tilde{{\sigma }_{P}^{2}}(x) & = & \displaystyle \frac{1}{\tilde{{\rm{\Sigma }}}(x)}{\displaystyle \int }_{x}^{\infty }\displaystyle \frac{{\tilde{M}}_{3{\rm{D}}}(z)\tilde{\rho }(z)}{{z}^{2}}\sqrt{{z}^{2}-{x}^{2}}{dz}.\end{array}\end{eqnarray} \tag{ A5 }$$
• 6.  
  Averaged projected velocity dispersion inside R:

  $$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{\rm{obs}}}^{2} & = & \langle {\sigma }_{P}^{2} \rangle (R)=\displaystyle \frac{4\left(1-{f}^{2}\right){\sigma }_{{\rm{mod}}}^{2}}{\pi \sqrt{q}}\tilde{ \langle {\sigma }_{P}^{2} \rangle }(x),\,\,{\rm{where}}\\ \tilde{ \langle {\sigma }_{P}^{2} \rangle }(x) & = & \displaystyle \frac{1}{{\tilde{M}}_{2{\rm{D}}}(x)}{\displaystyle \int }_{0}^{x}\tilde{{\sigma }_{P}^{2}}({x}^{{\rm{{\prime} }}})\tilde{{\rm{\Sigma }}}({x}^{{\rm{{\prime} }}}){x}^{{\rm{{\prime} }}}{{dx}}^{{\rm{{\prime} }}}.\end{array}\end{eqnarray} \tag{ A6 }$$

The dimensionless variables are $x\equiv R\sqrt{q}/{r}_{\mathrm{trun}}$, $z\,\equiv r\sqrt{q}/{r}_{\mathrm{trun}}$, and $f\equiv {r}_{\mathrm{core}}/{r}_{\mathrm{trun}}$ and q = 1 − e is the axis ratio. In our case, f = 0 and the slit width employed to measure the velocity dispersion was 04 (Inada et al. 2008) so we used the average dispersion inside R = 02. Taking into account the central values and errors for the ellipticity and the truncation radius, we obtain ${\sigma }_{\mathrm{mod},\mathrm{prior}}=325\pm 20$ km s⁻¹ as the prior for the model velocity dispersion.