© The Authors 2023

1. Introduction

The quasar J1004+4112 (SDSS J100434.80+411239.2 for component A) at redshift z = 1.734 is a gravitationally lensed quasar with four images separated by angular distances between 3.7″ and 14.6″, with the lens being a galaxy cluster at z = 0.68 (Inada et al. 2003; Oguri et al. 2004). Although the four components have the same redshift and display similar spectra, there are significant differences between the broad emission-line profiles observed in the different images. In particular, image A exhibits a strong, time-dependent, blue emission-line wing enhancement in the high-ionization lines, interpreted with a varying microlensing magnification of the broad emission-line region (BLR), rapid spectral changes intrinsic to the quasar and seen with a time delay between the different images being excluded (Richards et al. 2004b). However, the absence of apparent continuum microlensing in image A and the recurrence of some line profile deformations have been considered as a challenge for the microlensing interpretation (Gómez-Álvarez et al. 2006). Lamer et al. (2006) proposed a mixed scenario involving a stationary microlensing effect, with the time-dependent changes being due to intrinsic variability rather than to a moving caustic, under the assumption that image B was not microlensed itself. Based on X-ray spectra, these authors also ruled out the possibility that the line profile differences observed between the four images are due to different lines of sight through quasar outflows that lead to differential absorption (Green 2006). In parallel, theoretical simulations demonstrated that the microlensing of a biconical BLR can explain the observations (Abajas et al. 2007). More recent observations (Motta et al. 2012; Popović et al. 2020) showed that the blue emission-line wing enhancement in component A, although variable in strength, is a long-term effect. With the evidence of chromatic microlensing of the continuum (Motta et al. 2012), these observations provided additional support to the microlensing interpretation.

In the present paper, we revisit the microlensing effect in J1004+4112 using spectra obtained over 15 years including previously unpublished ones. Our analysis takes advantage of the fact that at one epoch, image B was not affected by microlensing and thus may be used as a reference to unambiguously characterize the microlensing effect at work in image A. It also benefits from recent radio measurements that provide the macro-magnification ratios of the different components (Hartley et al. 2021). The data are presented in Sect. 2. The microlensing-induced deformations of the C IV broad emission line profile are discussed and quantified in Sect. 3. In order to characterize the BLR geometry and kinematics, and estimate its size, as was recently done for other lensed quasars (Hutsemékers et al. 2019; Hutsemékers & Sluse 2021), we compare the observed C IV line profile deformations to microlensing simulations based on simple, representative BLR models magnified by a caustic network. The method is described in Sect. 4 and the results are given in Sect. 5. Our conclusions are given in the final section.

2. Data

For our analysis, we used a series of 12 spectra obtained from November 2003 to May 2006 with the ARC 3.5 m telescope at Apache Point Observatory. The spectra of images A and B were secured using the same instrument setup and data reduction for each epoch. Details on this data set are given in Richards et al. (2004b), with an analysis of the 2003 subsample.

We also considered the higher quality spectra of images A, B, C, and D, obtained with the Keck I telescope in May 2003. These spectra are described in Inada et al. (2003), Oguri et al. (2004), and Richards et al. (2004b).

Two series of good quality spectra of the four images obtained in 2008 and 2018 with the 6 m telescope of the Special Astrophysical Observatory (SAO) were also added to our sample. These spectra are described in Popović et al. (2020) and publicly available at the Strasbourg astronomical Data Center (CDS).

Table 1 lists the main characteristics of the considered spectra; Δλ is the spectral resolution. In total, there are 15 spectra unevenly distributed between 2003 and 2018. In the following, we focus on the C IV broad emission line which is strong, unblended with other emission lines over most of its profile, and present in all spectra.

3. Broad emission-line microlensing in J1004+4112

Figure 1 shows the continuum-subtracted C IV emission line profiles observed in 2003 and 2004 in images A and B. The continuum was measured in two windows on each side of the line profile, [3925, 3980] and [4660, 4715] Å, which correspond to [ − 22, −18] and [ + 30, +34] 10³ km s⁻¹. It was then interpolated by a straight line under the line profile. The line profiles are superimposed according to different normalizations, either with respect to the line peak or to the blue wing. The flux enhancement in the blue wing of component A in May 2003, followed by its decrease later in 2003, and its recovery in March 2004 (Richards et al. 2004a,b; Gómez-Álvarez et al. 2006) are best seen when normalizing the line profiles to the peak intensity. The enhancement and its variability were interpreted by a moving microlensing caustic (Richards et al. 2004b), although the recurrence of the same feature has cast doubt on this interpretation (Gómez-Álvarez et al. 2006). Normalizing the line profiles to the blue wing tells another story: while the net enhancement of the blue wing with respect to the blue wing seen in image B is still present, the profile variations appear as global changes that include the line core. This led Lamer et al. (2006) to propose a scenario in which the blue wing is affected by a long-term, stable microlensing effect in image A, and where the observed changes are due to intrinsic variability on a timescale of a few months. This scenario is supported by the fact that similar variations are seen in the image B spectrum (Fig. 1), as predicted by Lamer et al. (2006), although it does not fully explain the extreme behavior of the May 2003 spectrum. These results illustrate the complex intertwining of microlensing and intrinsic variations. In the remainder of this section, we try to disentangle these two effects.

