Issue 
A&A
Volume 691, November 2024



Article Number  A97  
Number of page(s)  7  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202347378  
Published online  31 October 2024 
Microlensing analysis of 14.5year light curves in SDSS J1004+4112: Quasar accretion disk size and intracluster stellar mass fraction
^{1}
Departamento de Astronomía y Astrofísica, Universidad de Valencia, 46100 Burjassot, Valencia, Spain
^{2}
Observatorio Astronómico, Universidad de Valencia, 46980 Paterna, Valencia, Spain
^{3}
Departamento de Física Teórica y del Cosmos, Universidad de Granada, Campus de Fuentenueva, 18071 Granada, Spain
^{4}
Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, 18071 Granada, Spain
^{5}
Instituto de Astrofísica de Canarias, Vía Láctea s/n, La Laguna, 38200 Tenerife, Spain
^{6}
Departamento de Astrofísica, Universidad de la Laguna, La Laguna, 38200 Tenerife, Spain
^{⋆} Corresponding author; raquel.fores@uv.es
Received:
6
July
2023
Accepted:
11
September
2024
Context. The gravitational lens system SDSS J1004+4112 was the first known example of a quasar lensed by a galaxy cluster. The interest in this system has been renewed following the publication of rband light curves spanning 14.5 years and the determination of the time delays between the four brightest quasar images.
Aims. We constrained the quasar accretion disk size and the fraction of the lens mass in stars using the signature of microlensing in the quasar image light curves.
Methods. We built the six possible histograms of microlensing magnitude differences between the four quasar images and compared them with simulated model histograms, using a χ^{2} test to infer the model parameters.
Results. We infer a quasar disk halflight radius of R_{1/2} = (0.70 ± 0.04)R_{E} = (6.4 ± 0.4) √M/0.3M_{⊙} lightdays at 2407 Å in the rest frame and stellar mass fractions at the quasar image positions of α_{A} > 0.059, α_{B} = 0.056^{+0.021}_{0.027}, α_{C} = 0.030^{+0.031}_{0.021}, and α_{D} = 0.072^{+0.034}_{0.016}.
Conclusions. The inferred disk size is broadly compatible with most previous estimates, and the stellar mass fractions are within the expected ranges for galaxy clusters. In the region where image C lies, the stellar mass fraction is compatible with a stellar contribution from the brightest cluster galaxy, galaxy cluster members, and intracluster light, but the values at images B, D, and especially A are slightly larger, possibly suggesting the presence of extra stellar components.
Key words: accretion / accretion disks / gravitation / gravitational lensing: micro / galaxies: clusters: intracluster medium / quasars: individual: SDSS J1004+4112
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Since Chang & Refsdal (1979) proposed that stars located near the light paths of gravitationally lensed quasars may cause flux variations in the quasar images, the field of gravitational microlensing has been developed and refined such that now we are able to estimate the sizes of quasar accretion disks (see e.g. Kochanek 2004; JiménezVicente et al. 2012; Hainline et al. 2013; Blackburne et al. 2015; MacLeod et al. 2015; Muñoz et al. 2016; Morgan et al. 2018; Fian et al. 2021a) and the abundance of microlenses in the lens (see e.g. Mediavilla et al. 2009; JiménezVicente et al. 2015; Awad et al. 2023).
Quasars also exhibit intrinsic flux variability that needs to be removed in order to obtain the fluctuations produced only by microlensing. If the time delay between images is known, we can subtract pairs of light curves shifted by their time delays to obtain light curves of the microlensing magnification differences. We also need sufficiently long light curves to statistically sample the magnifications and infer the system parameters properly. With the recent measurements of 14.5year light curves and time delays for the four brightest images of SDSS J1004+4112 (Muñoz et al. 2022), we are now able to carry out such studies.
SDSS J1004+4112 was the first discovered lens system in which a background quasar is lensed by a galaxy cluster instead of a single galaxy (Inada et al. 2003). This leads to a large 15″ separation between the images, and hence the light from the quasar travels mainly through the cluster instead of passing through an individual galaxy, as is typical for lensed quasars. With this configuration, the impact of microlensing in the quasar was expected to be small. However, shortly after its discovery, Richards et al. (2004) reported microlensing variability in the blue wing of image A in several emission lines corresponding to the broad line region, and this variability continued over the years (GómezÁlvarez et al. 2006; Lamer et al. 2006; Motta et al. 2012; Fian et al. 2018a, 2021b; Popović et al. 2020). The main explanation for these fluctuations is microlensing, but other possibilities have been raised, such as outflows (Green 2006; Abajas et al. 2007; Popović et al. 2020), because the stellar density in galaxy clusters is expected to be insufficient to drive such microlensing variability.
Microlensing variability in the continuum emission of the accretion disk is also observed (Fohlmeister et al. 2008; Chen et al. 2012; Fian et al. 2016). Fohlmeister et al. (2008) and Fian et al. (2016) used the light curves to determine the quasar disk size in the r band. However, Fohlmeister et al. (2008) used only images A and B, while Fian et al. (2016) used all four images but with the predicted time delay for image D from Oguri (2010), which now we know from Muñoz et al. (2022) was 413 days too short.
We used the newly measured time delays along with the light curves from Muñoz et al. (2022) to constrain the quasar disk size, and, given the quality and long coverage of the light curves, we simultaneously attempted to measure the fraction of the mass in stars in the regions where the quasar images are located. The paper is structured as follows. In Sect. 2 we present the data that we used for the analysis, Sect. 3 describes the methodology, the results are presented in Sect. 4, and we discuss and summarise them in Sect. 5.
2. Data
We used the rband light curves of the four brightest quasar images from Muñoz et al. (2022), which span 14.5 years. These light curves were shifted by the time delays Δt_{AC} = 825.99 ± 0.42, Δt_{BC} = 781.92 ± 0.44, and Δt_{DC} = 2456.99 ± 1.11 days, also from Muñoz et al. (2022). We assumed that the weighted means of the infrared (IR) 8.0 μm joint quasar and host (QSO+Host) magnitude measurements of Ross et al. (2009) are a solid proxy for the unmicrolensed flux ratios: m_{A, IR} = 12.262 ± 0.009, m_{B, IR} = 12.627 ± 0.015, m_{C, IR} = 12.906 ± 0.007, and m_{D, IR} = 13.483 ± 0.013. They are compatible with radio flux ratios (Jackson 2020; McKean et al. 2021; Hartley et al. 2021) and with the broad line region cores of Fian et al. (2024). Finally, we used the values of convergence and shear at the image positions from the mass model of ForésToribio et al. (2022) to simulate the magnification patterns.
3. Methodology
3.1. Observed histograms
To reduce the noise, we smoothed the light curves over a window of 10 days. The smoothing window is much smaller than the estimated microlensing sourcecrossing timescale of t_{S} = 0.28 yr (Mosquera & Kochanek 2011). The average of the data points inside the window was weighted by their uncertainties.
The aim was to remove the intrinsic variability of the quasar, so we subtracted the magnitudes of one light curve from those of another after shifting them by their corresponding time delays. In order to do so, we needed to fit one of the light curves to interpolate between points. We computed the six possible differences – A−B, C−B, D−B, A−C, D−C, and A−D – and we chose to fit the light curve that was going to be subtracted. We fitted this light curve with a fifthorder spline, which gives enough flexibility to account for the fluctuations in the light curves without introducing extra structure. We determined the fit to be unreliable if the gap between two points is larger than 100 days, and we also manually removed the peaks around 5200 and 5600 days (heliocentric Julian date of observation minus 2 450 000) in the light curve of B, which give residuals that can be associated with intrinsic variability. We show the spline fits of light curves B, C, and D in Fig. 1.
Fig. 1. Fits to the rband light curves B, C, and D from Muñoz et al. (2022) using fifthorder splines (solid black curves). Since gaps longer than 100 days were not considered in the analysis, the splines in these regions are not displayed. The days on the xaxis correspond to the observation dates. 
The differences between the shifted pairs of light curves still included the local macromagnification of the lens system, which had to be removed. We took the 8.0 μm measurements of Ross et al. (2009) as the macromagnification baseline assuming that the emitting region of QSO+Host in this band is large enough to not be affected by microlensing. We computed the microlensing differences as
$$\begin{array}{c}\hfill \mathrm{\Delta}{\mu}_{\mathit{XR}}^{i}=({m}_{X}^{i}{m}_{X,\mathrm{IR}})({m}_{R,\mathrm{fit}}^{i}{m}_{R,\mathrm{IR}}),\end{array}$$(1)
where m_{X}^{i} is an epoch of the X light curve that overlaps with the permitted regions of the R (reference) light curve, m_{R, fit}^{i} is the value of the fitted reference light curve shifted by the corresponding time delay between images, and m_{X, IR} and m_{R, IR} are the IR magnitudes. The six possible microlensing differences generated by this procedure are shown in Fig. 2 with the corresponding mean and standard deviation. The light curve with the largest variability is A−B, as expected given the previous detections of microlensing in image A. In the A−B and A−C light curves, a clear microlensing event can be seen between days 6500 and 7500.
Fig. 2. Microlensing differences between the six possible combinations of image pairs. The solid black lines are the mean of the micromagnification, and the dashed lines correspond to the standard deviation. All light curves are shifted to match the time reference of the leading image, C. 
We constructed the observed histograms, h_{XR}, from these light curves taking the observational errors of each point into account via Monte Carlo sampling, assuming that the errors of the data points are Gaussian. We generated 10^{6} histograms for each microlensing difference light curve. The histograms have a bin width of 0.05 mag and range from −5 to 5 mag. The value of each bin is the mean of the counts falling into each bin, and the associated error is the standard deviation; these histograms are shown in Fig. 3.
Fig. 3. Observed histograms of the six possible microlensing differences obtained after averaging the 10^{6} realisations from the Monte Carlo sampling. The error bars are the standard deviation of these realisations. The solid black lines represent the best model difference histograms, and the thin dotted lines are other model histograms used for the parameter inference. 
3.2. Model histograms
We computed magnification maps for the four images using the fast multipole methodinverse polygon mapping^{1} developed by JiménezVicente & Mediavilla (2022). We created 30 × 30 R_{E} magnification maps using stars of a constant mass of 0.3 M_{⊙}. For this mass, the Einstein radius in the source plane corresponds to 9.07 lightdays. Assuming an effective source velocity of 0.106 R_{E}/year (Mosquera & Kochanek 2011), the magnification map resolution was set to 1 pix = 0.003 R_{E} so that one pixel covers the smoothing window of 10 days applied to the observed light curves. We computed the magnification maps with the values of κ and γ reported in ForésToribio et al. (2022) and a variable fraction of the mass in stars, α = κ_{*}/κ.
Once the magnification maps were created, we convolved them with extended sources of different sizes. We modelled the sources as a Gaussian of width r_{s} since it is commonly accepted that the microlensing variability depends mainly on the halflight radius of the source rather than the particular source shape (Mortonson et al. 2005; Muñoz et al. 2016; Vernardos & Tsagkatakis 2019). After the convolution, the magnification maps were normalised in flux to remove the macromagnification, as also done for the data.
According to Mosquera & Kochanek (2011), the effective velocity on the source plane is dominated by the random velocity of the microlenses in the galaxy cluster (517 km/s) rather than by the peculiar velocities of the source, the lens, and the observer (82 km/s, 196 km/s, and 99 km/s, respectively). Hence, the orientation of the paths that the source describes over the magnifications maps is essentially uncorrelated between image pairs (see the discussion in Poindexter & Kochanek 2010). Given the lack of significant correlations, we generated randomly oriented straight tracks with the same length as the observed difference light curves for each image pair on its corresponding magnification map. Since the observed histograms may arise from regions with different mean magnifications, we only kept the pairs of tracks whose average after subtracting one from the other differs by less than 0.1 mag from the mean micromagnification difference of the observed histograms. This 0.1 mag is the error assumed for the unmicrolensed flux ratios given the average error budget in IR found by Ross et al. (2009) and in broad emission line cores by Fian et al. (2024) for image pairs. We collected 500 pairs of tracks that fulfil this condition (see Fig. 5 as an example for the A−B pair) and obtained their microlensing difference histograms with the same bounds and widths as the observed histograms^{2}. Then, a single model histogram for each combination of source size and stellar mass fraction was computed by averaging the histogram bins of the 500 selected pair tracks to obtain
$$\begin{array}{c}\hfill {\stackrel{\sim}{h}}_{\mathit{XR}}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R})=\frac{1}{N}{\displaystyle \sum _{n=1}^{N}}{\stackrel{\sim}{h}}_{\mathit{XR}}^{n}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R}),\end{array}$$(2)
where ${\stackrel{\sim}{h}}_{\mathit{XR}}^{n}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R})$ is the ith bin of the histogram obtained from the nth selected track pair of the difference X − R for a source size r_{s} and stellar mass fractions α_{X} and α_{R}. The term ${\stackrel{\sim}{h}}_{\mathit{XR}}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R})$ is the ith bin of the model microlensing difference histogram for these parameters, which is the average of N = 500 histograms.
3.3. Statistical test
For each pair of images, we compared the observed microlensing histograms with the model histograms that arise from the possible combinations of source sizes, r_{s}, and stellar mass fractions, α, after normalising their counts. The likelihood that a given model with r_{s}, α_{A}, α_{B}, α_{C}, and α_{D} reproduces the observed microlensing differences was computed using the χ^{2} test defined as
$$\begin{array}{c}\hfill {\chi}^{2}={\displaystyle \sum _{\mu}}{\displaystyle \sum _{i}}{\left(\frac{{h}_{\mu}(i){\stackrel{\sim}{h}}_{\mu}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R})}{{\u03f5}_{\mu}(i)}\right)}^{2}.\end{array}$$(3)
The sum in i runs over the histogram bins, and the sum in μ runs over the six distinct image pairs (i.e. the difference between image X and the image of reference, R). The h_{μ}(i) is the ith bin of the normalised observed histogram built as described in Sect. 3.1, ϵ_{μ}(i) is the error associated with that bin, and ${\stackrel{\sim}{h}}_{\mu}(i\u037e{r}_{s},{\alpha}_{X},{\alpha}_{R})$ is ith bin of the model histogram for the given set of parameters computed according to the procedure in Sect. 3.2. Given that lowcount bins can affect the χ^{2} statistic, we only considered bins of the observed histograms with at least five counts (see e.g. James 2006).
4. Results
We computed the χ^{2} test for all the parameter combinations within the following ranges: We varied the source size, r_{s}, from 3.2 lightdays to 7.6 lightdays using logarithmically spaced values and then added more points near the χ^{2} minimum. The stellar mass fractions at the locations of images A, B, C, and D were varied in a range from 0 to 0.2. This range was set to not exceed the average stellar mass fraction found for individual galaxies (e.g. JiménezVicente et al. 2015). Figure 4 shows the 68% and 95% maximum likelihood contours for two parameters (Δχ^{2} = 2.30 and 6.18) and the probability distributions for each variable with their 68% confidence limits (Δχ^{2} = 1). In this figure, a mild Gaussian smoothing is applied to obtain cleaner contours. The maximum likelihood histograms fit the observational data with a reduced χ^{2} of 0.997 (see the solid black lines in Figure 3). The radius is presented as the logarithm of the halflight radius in Einstein radii, which is related to the Gaussian width as R_{1/2} = 1.18 r_{s}.
Fig. 4. Joint Δχ^{2} = χ^{2} − χ_{min}^{2} contours for all pairs of parameters (R_{1/2}, α_{A}, α_{B}, α_{C}, and α_{D}). In the 2D likelihoods, the contours are drawn at the 68% (Δχ^{2} = 2.30) and 95% (Δχ^{2} = 6.18) confidence limits for two parameters. In the 1D distributions, the shaded regions are the 68% (Δχ^{2} = 1.00) confidence regions for one parameter. The maximum likelihood model is marked with a cross and a dashed line for the 2D and 1D plots, respectively. 
With this procedure, we obtain an accretion disk size of ${R}_{1/2}=(0.70\pm 0.04)\phantom{\rule{0.166667em}{0ex}}{R}_{\mathrm{E}}=(6.4\pm 0.4)\sqrt{M/0.3\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}}$ lightdays at 2407 Å in the rest frame. In Table 1 we summarise the stellar mass fractions and the convergence in stars at the four quasar positions. For α_{A}, we can only report a lower limit because of the poor convergence of the likelihood towards larger values, although these large values are unphysical. The magnification maps of images A and B along with the selected 500 pairs of tracks for the maximum likelihood model are shown in Fig. 5 as an example of the method.
Fig. 5. Magnification maps normalised to the mean macromagnification, ⟨μ⟩, for images A and B with the 500 pairs of tracks of the best fitting model. These tracks are randomly oriented, and their microlensing difference histograms were averaged to build the model histogram. The radius of the red circle at the bottom of each map corresponds to the inferred halflight radius of the source (i.e. 0.7 R_{E} or equivalently 6.4 lightdays). 
Inferred stellar mass fraction and convergence in stars at the locations of the four quasar images.
5. Discussion and conclusions
We compare our quasar accretion disk size estimate with previous determinations in Fig. 6. From spectroscopic measurements of images A and B, Fian et al. (2024) inferred a halflight radius of $7.{1}_{3.7}^{+7.4}$ lightdays for the continuum emission corresponding to the rband wavelengths, which fully covers our size estimate. Also from spectroscopic data, Hutsemékers et al. (2023) reported an upper limit on the continuum disk size of R_{1/2} < 2.9 lightdays for ⟨M⟩ = 0.3 M_{⊙} at 2400 Å, roughly the rest frame centre of the r band, which is significantly smaller than our size estimate. Fian et al. (2016) used light curves of the four quasar images and inferred a halflight radius of ${R}_{1/2}=8.{7}_{5.5}^{+18.5}$ lightdays using the histogram product method and ${R}_{1/2}=4.{2}_{2.2}^{+3.2}$ lightdays using the corrected χ^{2} method. They assumed in both cases a stellar mass of 0.3 M_{⊙} and that 10% of the mass is in stars for all images. Their results are totally compatible with our estimate, but we recover tighter constraints because of the longer light curves. JiménezVicente et al. (2014) estimated a Bayesian rband accretion disk size of ${R}_{1/2}=4.{9}_{3.3}^{+4.8}$ lightdays for ⟨M⟩ = 0.3 M_{⊙} using a wavelength dependence of R ∝ λ^{1.3 ± 0.6}, which is inferred from multiwavelength microlensing of images A and B. This value is compatible within errors with the size inferred in this work. Motta et al. (2012) inferred a disk size at λ_{rest} = 3363 Å of ${r}_{s}={6}_{3}^{+4}$ lightdays and a wavelength dependence of R ∝ λ^{1.0 ± 0.4}. This implies a halflight radius in the r band of ${R}_{1/2}=2.{8}_{1.4}^{+1.9}$ lightdays for ⟨M⟩ = 0.3 M_{⊙}, which is in a 1.9σ tension with our estimate. We infer significantly larger sizes than Mosquera & Kochanek (2011) and Fohlmeister et al. (2008), who respectively found R_{1/2} = 0.35 lightdays for a faceon disk in the r band based on the thindisk scaling relation between flux and size and R_{1/2} = 0.6 ± 0.4 lightdays based on the microlensing in images A and B.
Fig. 6. Comparison with previous quasar accretion disk size estimates for SDSS J1004+4112. The results are reported as the rband (λ_{obs} = 6580 Å, λ_{rest} = 2407 Å) halflight radius in lightdays for a mean stellar mass of 0.3 M_{⊙} to homogenise the comparison. 
The Xray monitoring of the Fe Kα lines in SDSS J1004+4112 and another lens systems (Chen et al. 2012; Chartas et al. 2017; Dai et al. 2019) has enabled the size inference of this emission region via microlensing. Chartas et al. (2017) placed an order of magnitude of ∼20 gravitational radii (r_{g}) for the emission region by combining three gravitational lens systems. Dai et al. (2019) constrained the Fe Kα emission region to 5.9 − 7.4 r_{g} at the 68% confidence level for a joint sample of four lenses, including SDSS J1004+4112. These size estimates for the SDSS J1004+4112 black hole mass of 3.98 × 10^{8} M_{⊙} (Chartas et al. 2017) correspond to ∼0.45 and 0.13 − 0.17 lightdays, respectively. Future individual studies of Xray emission lines and continuum emission in other optical bands will be useful to constrain the structure of this quasar accretion disk. Nonetheless, these Xray determinations indicate that the source Xrayemitting region is more compact than the rbandemitting region, as expected.
It is not straightforward to compare our stellar mass fraction results with other determinations. The stellar mass fraction depends not only on the radial distance from the brightest cluster galaxy (BCG) but also on additional cluster galaxy members along the line of sight where the quasar images are located and on the diffuse matter of the intracluster medium. We can take the baryonic mass fraction of the Universe as a reference for an upper estimate on the mean stellar mass fraction. From Planck Collaboration VI (2020) results, the mean fraction of baryons in the Universe is f_{b} = Ω_{b}/(Ω_{b} + Ω_{c}) = 0.1573 ± 0.0013 (Ω_{b} being the baryon density and Ω_{c} the cold dark matter density, both in terms of the critical density). Except for α_{A}, where we only obtain a lower limit, the α values are below this upper limit (see Table 1). Another upper estimate can be the fraction between the surface density of galaxy cluster members and the total surface density of the cluster at the quasar positions from the ForésToribio et al. (2022) model. These limits are defined as α_{i}^{max} = ∑_{n}κ_{pJaffe, n}(r_{i})/κ_{tot}(r_{i}), where n runs over all the modelled galaxies, including the BCG, and r_{i} is the position of the ith quasar image. This should be a fair stellar mass fraction limit since the galaxies were modelled as pseudoJaffe ellipsoids, which include both the dark matter and baryonic mass in galaxies, and the above definition assumes that all the mass in these galaxies is in stars. The maximum value for each image corresponds to ${\alpha}_{\mathrm{A}}^{\mathrm{max}}=0.10$, ${\alpha}_{\mathrm{B}}^{\mathrm{max}}=0.08$, ${\alpha}_{\mathrm{C}}^{\mathrm{max}}=0.03$, and ${\alpha}_{\mathrm{D}}^{\mathrm{max}}=0.12$. The stellar mass fractions’ lower confidence levels are at least 30% of the α_{i}^{max} estimates. Even though the central values do not surpass the mass fraction in pseudoJaffe components, these percentages are rather high because one would expect, on average, ∼20% of the mass in galaxies to be contained in stars (JiménezVicente et al. 2015). Hence, it is likely that stars in the intracluster medium are contributing to increase the total stellar mass fraction.
A lower bound could be the mean stellar mass fraction of entire galaxy clusters given that in this system the quasar images are located within R < 50 kpc/h of the centre. Following the stellar mass fraction relation of Giodini et al. (2009) obtained from observed galaxy groups and poor clusters, we find ${f}_{500}^{\mathrm{stars}}=0.014\pm 0.002$ for a M_{500} = 9.89 × 10^{14} M_{⊙}/h (derived from the mass model of ForésToribio et al. 2022). Gonzalez et al. (2013) derived a similar relation for the projected stellar mass with respect to M_{500}. Dividing this relation by M_{500}, we find a mean stellar mass fraction of 0.0119 ± 0.0014 for our cluster mass. Kravtsov et al. (2018) added more massive clusters to the Gonzalez et al. (2013) sample and derived a M_{*, tot} − M_{500} that resulted in a stellar mass fraction of 0.0088 ± 0.0012. From TNG100 and TNG300^{3} galaxy group and cluster simulations, Pillepich et al. (2018) derived the same type of relation as above, obtaining a stellar mass fraction of around 0.0149. All these observational or numerical works infer a mean stellar mass fraction of around 1% for entire galaxy clusters. This value is within the confidence interval estimated for α_{C}, but for α_{A}, α_{B}, and α_{D} our estimations are larger.
On the other hand, the radial dependence of the stellar mass due to the BCG and the intracluster light (ICL) offers a minimum stellar mass estimate at the specific image positions. They are usually reported as the enclosed stellar mass within a specified radius. This can be translated to a stellar surface density, and, when divided by the critical surface density of the lens system Σ_{crit} = 3.23 × 10^{9} h M_{⊙} kpc^{−2} and by the total convergence at the quasar positions, an estimate of α for each quasar image can be obtained. Kravtsov et al. (2018) reported the enclosed stellar mass of the BCG plus the ICL in the inner 30, 50, and 70 kpc for nine observed galaxy clusters. DeMaio et al. (2018) also inferred the stellar masses within 10, 50, and 100 kpc for 23 galaxy clusters; for 16 of them they found M_{500} > 10^{14} M_{⊙}, which is more similar to SDSS J1004+4112. Henden et al. (2020) reported the radial stellar surface density profile of galaxy clusters from FABLE^{4} simulations. The resulting mean stellar mass fractions and the standard deviations at the locations of each quasar image from these works are reported in Table 2.
Estimates based on previous works for the stellar mass fraction.
These BCG+ICL estimates are lower limits because the contributions from the galaxy cluster members are excluded. To obtain more realistic estimates, the 20% of the mass fraction from galaxy members (excluding the BCG) in the ForésToribio et al. (2022) model (0.008, 0.001, 0.002, and 0.008, for images A, B, C, and D, respectively) can be added to the values estimated by Kravtsov et al. (2018), DeMaio et al. (2018), and Henden et al. (2020). The average discrepancies for the stellar mass fractions at the positions of images A, B, C, and D are 2.16, 1.47, 0.69, and 1.65σ, respectively. For the discrepancy of image A, we used its central value, 0.090 (where the minimum χ^{2} is found) and the lower error, 0.031 (the difference between the central value and the lower estimate). Hence, the stellar mass fraction at image C is completely consistent with the estimates, at images B and D a slight tension is found, and at image A our inference is 2σ discrepant with these estimates. The Fohlmeister et al. (2008) models also preferred larger values of κ_{*} for images A and B, and the authors suggested that the microlensing was probably due to a satellite galaxy rather than stars in the ICL. Recently, Perera et al. (2024) built freeform and hybrid mass models for SDSS J1004+4112 and found in all models a mass clump of around 2″ located at ∼(5, −1)″ from quasar image A. The nature and properties of this possible clump are still undefined, but it remains present after the matter is redistributed following the procedure described in Liesenborgs et al. (2024). Although this substructure may not be responsible for the slightly larger stellar mass fractions that we estimated, this finding indicates that the mass structure of SDSS J1004+4112 may be complicated and composed of substructures yet to be detected. Additional stellar components may help explain the inferred stellar mass fractions being higher than the fractions expected from galaxy members and from ICL studies, but observational evidence should be provided to confirm this hypothesis. If microlensing is being produced by an undetected individual galaxy member, it must be aligned with the line of sight of images and, hence, the quasar emission can mask it, making it difficult to detect. On the other hand, the stellar mass fraction can be explained if diffuse components (mainly in the form of stars) are present in these regions. The origin of this component could be a substructure in the cluster galaxy or a stellar tidal stream from the infall of satellite galaxies into the BCG.
In this work we attempted to estimate the microlensing parameters using magnification fluctuation statistics. This kind of microlensing study in which the time dependence is not included seems to obtain results consistent with other techniques while keeping the computational cost low (see Fian et al. 2016, 2018b, 2021a). Nonetheless, the temporally correlated microlensing, especially the event detected between days 6500 and 7500 in the A−B and A−C difference light curves of Fig. 2, may be worthy of study in a future work with direct light curve fitting (Kochanek 2004) to further test the validity of the method employed in this work.
In conclusion, we have determined the quasar accretion disk size and the stellar mass fraction at the four brightest quasar image positions in SDSS J1004+4112 using the 14.5year light curves and time delays reported by Muñoz et al. (2022). We were able to constrain, for the first time, the stellar mass fraction at the individual image positions. Our quasar disk size is broadly compatible with previous works, and the stellar mass fractions we find lie between the expected upper and lower bounds extracted from the literature. The stellar mass fraction at the position of image C can be mostly explained by the contributions of the BCG, the ICL, and the known galaxy cluster members. At images B and D, the inferred mass fractions are slightly above expectations, with discrepancies around 1.5σ. For image A, only a lower limit has been obtained, and the central value is in a 2σ tension with the contribution of the known stellar components. These inferred stellar mass fractions suggest that other stellar components are present in the regions where images B, D, and especially image A are located.
These computations were performed at the computing service PROTEUS (https://proteus.ugr.es) of the Carlos I Institute of Theoretical and Computational Physics because of the large computational time required.
Acknowledgments
We thank C. S. Kochanek for his useful comments and suggestions throughout this work. We are also grateful for the anonymous referee’s reviews, which have improved the final version of the article. This research was supported by the grants PID2020118687GBC31, PID2020118687GBC32 and PID2020118687GBC33 financed by the Spanish Ministerio de Ciencia e Innovación through MCIN/AEI/10.13039/501100011033. J.A.M. is also supported by the project of excellence PROMETEO CIPROM/2023/21 of the Conselleria de Educación, Universidades y Empleo (Generalitat Valenciana). J.J.V. is also supported by projects FQM108, P20_00334, and AFQM510UGR20/FEDER, financed by Junta de Andalucía.
References
 Abajas, C., Mediavilla, E., Muñoz, J. A., GómezÁlvarez, P., & GilMerino, R. 2007, ApJ, 658, 748 [Google Scholar]
 Awad, P., Chan, J. H. H., Millon, M., Courbin, F., & Paic, E. 2023, A&A, 673, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blackburne, J. A., Kochanek, C. S., Chen, B., Dai, X., & Chartas, G. 2015, ApJ, 798, 95 [Google Scholar]
 Chang, K., & Refsdal, S. 1979, Nature, 282, 561 [Google Scholar]
 Chartas, G., Krawczynski, H., Zalesky, L., et al. 2017, ApJ, 837, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, B., Dai, X., Kochanek, C. S., et al. 2012, ApJ, 755, 24 [Google Scholar]
 Dai, X., Steele, S., Guerras, E., Morgan, C. W., & Chen, B. 2019, ApJ, 879, 35 [NASA ADS] [CrossRef] [Google Scholar]
 DeMaio, T., Gonzalez, A. H., Zabludoff, A., et al. 2018, MNRAS, 474, 3009 [Google Scholar]
 Fian, C., Mediavilla, E., Hanslmeier, A., et al. 2016, ApJ, 830, 149 [Google Scholar]
 Fian, C., Guerras, E., Mediavilla, E., et al. 2018a, ApJ, 859, 50 [Google Scholar]
 Fian, C., Mediavilla, E., JiménezVicente, J., Muñoz, J. A., & Hanslmeier, A. 2018b, ApJ, 869, 132 [Google Scholar]
 Fian, C., Mediavilla, E., JiménezVicente, J., et al. 2021a, A&A, 654, A70 [Google Scholar]
 Fian, C., Mediavilla, E., Motta, V., et al. 2021b, A&A, 653, A109 [Google Scholar]
 Fian, C., Muñoz, J. A., ForésToribio, R., et al. 2024, A&A, 682, A57 [Google Scholar]
 Fohlmeister, J., Kochanek, C. S., Falco, E. E., Morgan, C. W., & Wambsganss, J. 2008, ApJ, 676, 761 [Google Scholar]
 ForésToribio, R., Muñoz, J. A., Kochanek, C. S., & Mediavilla, E. 2022, ApJ, 937, 35 [Google Scholar]
 Giodini, S., Pierini, D., Finoguenov, A., et al. 2009, ApJ, 703, 982 [Google Scholar]
 GómezÁlvarez, P., Mediavilla, E., Muñoz, J. A., et al. 2006, ApJ, 645, L5 [Google Scholar]
 Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14 [Google Scholar]
 Green, P. J. 2006, ApJ, 644, 733 [Google Scholar]
 Hainline, L. J., Morgan, C. W., MacLeod, C. L., et al. 2013, ApJ, 774, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Hartley, P., Jackson, N., Badole, S., et al. 2021, MNRAS, 508, 4625 [Google Scholar]
 Henden, N. A., Puchwein, E., & Sijacki, D. 2020, MNRAS, 498, 2114 [Google Scholar]
 Hutsemékers, D., Sluse, D., Savić, D., & Richards, G. T. 2023, A&A, 672, A45 [Google Scholar]
 Inada, N., Oguri, M., Pindor, B., et al. 2003, Nature, 426, 810 [Google Scholar]
 Jackson, N. 2020, ApJ, 900, L15 [NASA ADS] [CrossRef] [Google Scholar]
 James, F. 2006, Statistical Methods in Experimental Physics, 2nd edn. (World Scientific Publishing Co., Pte. Ltd.) [Google Scholar]
 JiménezVicente, J., & Mediavilla, E. 2022, ApJ, 941, 80 [Google Scholar]
 JiménezVicente, J., Mediavilla, E., Muñoz, J. A., & Kochanek, C. S. 2012, ApJ, 751, 106 [Google Scholar]
 JiménezVicente, J., Mediavilla, E., Kochanek, C. S., et al. 2014, ApJ, 783, 47 [Google Scholar]
 JiménezVicente, J., Mediavilla, E., Kochanek, C. S., & Muñoz, J. A. 2015, ApJ, 799, 149 [Google Scholar]
 Kochanek, C. S. 2004, ApJ, 605, 58 [Google Scholar]
 Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astron. Lett., 44, 8 [Google Scholar]
 Lamer, G., Schwope, A., Wisotzki, L., & Christensen, L. 2006, A&A, 454, 493 [Google Scholar]
 Liesenborgs, J., Perera, D., & Williams, L. L. R. 2024, MNRAS, 529, 1222 [Google Scholar]
 MacLeod, C. L., Morgan, C. W., Mosquera, A., et al. 2015, ApJ, 806, 258 [Google Scholar]
 McKean, J. P., Luichies, R., Drabent, A., et al. 2021, MNRAS, 505, L36 [Google Scholar]
 Mediavilla, E., Muñoz, J. A., Falco, E., et al. 2009, ApJ, 706, 1451 [Google Scholar]
 Morgan, C. W., Hyer, G. E., Bonvin, V., et al. 2018, ApJ, 869, 106 [Google Scholar]
 Mortonson, M. J., Schechter, P. L., & Wambsganss, J. 2005, ApJ, 628, 594 [Google Scholar]
 Mosquera, A. M., & Kochanek, C. S. 2011, ApJ, 738, 96 [Google Scholar]
 Motta, V., Mediavilla, E., Falco, E., & Muñoz, J. A. 2012, ApJ, 755, 82 [Google Scholar]
 Muñoz, J. A., VivesArias, H., Mosquera, A. M., et al. 2016, ApJ, 817, 155 [Google Scholar]
 Muñoz, J. A., Kochanek, C. S., Fohlmeister, J., et al. 2022, ApJ, 937, 34 [Google Scholar]
 Oguri, M. 2010, PASJ, 62, 1017 [NASA ADS] [Google Scholar]
 Perera, D., Williams, L. L. R., Liesenborgs, J., Ghosh, A., & Saha, P. 2024, MNRAS, 527, 2639 [Google Scholar]
 Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648 [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [Google Scholar]
 Poindexter, S., & Kochanek, C. S. 2010, ApJ, 712, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Popović, L. Č., Afanasiev, V. L., Moiseev, A., et al. 2020, A&A, 634, A27 [EDP Sciences] [Google Scholar]
 Richards, G. T., Keeton, C. R., Pindor, B., et al. 2004, ApJ, 610, 679 [Google Scholar]
 Ross, N. R., Assef, R. J., Kochanek, C. S., Falco, E., & Poindexter, S. D. 2009, ApJ, 702, 472 [Google Scholar]
 Vernardos, G., & Tsagkatakis, G. 2019, MNRAS, 486, 1944 [Google Scholar]
All Tables
Inferred stellar mass fraction and convergence in stars at the locations of the four quasar images.
All Figures
Fig. 1. Fits to the rband light curves B, C, and D from Muñoz et al. (2022) using fifthorder splines (solid black curves). Since gaps longer than 100 days were not considered in the analysis, the splines in these regions are not displayed. The days on the xaxis correspond to the observation dates. 

In the text 
Fig. 2. Microlensing differences between the six possible combinations of image pairs. The solid black lines are the mean of the micromagnification, and the dashed lines correspond to the standard deviation. All light curves are shifted to match the time reference of the leading image, C. 

In the text 
Fig. 3. Observed histograms of the six possible microlensing differences obtained after averaging the 10^{6} realisations from the Monte Carlo sampling. The error bars are the standard deviation of these realisations. The solid black lines represent the best model difference histograms, and the thin dotted lines are other model histograms used for the parameter inference. 

In the text 
Fig. 4. Joint Δχ^{2} = χ^{2} − χ_{min}^{2} contours for all pairs of parameters (R_{1/2}, α_{A}, α_{B}, α_{C}, and α_{D}). In the 2D likelihoods, the contours are drawn at the 68% (Δχ^{2} = 2.30) and 95% (Δχ^{2} = 6.18) confidence limits for two parameters. In the 1D distributions, the shaded regions are the 68% (Δχ^{2} = 1.00) confidence regions for one parameter. The maximum likelihood model is marked with a cross and a dashed line for the 2D and 1D plots, respectively. 

In the text 
Fig. 5. Magnification maps normalised to the mean macromagnification, ⟨μ⟩, for images A and B with the 500 pairs of tracks of the best fitting model. These tracks are randomly oriented, and their microlensing difference histograms were averaged to build the model histogram. The radius of the red circle at the bottom of each map corresponds to the inferred halflight radius of the source (i.e. 0.7 R_{E} or equivalently 6.4 lightdays). 

In the text 
Fig. 6. Comparison with previous quasar accretion disk size estimates for SDSS J1004+4112. The results are reported as the rband (λ_{obs} = 6580 Å, λ_{rest} = 2407 Å) halflight radius in lightdays for a mean stellar mass of 0.3 M_{⊙} to homogenise the comparison. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.