Fig. 1. Continuum-subtracted C IV emission line profiles observed in 2003 and 2004. Top: profiles in image A normalized to the peak. Middle: profiles in image A normalized to the blue wing. Bottom: profiles in image B normalized to the blue wing.

To characterize and quantify the line profile deformations induced by the microlensing effect, we considered the magnification profile μ(v), which is the macro-magnification-corrected ratio of the continuum-subtracted emission line flux densities observed in two different images. We also considered three indices: 1) μ^BLR, integrated over the line profile; it essentially quantifies the total magnification of the line, 2) WCI, the wings-core index integrated over the μ(v) profile; when it is above or below a value of 1, it indicates whether the whole emission line is, on average, more or less affected by microlensing than its center, and 3) RBI, the red-blue index integrated over the μ(v) profile; it takes non-null values when the effect of microlensing on the blue and red parts of the line is asymmetric. A fourth index, μ^cont, measures the microlensing magnification of the continuum underlying the emission line. All these quantities were computed using the spectrum of one image affected by microlensing and the spectrum of a second image not or much less affected. When the spectra of both images are simultaneously recorded, μ(v) and the indices are independent of quasar intrinsic variations that occur on timescales longer that the time delay between the two images. The indices μ^cont and μ^BLR are corrected for the image macro-magnification ratio; RBI and WCI are independent of this ratio. An exact definition of these quantities can be found in Braibant et al. (2017) or Hutsemékers & Sluse (2021).

Since microlensing does not affect the radio emission, which originates from a more extended source than the optical emission, in the following, we use the macro-magnification ratio, M_A/M_B = 1.60 ± 0.15, computed from the radio fluxes given in Hartley et al. (2021). Figure 2 shows examples of μ(v) profiles computed at four epochs from images A and B. As a flux ratio, μ(v) can be extremely noisy in the wings of the emission lines where the flux density reaches zero, so that it is necessary to cut the faintest parts of the line wings. Thus, we only considered the parts of the line profiles whose flux density is above l_cut × F_peak, where F_peak is the maximum flux in the line profile and l_cut is fixed around 0.1. In practice, this corresponds to cut the profiles between −10⁴ and +10⁴ km s⁻¹ at all epochs. To increase the signal-to-noise ratio, we binned μ(v) into 20 spectral elements, which also corresponds to the spectral resolution of the line profiles produced by the microlensing simulations (see Sect. 4).

Fig. 2. Examples of μ(v) magnification profiles of C IV at four epochs (in red) computed from spectra of images A and B that were simultaneously recorded. The μ(v) profiles are binned into 20 spectral elements. The superimposed line profiles from image A (in blue) and B (in black) are continuum-subtracted, corrected for the macro-magnification ratio, and arbitrarily rescaled for convenience. The thicker line indicates the part used in the computation of μ(v). The error of μ(v) is shown in green. The observation date is given as yymmdd. The zero-velocity corresponds to the C IVλ1549 wavelength at the redshift z = 1.734.

A proper interpretation of the microlensing effect in one image requires a reference image that is not affected by microlensing, as well as a short time delay between the two images, so as to rule out the possibility that observed spectral differences could be attributed to a delayed, short timescale intrinsic variability. The time delay between images B and A is 40.6 ± 1.8 days (Fohlmeister et al. 2008), which is small with respect to the timescale of the intrinsic variations shown in Fohlmeister et al. (2008) and Muñoz et al. (2022). We may therefore assume that there is only a negligible delay between images A and B when a significant intrinsic variation occurs; however, we ought to keep in mind that for some rare events, this might not be the case. This assumption is supported by simulations that show that (statistically) intrinsic variability is unlikely to mimic microlensing-induced line deformations if the time delay is shorter than about 50 days (Sluse et al. 2012).

In 2018 (and only in 2018), the C IV continuum and line profile of images B and C are nearly ideally superimposed with a constant factor, indicating the absence of microlensing in both images B and C (Fig. 3). The constant factor M_B/M_C = 1.40 ± 0.14 is in agreement, within the uncertainties, with the macro-magnification ratio derived from radio observations, M_B/M_C = 1.10 ± 0.13 (Hartley et al. 2021), keeping in mind that the optical ratio could be slightly off because the A-B and C-D spectra were not obtained simultaneously, making the relative flux calibration of B and C more uncertain (Popović et al. 2020). The B/C magnification profile μ(v) is flat (Fig. 3). The microlensing indices μ^cont, μ^BLR, and WCI measured between images B and C are around one, and RBI around zero, consistently indicating the absence of any differential microlensing affecting the C IV broad emission line in images B and C. Without microlensing, the line profiles appear quasi-symmetric.

Fig. 3. μ(v) magnification profile of C IV (in red) computed from the spectra of images B and C obtained in May 2003 and February 2018, and binned into 20 spectral elements, as in Fig. 2. The line profiles from image B (in blue) and C (in black) are continuum-subtracted, corrected for the macro-magnification ratio and arbitrarily rescaled for convenience. The error of μ(v) is shown in green.

By comparing the spectra of images A and B at that epoch using B as the non-microlensed reference, we can thus unambiguously characterize the microlensing effect observed in image A. We find that the continuum of image A is significantly demagnified with μ^cont = 0.59 ± 0.10, differential dust extinction being negligible (Motta et al. 2012). The C IV emission line profile is strongly distorted (Fig. 2): the blue part of the profile is magnified by up to a factor of two, the red part demagnified by almost the same factor, reaching the continuum demagnification μ^cont at the highest velocities, while the line core is essentially unaffected by the microlensing effect.

As seen in Figs. 2 and 4, the μ(v) profiles derived from images A and B retain roughly the same shape at the different epochs, with significant variability in the blue and red parts of the profile. The fact that μ(v = 0) remains stable at all epochs confirms that these variations are not due to delayed, short-timescale intrinsic variations that would otherwise mimic the microlensing magnification pattern. Since μ(v) is independent of long-term intrinsic profile variations, the μ(v) profile variability is thus likely due to microlensing variations in either image A or image B, or both. Time series of the four indices are shown in Fig. 5. Consistently, at all epochs, WCI > 1 indicates higher magnification of the line wings than the peak, and RBI < 0 indicates higher magnification of the blue part of the line than the red part. μ^BLR remains nearly constant, namely, around one, indicating no significant magnification of the line profile as a whole, the demagnification of the red part of the line profile counterbalancing the magnification of the blue part. μ^cont shows significant variation, reaching μ^cont ≃ 0.6 − 0.7 after 2006. As for the μ(v) profile, the variability of the indices is likely due to microlensing variability in component A, B, or both.

Fig. 4. Time series of μ(v) profiles of C IV computed from the spectra of images A and B obtained at the same epoch. Four epochs are emphasized: May 2003 (blue), May 2006 (cyan), October 2008 (magenta), and February 2018 (red); the other μ(v) profiles computed from spectra obtained in 2003 and 2004 are shown in green and black, respectively. The half error bars of μ(v) are drawn from zero for clarity.

Fig. 5. Time series of the four indices computed from the spectra of images A and B obtained at the same epoch. The two panels have different scales for the abscissae. At all epochs, WCI > 1 indicates higher magnification of the line wings than the peak, and RBI < 0 indicates a higher magnification of the blue part of the line than the red part.

In Fig. 6, we show μ(v) computed from the spectra of image A at all epochs with respect to the spectrum of image B obtained in 2018 (hereafter, B2018), which is unaffected by microlensing. We also compute μ(v) from the spectra of image B at all epochs with respect to its 2018 spectrum. As expected, variations of the B/B2018 μ(v) profile are observed. They are due to intrinsic variations between 2003 and 2018, microlensing in image B, and absolute flux calibration errors. We then normalize the μ(v) profiles by the factor f_n = F_BEL(B2018)/F_BEL(B), where F_BEL is the flux integrated over the continuum-subtracted C IV line profile. This normalization corrects for calibration errors as well as global intrinsic and microlensing variations, but not for line profile deformations. After normalization, time-dependent μ(v) profile deformations are still present in image B, due to intrinsic or microlensing velocity-dependent variations in image B. Similarly, time-dependent μ(v) profile deformations are present in image A with respect to B2018; in this case, they are due to intrinsic or velocity-dependent microlensing variations in image A.

Fig. 6. Time series of μ(v) profiles of C IV computed from the spectra of images A and B with respect to the 2018 spectrum of image B, unaffected by microlensing (top). Same profiles are shown below, but normalized by the factor, fn, defined in the text. The color-coding of the profiles and the error bars are the same as in Fig. 4. The A/B2018 μ(v) profiles are not affected by microlensing in B, while the B/B2018 μ(v) profiles are not affected by microlensing in A. The normalized A/B2018 μ(v) profiles have similar shapes at the different epochs.

If only of intrinsic origin, these variations would disappear from the A/B μ(v) profiles derived from images A and B observed at the same epoch (Fig. 4), resulting in a single A/B μ(v) profile for all epochs. Since this is not observed, microlensing variations are thus present in both images A and B, affecting the A/B2108 and B/B2018 μ(v) profiles separately and contaminating the A/B μ(v) profiles seen in Fig. 4. The normalized A/B2018 μ(v) profiles, which are not affected by microlensing in B, show less scatter than the A/B μ(v) profiles as well as more similar shapes through the different epochs, suggesting that they are closer to the true microlensing effect that occurs in image A.

In summary, we can see that although they are variable, the microlensing-induced distortions of the C IV line profile seen in image A appear remarkably similar over a period of 15 years. The distortions are characterized by a magnification of the blue part of the profile and a demagnification of the red part, with only a small demagnification of the line core. The magnification ⟨μ′⟩ = ⟨μ × f_n⟩ averaged over the different epochs (Fig. 6) approximately decreases as log₁₀⟨μ′⟩ ≃ −0.4 × (v/10 000 km s⁻¹). The stronger demagnification observed in the velocity range 1000–5000 km s⁻¹ in the two spectra obtained after 2008 (in magenta and red) may indicate a long-term change (see also RBI in Fig. 5). Microlensing is also at work in image B, albeit more weakly. In the following, we compare the μ(v) profile and the indices to simulations, focusing on the 2018 data for which the microlensing effect can be unambiguously characterized.

4. Microlensing simulations

4.1. Methodology and model parameters

We computed the effect of gravitational microlensing on the broad emission line profiles by convolving in the source plane the emission from representative BLR models with microlensing magnification maps. For the BLR models, we considered a rotating Keplerian disk (KD), as well as biconical polar (PW) and equatorial (EW) radially accelerated winds characterized by inclinations with respect to the line of sight of i = 22°, 34°, 44°, 62°. We considered 11 BLR inner radius values: r_in = 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.25, 0.35, and 0.5r_E, where r_E is the microlensing Einstein radius in the source plane. The outer radius is fixed to r_out = 10 r_in. The BLR models are assumed to have an emissivity ϵ = ϵ₀ (r_in/r)^q that either sharply decreases with increasing radius, q = 3, or more slowly, q = 1.5. Twenty BLR monochromatic images forming a data cube and corresponding to twenty spectral bins in the line profile are produced, using the radiative transfer code STOKES (Goosmann & Gaskell 2007; Marin et al. 2012; Goosmann et al. 2014). Each model also contains a continuum-emitting uniform disk seen under the same inclination as the BLR, and with an outer radius fixed at 11 different values: r_s = 0.02, 0.03, 0.05, 0.07, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, and 0.5r_E. Only models with r_in ≥ r_s are considered. For a full description of the models and additional details, we refer to Braibant et al. (2017).

Modeling the effect of microlensing on the BLR is achieved using a magnification map specific to the image A of J1004+4112. The magnification map is computed using the microlens ray-tracing code (Wambsganss 1999), considering a total convergence of κ= 0.763 with a shear γ= 0.3. The fraction of matter in compact objects is arbitrarily taken to be lower than usual because the system is lensed by a cluster: κ_⋆/κ = 3%. These values are from Guerras et al. (2017) based on the model of Oguri (2010) and they are in agreement with more recent modeling by Forés-Toribio et al. (2022). Since κ_⋆/κ is poorly known, we tested the robustness of the results using κ_⋆/κ = 1% and κ_⋆/κ = 7%, the latter value being more typical of the fraction of compact objects in lens galaxies (Mediavilla et al. 2009). The map extends over a 200 r_E × 200 r_E area of the source plane and is sampled by 20 000 × 20 000 pixels. To investigate the impact of preferential alignment between the symmetry axis of the BLR models and the caustic network, the maps are rotated clockwise by θ = 15°, 30°, 45°, 60°, 75°, and 90° with respect to the BLR model axis. Here, θ = 0° corresponds to the caustic elongation and shear direction perpendicular to the BLR model axis. After rotation, only the central 10 000 × 10 000 pixel part of the map is used.

Distorted line profiles are obtained from the convolution of the monochromatic images of the BLR by the magnification maps and computed for each position of the BLR on the magnification maps, generating ∼10⁸ simulated profiles per map and BLR model. In previous works, we only considered the indices μ^cont, μ^BLR, RBI, and WCI for the comparison with the observations. In the present work, we modified the code to accelerate the convolution and save disk space so that we can also use the full μ(v) magnification profiles with 20 spectral elements. For a proper comparison, the wings of the simulated line profiles are truncated as the observed ones, so that only the parts of the profiles with a flux density larger than 0.1 × F_peak are considered. After truncation, the simulated line profiles, and thus the data cubes, are finally rescaled to match the observed μ(v) profile in velocity. Simulations show that the final probabilities of the different BLR models are robust with respect to small changes of the adopted flux threshold.

4.2. Probability of the different BLR models

In order to identify the BLR models that provide the best fit to the observations, we computed the relative posterior probability that a given model (G, i) (where G = KD, PW, or EW, and i = 22°, 34°, 44°, or 62°) can reproduce the four observables μ^cont, μ^BLR, RBI, and WCI. The method is described in detail in Hutsemékers et al. (2019). In practice, we computed the likelihood of the observables for each set of parameters characterizing the simulations. We then marginalized the likelihood over r_s, r_in, and q, as well as over the microlensing parameters. Since the different BLR models share the same parameters and associated priors, we can quantify their relative efficiency to reproduce the data by comparing their likelihoods, that is, by normalizing the marginalized likelihood by the sum of the likelihoods associated to each model, G, for each inclination, i. This procedure yields the relative probability of the different models.

The probabilities of the different models are similarly computed using μ^cont and the 20 spectral elements of μ(v) instead of the three indices used to characterize the line profile magnification. However, given the higher number of degrees of freedom and the non-continuous nature of the model parameters, the previously used likelihood estimator L_G = exp (−0.5  χ²), which is very sensitive to small, rare χ² values, may provide unreliable results. In particular, observed μ(v) profiles comparable within the uncertainties can result in very different probabilities of the models. To circumvent this problem, we considered the likelihood estimator used by Kochanek (2004):

$$\begin{array}{c}\hfill L_K=\mathrm{\Gamma }[\frac{n_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}-2}{2},\frac{{\chi }^2}{2f_0^2}],\end{array}$$(1)

where n_d.o.f. = 19 is the number of degrees of freedom, and f₀ a free parameter. This likelihood estimator is essentially flat for small χ²/n_d.o.f. values, and behaves like the L_G estimator for higher χ² values. With the L_K estimator and f₀ fixed to one, the resulting model probabilities were found to be stable with respect to small changes of the observed μ(v) profile.

5. Results

For the comparison with the microlensing simulations, we used the μ(v) magnification profile measured from the spectra of images A and B obtained in 2018, with image B being unaffected by microlensing at that epoch. The indices measured from this data set are μ^cont = 0.59 ± 0.10, μ^BLR = 1.01 ± 0.08, WCI = 1.24 ± 0.11, and RBI = −0.50 ± 0.03.

In Fig. 7, we show examples of simulated μ(v) profile for the KD and EW models (the most likely BLR models; see Table 2). As can be seen, the observed μ(v) profile is reproduced by a number of simulated profiles, without the need for more sophisticated BLR models, in particular, models involving two BLRs of different sizes. Figure 8 shows typical locations of the BLR on the magnification map, at which the observed line profile deformations and the continuum demagnification are reproduced. As anticipated, one of the high-velocity parts of the BLR coincides with high-magnification caustics, with the other part and the continuum source lying in a demagnification region. The microlensing effect appears dominated by the distance of the different parts of the BLR with respect to the caustics rather than by the velocity-dependent BLR size. Indeed, the effective radius of the BLR is smaller at high velocities for the KD model, while it is larger for the EW model (at least for inclinations lower than about 60°; Braibant et al. 2017). If the microlensing magnification were dominated by the source size, we would expect higher magnification in the line wings for the KD models and in the line core for the EW models. On the contrary, both models reproduce the observed μ(v) profile that shows high (de)magnification in the wings only. Moreover, the change of the effective BLR radius with the velocity remains smaller than a factor three (Braibant et al. 2017). Since the magnification varies as $1/\sqrt{r}$ when close to a caustic, with r denoting the source radius (Schneider et al. 1992), the change in radius would induce at most a factor of$\sqrt{3}$ in magnification, which is lower than the magnification ratio observed between the different parts of the line profile.

Fig. 7. Examples of 20 simulated μ(v) profiles (in color) that fit, with χ2/nd.o.f. ≤ 1.2, the μ(v) magnification profile measured in 2018 for the C IV emission line (in black). The simulated profiles are computed with i = 34°, q = 3, rin = 0.15 rE, and rs = 0.03 rE. The left and right panels show profiles from the KD and EW models computed with the magnification map oriented at θ = 60° and θ = 30°, respectively.

Similar microlensing effects are obtained from close positions on the maps, which may extend over regions larger than one r_E (Fig. 8). The timescale to cross one r_E is about 10 years (Mosquera et al. 2011), so that the observation of similar line-profile deformations over 15 years is plausible. Since crossing caustics would typically produce blue-red sequences of magnification in the line profiles (Abajas et al. 2007; Braibant et al. 2017), the similarity of the magnification profiles indicates the absence of caustic crossing during this time interval. Interestingly, tracks that do not cross caustics are more likely to occur in lens systems with stretched caustic networks like J1004+4112 (Fig. 8) than in systems with more randomly oriented caustics such as Q2237+0305 (see Hutsemékers & Sluse 2021).

Fig. 8. Examples of positions on the magnification map where the computed line profile deformations fit the observed one with χ2/nd.o.f. ≤ 1.2 and where the computed indices match the observed ones within the error bars. The exact positions are at the center of the circles, the radius of which represents the BLR half-light radius r1/2 = 1.64 rin = 0.25 rE. The left panel shows the map convolved by the part of the BLR at the origin of the blue part of the line profile, while the right panel shows the map convolved by the part of the BLR at the origin of the red part of the line profile. The unconvolved magnification map is superimposed and seen as narrow caustics. The considered model is KD with i = 34°, q = 3, rin = 0.15 rE, and rs = 0.03 rE. The magnification map, computed with κ⋆/κ = 3%, is oriented at θ = 60°. The size of the illustrated part of the map is 8 rE × 8 rE. Similar patterns are observed throughout the full map, as well as for the EW model with θ ≤ 30°.

5.1. Geometry and kinematics of the BLR

The probabilities of the different models of the C IV BLR are given in Table 2. We found no significant differences using the magnification maps computed with the different κ_⋆/κ values. The results obtained using the full μ(v) profile are in good agreement with those derived from the integrated indices. The model selection is mostly driven by the large RBI value, as seen in WCI − RBI diagrams (Braibant et al. 2017; Hutsemékers et al. 2019; Hutsemékers & Sluse 2021). Since the probablities were found to strongly depend on the orientation of the magnification maps with respect to the BLR axis, we computed them individually for the map orientations θ ≤ 30° and θ ≥ 60°. In the first case, EW models are strongly favored, while in the second case, KD models are the only possibility. In all cases, the PW models, within the framework of which high RBI values are more difficult to reproduce (Braibant et al. 2017; Hutsemékers et al. 2019; Hutsemékers & Sluse 2021), are much less likely and can thus be rejected. The dichotomy in the probabilities of the KD and EW models with respect to the magnification map orientation is due to the fact that the caustic network is especially stretched in one direction (Fig. 8); in addition, in the plane of the sky, the high-velocity regions of the BLR are perpendicular to the BLR axis in the KD models, whereas they are parallel in the EW models (Braibant et al. 2017). Unfortunately, no constraint on the radio jet orientation can be derived from the Very Large Array maps shown in Hartley et al. (2021). Using the jet position angle as a proxy of the BLR axis orientation, it would have been possible to disentangle the KD and EW models.

5.2. Size of the BLR

By marginalizing over all parameters but r_in, we can estimate the most likely BLR radius. Since r_in does not properly represent the size of the BLR, which also depends on the light distribution, we computed the half-light radii for the different models as in Hutsemékers & Sluse (2021). The resulting probabilities of the half-light radii, r_1/2, are shown in Fig. 9. The probability distributions depend on κ_⋆/κ, but the most likely value appears to be independent of it. This value is equal to $r_{1/2}=0.3_{-0.18}^{+0.22}{r}_{\mathrm{E}}$, with the uncertainties corresponding to a 68% confidence interval estimated for κ_⋆/κ = 3%¹. Similar values are obtained when marginalizing over all values of κ_⋆/κ. To compute the Einstein radius, we adopted a flat ΛCDM cosmology with H₀ = 68 km s⁻¹ Mpc⁻¹ and Ω_m = 0.31. The source and lens redshifts are z_s = 1.734 and z_L = 0.68, respectively (Inada et al. 2003), so that the Einstein radius is $r_{\mathrm{E}}=9.17\sqrt{M/0.3M_{\odot }}$ light-days in the source plane. For an average microlens mass of 0.3 M_⊙, we obtained: $r_{1/2}=2.8_{-1.7}^{+2.0}$ light-days. This value is unchanged when considering the map orientations θ ≤ 30° or θ ≥ 60° separately or when considering specific model inclinations. Similarly, we can estimate the half-light radius of the continuum source using $r_{1/2}=r_{\mathrm{s}}/\sqrt{2}$ for a uniform disk. The probability distribution is shown in Fig. 10, only providing an upper limit, r_1/2 < 1.8 light-days (95% confidence level) at λ_rest = 1550 Å, consistently smaller than the BLR size. This upper limit translates to r_1/2 < 2.9 light-days at 2400 Å, assuming r_1/2 ∝ λ^p with p = 1.1 (Motta et al. 2012). It is in agreement with the continuum source size estimated by Fohlmeister et al. (2008) using the light curve fitting method, $r_{1/2}=0.8_{-0.4}^{+0.8}$ light-days at λ_rest = 2300 Å, but smaller than the size computed by Fian et al. (2016) based on the single epoch method, $r_{1/2}=4.2_{-2.2}^{+3.2}$–$8.7_{-5.5}^{+18.5}$ light-days at λ_rest = 2407 Å.

Fig. 9. Relative probabilities (in percent) of the BLR half-light radius r1/2 in rE and light-day units, considering the magnification maps computed with κ⋆/κ = 1, 3, and 7% (in blue). The probability distribution obtained by marginalizing over all κ⋆/κ is shown in red.

Fig. 10. Same as Fig. 9 but for the continuum half-light radius.

Rakshit et al. (2020) reported the luminosity λL_λ(1350 Å) for image A in January 2016. Divided by the macro-model magnification factor 29.7 (Oguri 2010), it amounts to λL_λ(1350 Å) = 1.9 × 10⁴⁴ ergs s⁻¹. Popović et al. (2020) measured this luminosity in February 2018, which, corrected for the macro-model magnification factor 29.7, amounts to λL_λ(1350 Å) = 3.2 × 10⁴⁴ ergs s⁻¹. Since the object does not seem to be brighter in 2018 than in 2016 (Muñoz et al. 2022), hereafter, we use the average value λL_λ(1350 Å) = 3.6 ± 0.9 × 10⁴⁴ ergs s⁻¹, corrected for the microlensing demagnification μ^cont ≃ 0.7 that is likely to occur for the 2016–2018 data (Fig. 5). With this luminosity, we derived the expected C IV BLR radius R_BLR from the radius-luminosity (R − L) relations of Kaspi et al. (2021) based on reverberation mapping: $R_{\mathrm{BLR}}={15}_{-7}^{+15}$ light days. The reverberation mapping BLR radius is the rest-frame time lag measured as the centroid of the cross-correlation function between the continuum and line light curves. It corresponds (in theory) to the luminosity-weighted radius of the BLR (Koratkar & Gaskell 1991) and can thus be compared to the half-light radius derived from microlensing. Our estimate of the BLR size, $r_{1/2}=2.8_{-1.7}^{+2.0}$ light-days, is smaller by a factor of 5. However, the errors are large: a radius of 2.8 light-days represents only a 1.7σ departure from the R − L relation, which shows several outliers with comparable differences (Fig. 11). Given that the range of BLR sizes involved in the R − L relation covers four orders of magnitude, our BLR radius appears to be in reasonable agreement, as it is only marginally smaller than expected. The BLR radius estimated for Q2237+0305 (Hutsemékers & Sluse 2021) is in excellent agreement with the R − L relations, indicating that the microlensing values are not biased². On the other hand, they depend on the microlens mass, which could differ from the average value of 0.3 M_⊙, especially in a system like J1004+4112, which is lensed by a cluster. For the microlensing BLR radius to precisely fit the R − L relation, the microlens mass should be changed from the average value of 0.3 M_⊙ to 7.5 M_⊙.

Fig. 11. Relation between the radius of the C IV BLR, that is, the rest-frame time lag from reverberation mapping and the continuum luminosity at 1350 Å, with different fits superimposed as continuous lines (from Kaspi et al. 2021). The BLR half-light radii measured from our microlensing analysis of the quasars J1004+4112 and Q2237+0305 are superimposed in red.

5.3. Black hole mass

Assuming virial motion, the mass of the black hole can be estimated based on:

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=f\frac{R_{\mathrm{BLR}}{V}^2}{G},\end{array}$$(2)

where V is the velocity full width at half-maximum (FWHM) of the broad emission line, f is the virial factor, and G is the gravitational constant. For a thin disk, Decarli et al. (2008) give³:

$$\begin{array}{c}\hfill f={(4{sin}^{2}i)}^{-1}.\end{array}$$(3)

For the broad C IV line observed in the 2018 spectra of image B (unaffected by microlensing), we measured FWHM = 4900 ± 300 km s⁻¹. Using i = 35 ± 10° and R_BLR = $r_{1/2}=2.8_{-1.7}^{+2.0}$ light-days, we finally obtained $M_{\mathrm{BH}}=1.0_{-0.7}^{+1.1}\times {10}^7$ M_⊙. To correct the bias that affects the C IV-based virial black hole masses, we used the prescriptions of Coatman et al. (2017). With a C IV blueshift of 350 ± 50 km s⁻¹ measured in the image B spectrum, the corrected black hole mass amounts to $M_{\mathrm{BH}}=1.7_{-1.2}^{+2.0}\times {10}^7$ M_⊙. As the BLR radius, this black hole mass is a factor of 4-5 smaller than the value expected from the black hole mass-luminosity (M_BH − L) relations of Kaspi et al. (2021), $M_{\mathrm{BH}}=7.5_{-4.2}^{+11.2}\times {10}^7$ M_⊙. It is nevertheless in agreement within the uncertainties, due to the large dispersion of the M_BH − L  correlation (Fig. 12). The black hole mass of Q2237+0305 derived in Hutsemékers & Sluse (2021) is in better agreement with the M_BH − L relation.

Fig. 12. Relation between the black hole mass and the continuum luminosity at 1350 Å, with different fits superimposed as continuous lines (from Kaspi et al. 2021). The black hole masses measured from our microlensing analysis of the quasars J1004+4112 and Q2237+0305 are superimposed in red.

6. Conclusions

Based on 15 spectra obtained from 2003 to 2018, we revisited the microlensing effect that distorts the C IV broad emission line profile in the lensed quasar J1004+4112. We took advantage of recent measurements of the image macro-magnification ratios and of the fact that, at one epoch, image B was not microlensed, thus constituting a reference spectrum that allowed us to unambiguously characterize the microlensing effect observed in image A. After disentangling the microlensing in images A and B, we showed that the microlensing-induced line profile distortions in image A, although variable, are remarkably similar over a period of 15 years. They are characterized by a strong magnification of the blue part of the line profile, a strong demagnification of the red part of the line profile, and a small to negligible demagnification of the line core. Microlensing also affects the C IV line profile in image B, but to a lesser extent.

We then compared the microlensing effect observed in image A, characterized by the magnification profile of the C IV emission line or by a set of four integrated indices, to simulated line-profile deformations that result from the selective magnification of BLR models by caustic networks. Our main conclusions are as follow.

• The observed μ(v) magnification profile of C IV in J1004+4112 can be reproduced with the simple BLR models we considered, without the need for more complex BLR features. The magnification pattern is dominated by the position of the BLR with respect to the caustic network – and not by the velocity-dependent size of the BLR. The observation of such a microlensing effect over 15 years has been deemed plausible.

• The favored models for the C IV BLR are either the Keplerian disk or the equatorial wind, depending on the orientation of the caustic map with respect to the BLR axis. The polar wind model is much less likely and can be rejected. In Q2237+0305, the most likely BLR model is that of a Keplerian disk (Hutsemékers & Sluse 2021).

• We measured the C IV BLR half-light radius as $r_{1/2}=2.8_{-1.7}^{+2.0}$ light-days. This result is marginally lower than the BLR radius expected from empirical R − L relations, while the C IV BLR half-light radius measured in Q2237+0305 (Hutsemékers & Sluse 2021) is in excellent agreement with these relations. More measurements are needed to compare the BLR sizes estimated from reverberation mapping with those derived from microlensing.

Finally, it is interesting to note that similar conclusions on the preferred BLR models can be reached using either the four indices or the full μ(v) magnification profile. While both can be used for single-epoch microlensing analyses, this result validates the use of indices, which can ultimately prove more convenient for interpreting time series of line-profile deformations.

Acknowledgments

D.H. and D.S acknowledge support from the Fonds de la Recherche Scientifique – FNRS (Belgium) under grants PDR T.0116.21 and No 4.4503.19.

References

All Tables

All Figures

	Fig. 1. Continuum-subtracted C IV emission line profiles observed in 2003 and 2004. Top: profiles in image A normalized to the peak. Middle: profiles in image A normalized to the blue wing. Bottom: profiles in image B normalized to the blue wing.
In the text

Fig. 2. Examples of μ(v) magnification profiles of C IV at four epochs (in red) computed from spectra of images A and B that were simultaneously recorded. The μ(v) profiles are binned into 20 spectral elements. The superimposed line profiles from image A (in blue) and B (in black) are continuum-subtracted, corrected for the macro-magnification ratio, and arbitrarily rescaled for convenience. The thicker line indicates the part used in the computation of μ(v). The error of μ(v) is shown in green. The observation date is given as yymmdd. The zero-velocity corresponds to the C IVλ1549 wavelength at the redshift z = 1.734.

In the text

	Fig. 3. μ(v) magnification profile of C IV (in red) computed from the spectra of images B and C obtained in May 2003 and February 2018, and binned into 20 spectral elements, as in Fig. 2. The line profiles from image B (in blue) and C (in black) are continuum-subtracted, corrected for the macro-magnification ratio and arbitrarily rescaled for convenience. The error of μ(v) is shown in green.
In the text

Fig. 4. Time series of μ(v) profiles of C IV computed from the spectra of images A and B obtained at the same epoch. Four epochs are emphasized: May 2003 (blue), May 2006 (cyan), October 2008 (magenta), and February 2018 (red); the other μ(v) profiles computed from spectra obtained in 2003 and 2004 are shown in green and black, respectively. The half error bars of μ(v) are drawn from zero for clarity.

In the text

	Fig. 5. Time series of the four indices computed from the spectra of images A and B obtained at the same epoch. The two panels have different scales for the abscissae. At all epochs, WCI > 1 indicates higher magnification of the line wings than the peak, and RBI < 0 indicates a higher magnification of the blue part of the line than the red part.
In the text

Fig. 6. Time series of μ(v) profiles of C IV computed from the spectra of images A and B with respect to the 2018 spectrum of image B, unaffected by microlensing (top). Same profiles are shown below, but normalized by the factor, fn, defined in the text. The color-coding of the profiles and the error bars are the same as in Fig. 4. The A/B2018 μ(v) profiles are not affected by microlensing in B, while the B/B2018 μ(v) profiles are not affected by microlensing in A. The normalized A/B2018 μ(v) profiles have similar shapes at the different epochs.

In the text

Fig. 7. Examples of 20 simulated μ(v) profiles (in color) that fit, with χ2/nd.o.f. ≤ 1.2, the μ(v) magnification profile measured in 2018 for the C IV emission line (in black). The simulated profiles are computed with i = 34°, q = 3, rin = 0.15 rE, and rs = 0.03 rE. The left and right panels show profiles from the KD and EW models computed with the magnification map oriented at θ = 60° and θ = 30°, respectively.

In the text

Fig. 8. Examples of positions on the magnification map where the computed line profile deformations fit the observed one with χ2/nd.o.f. ≤ 1.2 and where the computed indices match the observed ones within the error bars. The exact positions are at the center of the circles, the radius of which represents the BLR half-light radius r1/2 = 1.64 rin = 0.25 rE. The left panel shows the map convolved by the part of the BLR at the origin of the blue part of the line profile, while the right panel shows the map convolved by the part of the BLR at the origin of the red part of the line profile. The unconvolved magnification map is superimposed and seen as narrow caustics. The considered model is KD with i = 34°, q = 3, rin = 0.15 rE, and rs = 0.03 rE. The magnification map, computed with κ⋆/κ = 3%, is oriented at θ = 60°. The size of the illustrated part of the map is 8 rE × 8 rE. Similar patterns are observed throughout the full map, as well as for the EW model with θ ≤ 30°.

In the text

	Fig. 9. Relative probabilities (in percent) of the BLR half-light radius r1/2 in rE and light-day units, considering the magnification maps computed with κ⋆/κ = 1, 3, and 7% (in blue). The probability distribution obtained by marginalizing over all κ⋆/κ is shown in red.
In the text

	Fig. 10. Same as Fig. 9 but for the continuum half-light radius.
In the text

	Fig. 11. Relation between the radius of the C IV BLR, that is, the rest-frame time lag from reverberation mapping and the continuum luminosity at 1350 Å, with different fits superimposed as continuous lines (from Kaspi et al. 2021). The BLR half-light radii measured from our microlensing analysis of the quasars J1004+4112 and Q2237+0305 are superimposed in red.
In the text

	Fig. 12. Relation between the black hole mass and the continuum luminosity at 1350 Å, with different fits superimposed as continuous lines (from Kaspi et al. 2021). The black hole masses measured from our microlensing analysis of the quasars J1004+4112 and Q2237+0305 are superimposed in red.
In the text