© The Authors 2024

1. Introduction

Emission regions in distant active galactic nuclei (AGN) remain inaccessible with the current capabilities of existing telescopes owing to their small dimensions and the vast distances involved. Microlensing variability induced by stellar mass objects in the lens galaxy allows us to infer the sizes of these regions indirectly, thereby overcoming the limitations of direct observations. Microlensing particularly affects compact sources such as the X-ray emitting regions, the accretion disk, and the innermost broad-line emitting clouds. This phenomenon has emerged as an exceptionally powerful tool for probing the structure of lensed quasars (Jiménez-Vicente et al. 2012, 2014; Guerras et al. 2013a,b; Motta et al. 2012, 2017; Muñoz et al. 2016; Rojas et al. 2020). Specifically, the optical continuum light curves of lensed quasars at redshifts of z ∼ 1 − 2 can yield measurements of the radii for ultraviolet (UV) continuum-emitting sources surrounding the central supermassive black holes (SMBH; see Cornachione et al. 2020; Fian et al. 2016, 2018a, 2021a; Shalyapin et al. 2021; Rivera et al. 2023). Additionally, the wings of broad-emission lines (BELs) observed in quasar rest-frame UV spectra offer valuable insights into the structure and kinematics of the broad-line regions (BLR), as detailed in Hutsemékers & Sluse (2021), Hutsemékers et al. (2023), Fian et al. (2018b, 2021b, 2023a, 2024), and Popović et al. (2020). Unfortunately, many lensed quasars either lack sufficient high-quality data or only exhibit weak microlensing effects in their spectra and light curves. The doubly-lensed quasar SDSS J1339+1310, discovered by the Sloan Digital Sky Survey (Inada et al. 2009), is a notable exception. It is often referred to as a ‘microlensing factory’ (Shalyapin & Goicoechea 2014; Goicoechea & Shalyapin 2016) owing to its pronounced microlensing variations on short timescales in the optical light curves and substantial spectral distortions in several emission lines. The two optically bright images A and B of this quasar, with magnitudes of approximately r ∼ 18 − 19, are separated by $1\stackrel{″}{.}7$ and are located at a redshift of z_s = 2.231 (Shalyapin et al. 2021). An early-type galaxy G at redshift z_l = 0.607 (Shalyapin & Goicoechea 2014; Goicoechea & Shalyapin 2016) is located $0\stackrel{″}{.}63$ from image B (Inada et al. 2009) and serves as gravitational lens. It is believed that stellar mass objects within galaxy G act as microlenses and strongly affect image B (Shalyapin & Goicoechea 2014). Systems that exhibit highly variable microlensing events, like this one, are excellent tools for in-depth studies of the central engines of AGN and represent the best targets for subsequent follow-up investigations using suitable facilities.

In this paper, we conduct a thorough analysis of microlensing signals in both the wavelength and time domains for SDSS J1339+1310 to gain deeper insights into the quasar’s inner workings. In Section 2, we show the photometric and spectroscopic data compiled from various literature sources. Section 3 is dedicated to the analysis and interpretation of the data. Section 4 describes our microlensing simulations and explains the Bayesian methods we use to infer the sizes of emitting regions. In Section 5, we present and discuss our results, encompassing the structure of the continuum- and broad-line-emitting regions, as well as the mass estimation of the central SMBH. Our conclusions are summarized in Section 6. Throughout this work we assume the cosmological parameters Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 72 km s⁻¹ Mpc⁻¹.

2. Data and observations

From the literature, we have compiled a set of seven rest-frame UV spectra for images A and B of the gravitationally lensed quasar SDSS J1339+1310, covering a decade from May 2007 to April 2017. These spectra encompass a range of high- and low-ionization lines typical in quasars, such as Lyαλ1216, Si IV λ1397, C IV λ1549, C III] λ1909, and Mg II λ2798. The data, which have been fully reduced, are detailed in Table 1, along with information about the observations and corresponding references. In Figure 1, we display superpositions of the emission line profiles for Lyα, Si IV, C IV, C III], and Mg II corresponding to images A and B across various epochs. In Figure 2, we present the average profiles of these BELs for each lensed image. Individual epochs exhibiting noisy data within the emission line range were omitted before calculating the average spectra. We want to emphasize that in Figures 1 and 2, the emission line profile differences are highlighted by subtracting a linear fit of the continuum and normalizing the spectra to the peak of each respective emission line.

Fig. 1. Emission line profiles of images A (black) and B (blue) at seven different epochs in the rest-frame after the continuum has been subtracted and the line cores have been normalized. The observations reveal significant deformations in the red wings of image B for all emission lines.

Fig. 2. Average Lyα, Si IV, C IV, C III], and Mg II emission line profiles (from left to right) of the images A (black) and B (blue). Blue and red shaded regions indicate the selected windows for magnitude difference calculations, spanning 3000 − 8000 km s−1. We note that the blue wings of Lyα and Si IV are affected by absorption features and have been excluded from our analysis. The observations reveal significant microlensing-induced differences between the images in the red wings of all emission lines. The x-axis is represented on a velocity scale.

A visual inspection of the data reveals significant deformations in the emission line wings when comparing images A and B. While the blue wings of the high-ionization lines coincide well with each other, a pronounced enhancement is evident in the red wing of B compared to A across all emission lines. These observed anomalies are likely attributable to microlensing and indicate that the broad-line emission region has an anisotropic structure, as supported by findings in previous studies (e.g., Schneider & Wambsganss 1990; Abajas et al. 2002; Lewis & Ibata 2004; Sluse et al. 2012; Fian et al. 2018b, 2021b, 2024). We also observe time-varying distortions in the shelf-like feature located blueward of He II (around ∼λ1610), manifesting with varying intensities across different epochs. The origin of this feature remains uncertain; potential explanations include an extreme C IV red wing or He II blue wing, or possibly an unidentified species as discussed in Fine et al. (2010) and references therein. The blue wings of the low-ionization lines C III] and Mg II are contaminated by other emission features – Al III and Si III] in the case of C III], and the ‘small blue bump’ formed by the blending of various iron lines in the case of Mg II – both exhibiting indications of microlensing.

Following the discovery of SDSS J1339+1310 (Inada et al. 2009), the system was also monitored in the r-band within the framework of the Gravitational LENses and DArk MAtter (GLENDAMA) project¹ (Gil-Merino et al. 2018). The photometric monitoring was conducted with the 2.0 m Liverpool Telescope from 2009 to 2019, comprising 241 observing epochs. The light curves, as shown in the upper panel of Figure 3, are publicly available on the GLENDAMA website. For a comprehensive analysis of the data, we refer the reader to the detailed study conducted by Shalyapin et al. (2021).

Fig. 3. Light curves and differential microlensing variability of the quasar SDSS J1339+1310. Upper panel: Light curves of the lensed images A (black) and B (light blue), as obtained by the GLENDAMA project (refer to Gil-Merino et al. 2018). The light curve of image B is shifted by its time delay (−48 days) relative to image A. Lower panel: Differential microlensing variability in light curve B, compared against a polynomial fit to light curve A. The differential microlensing variability has been corrected for the magnitude difference between the images, using the Mg II line core flux ratio as a baseline for no microlensing. The differential microlensing measurements in the continuum at rest-frame 2000 Å, chosen for its close wavelength proximity to the r-band at 1930 Å, are depicted as blue triangles for comparison. We note that one of the five available spectroscopic epochs has been discarded due to excessive noise in the data. The horizontal axes display Julian (bottom) and Gregorian (top) dates.

Lensed quasar images exhibit a time delay in their arrival as they traverse different paths through the universe. Shalyapin et al. (2021) determined the time delay between the images to be Δt_AB = 48 ± 2 days, with image A leading. This measurement is crucial for effectively distinguishing microlensing-induced variability from the intrinsic variability of the source quasar. From Figure 3 we see that the quasar images display considerable time variability, with changes reaching up to ∼1 mag within periods of approximately two years. Images A and B exhibit parallel variations over long timescales, indicating that these long-term fluctuations primarily originate from the distant quasar itself, suggesting an intrinsic origin (see Goicoechea & Shalyapin 2016). Furthermore, there is evidence of extrinsic variability (i.e., microlensing) occurring on various timescales, as illustrated by several sharp events in image B that are not observed in image A.

3. Data analysis

3.1. Core ratio measurements

The cores of emission lines presumably originate from spatially very extended regions (narrow-line region and outer parts of the BLR) that remain largely unaffected by microlensing and intrinsic variability, as demonstrated in studies by Guerras et al. (2013a) and Fian et al. (2018b). The methodology of quantifying the flux ratios of the line cores between lensed images presents a valid approach for establishing a baseline that is relatively free of microlensing effects. This baseline is essential for accurately assessing (microlensing-based) sizes of emitting regions in quasars. Our study focuses on analyzing three high-ionization lines (including Lyα, Si IV, and C IV) alongside two low-ionization lines (C III] and Mg II). For each image and emission line, we remove the continuum by fitting and then subtracting a straight line from the continuum on both sides of the emission line. Given the variable widths of the emission lines, we employ windows of different widths to estimate the continuum for each line and quasar image, deliberately avoiding areas with known emission features. We then define the core fluxes using a narrow interval spanning between ±5 to 8 Å centered around the lines’ peaks (i.e., λ1216, λ1397, λ1549, λ1909, and λ2798). In Table 2, we present the average core flux ratios (±1σ standard deviation) between images A and B calculated from all spectra. Notably, our inferred core flux ratios are in excellent agreement with the values reported by Shalyapin & Goicoechea (2014) and Goicoechea & Shalyapin (2016). While it is presumed that the central features of the lines are dominated by photons arising from the narrow-line emission region and the outer parts of the BLR, some of them may also contain a significant number of photons emitted by clouds within the inner zones of the BLR. These clouds are moving almost perpendicularly to the line of sight, exhibiting small projected motions (see, e.g., Goicoechea & Shalyapin 2016). Goicoechea & Shalyapin (2016) reported no evidence of microlensing effects in the cores of the Mg II lines, a finding corroborated by our study. However, the core of C IV is clearly affected by microlensing, and to a certain extent, microlensing also impacts the cores of the other high-ionization lines studied in this work.

3.2. Continuum variability in spectroscopic data

The continuum emission of quasars primarily originates from a nuclear region, although recently several studies (see, e.g., Chelouche et al. 2019; Korista & Goad 2019; Fian et al. 2023b) offer a slightly different perspective. Microlensing magnification ratios of the continuum at various wavelengths yield crucial insights into the structure and kinematics of the continuum-emitting region. This is because microlensing is sensitive to the size of the source region, with smaller regions leading to higher magnifications. To accurately assess the impact of microlensing on the continuum, it is essential to distinguish this effect from the influences of macro-lensing magnification, which arises from the lens galaxy’s smooth potential, and extinction. While each lensed quasar image shares an intrinsic spectrum, they undergo different levels of extinction. This variation occurs because the light from each image traverses distinct paths through the lens galaxy, encountering different concentrations of dust and gas (Falco et al. 1999; Motta et al. 2002; Muñoz et al. 2004, 2011). Moreover, the macro-magnification and the differential extinction between the lensed quasar’s components are not influenced by the source size and affect both the continuum and emission line flux ratios (Motta et al. 2002; Mediavilla et al. 2005, 2009, 2011a). We attempt to correct for these effects by estimating the offsets between the continuum adjacent to the emission lines, (m_B − m_A)_cont, and the magnitude differences in the cores of the emission lines, (m_B − m_A)_core, between images A and B, μ_cont = (m_B − m_A)_cont − (m_B − m_A)_core. The cores of the emission lines, originating from material distributed across a broad area (narrow-line region and outer regions of the BLR), are typically extensive enough to remain unaffected by microlensing caused by solar mass objects (Guerras et al. 2013a; Fian et al. 2018b). It is important to note that the BLR is likely stratified, meaning different atoms or ions are located at varying distances from the nuclear continuum-emitting regions. For instance, there is evidence of a correlation between the degree of ionization and proximity to the continuum source, with high-ionization lines located much closer than low-ionization ones (e.g., Krolik et al. 1991; Peterson & Wandel 2000; Gaskell et al. 2007). In this stratification model of the BLR, the cores of the high-ionization lines might be subject to microlensing. Therefore, caution is advised when interpreting the measured signals in those lines. Given the relatively minor wavelength differences between the line cores and the selected wavelength intervals for the continuum (see Table 3), this estimator effectively eliminates the influences of macro-magnification and extinction, as demonstrated by Guerras et al. (2013a).

In Table 3, we present the central values employed for fitting the continuum, along with the average magnitude difference (including the ±1σ standard deviation) at each respective wavelength. We chose image A as the reference image as it exhibits less microlensing variability when compared to image B (see Shalyapin & Goicoechea 2014; Goicoechea & Shalyapin 2016). The values of μ_cont indicate that B is amplified relative to A by a factor of approximately 3 − 6, with greater amplifications at shorter wavelengths. This is in excellent agreement with the findings of Shalyapin & Goicoechea (2014), who reported a factor of 3 − 5, and Goicoechea & Shalyapin (2016), who estimated a mean magnification ratio of ∼5.7. The microlensing signal observed in the continuum of SDSS J1339+1310 is notably more pronounced than those typically seen in the majority of lensed quasars (Fian et al. 2016, 2018a, 2021a; Rojas et al. 2020; Cornachione et al. 2020; Muñoz et al. 2022).

3.3. Continuum variability in photometric data

Intrinsic variability of the source, combined with the light path time delay between the quasar images, could erroneously mimic anomalies in the optical flux ratios. Photometric monitoring of continuum variability in gravitationally lensed quasars offers the benefit of distinguishing between intrinsic and extrinsic variability, once the time delay is known. To analyze extrinsic variations in the light curves of SDSS J1339+1310, we initially model and subtract the variability intrinsic to the quasar itself. We applied the latest time delay estimates of Shalyapin et al. (2021) to shift the light curves of the lagging image B by −48 days. Subsequently, we corrected for the magnitude difference between the images by using the core ratio of the low-ionization line Mg II. We note that while the C III] emission line core is proximal to the r-band, it is likely contaminated by adjacent emission lines. Consequently, we selected the Mg II line core (B/A = 1.5 ± 0.2) as a reference, which is assumed to offer an uncontaminated baseline for assessing the true magnification ratio of the images in the absence of microlensing. We performed a single spline fitting on the light curve of image A, which is less susceptible to microlensing, to model the quasar’s intrinsic variability. Subsequently, by subtracting this spline fitting from the time-shifted light curve of image B, we constructed a difference light curve in which only the microlensing variability remains. During the decade of photometric monitoring, image B appears to have experienced several microlensing events with timescales of 50 − 100 days, as observable in the lower panel of Figure 3. While we rely on the spline fit to approximate the intrinsic variability of the source, we note that some microlensing variability in image A is likely still present. However, this is irrelevant in our analysis, as we account for the contributions of microlensing of image A in the simulated microlensing difference histograms.

3.4. BLR variability measurements

To assess the (minimal) size of the BEL emitting regions, we normalize the continuum-subtracted spectra across all images and epochs to match the core of the emission line, determined by the flux within a narrow range centered around the line’s peak. Assuming that the line cores serve as a reference that is little affected by microlensing and intrinsic variability, comparing the line wing fluxes between image pairs enables us to estimate the lower limit of the emitting region’s size. Distinguishing microlensing from intrinsic variability with certainty requires observations separated by the time delay between images. Nevertheless, given the short time delay between images A and B (48 ± 2 days), we posit that the impact of intrinsic variability on the observed magnitude differences is minimal. This assumption makes it reasonable to attribute most of the observed B–A magnitude differences in the continuum and BEL wings primarily to microlensing. We estimate the microlensing in the line wings, μ_wing = (m_B − m_A)_wing − (m_B − m_A)_core, on either side of the emission line peak, corresponding to a velocity range of 3000 − 8000 km s⁻¹. The velocity ranges are marked in Figure 2 and the differential microlensing measurements are listed in Table 4. We exclude calculations in line wings that show absorption features, such as the Lyα forest affecting the blue wing of the Lyα line and the pronounced absorption lines in the blue wing of Si IV.

Significant microlensing-induced variations are evident in the red wings of all emission lines in image B, as well as in the blue wing of C III]. The findings related to C III] warrant careful interpretation, considering that its blue wing is blended with Si III], Al III, and Ne III/Si II. It is noteworthy that the microlensing magnification of image B relative to image A achieves remarkable values of 2 − 3 in the enhanced red wing of C IV, in excellent agreement with the findings of Goicoechea & Shalyapin (2016).

3.5. UV Fe III gravitational redshift

Recently, Mediavilla et al. (2018) observed that the Fe III λλ2039−2113 feature in quasars is systematically redshifted, a phenomenon that can be explained by the gravitational redshift induced by the central SMBH. This feature is relatively free of contamination by lines from other species, and based on the observed impact of microlensing magnification, it is inferred to originate from an inner region of the BLR (Fian et al. 2018b, 2021b). We have modeled the Fe III λλ2039 − 2113 emission feature in the average spectra of images A and B, constructed from five different epochs encompassing the C III] to Mg II wavelength range. First, we fit straight lines to the continuum, one on each side of the C III] emission line and another across two windows flanking the blue and red sides of the Fe III blend. To accommodate variations in the continuum’s shape, which can change from epoch to epoch, we have adjusted the windows used for the continuum fitting where necessary. After subtracting the continuum from the spectrum, we normalize the resulting continuum-subtracted spectra to match the core of the C III] line, defined by the flux within a narrow interval centered on the line’s peak. The processed spectra for images A and B in the wavelength regions around the C III] line and the Fe III feature are displayed in Figure 4. Significant spectral shape deformations are present in the Fe III λλ2039−2113 blend redward of C III], where image B appears to be globally affected by high-magnification microlensing.

Fig. 4. Average spectra for images A (in black) and B (in blue), focusing on the regions around the C III] emission line and the Fe III λλ2039−2113 emission feature. The continuum has been subtracted, and the spectra normalized to the core of the C III] line. Vertical dotted blue lines mark the wavelengths of the Fe III lines as per the Vestergaard & Wilkes (2001) template. Dashed red and orange curves represent fits to the Fe III λλ2039−2113 blend in images A and B, respectively, shifted by approximately ∼8 Å and ∼12 Å to match the template rest-frame. The ordinate is in arbitrary flux units.

To model the Fe III λλ2039 − 2113 blend, we utilized a template comprising 19 distinct Fe III lines with fixed relative amplitudes, as provided by Vestergaard & Wilkes (2001). Each line in the template is shaped as a Gaussian, broadened uniformly to the same width, defined by σ (equivalent to FWHM/2.35). Additionally, the lines are collectively shifted in wavelength by a value of Δλ and scaled in amplitude by a factor of K_scale. These three parameters – σ, Δλ, and K_scale – are the free variables employed to fit the blend in both images. Figure 4 demonstrates that the Vestergaard & Wilkes (2001) template effectively reproduces the observed shape of the Fe III λλ2039−2113 blend in both images when a redshift of 8.0 Å for image A and 11.7 Å for image B, respectively, is applied. The average systematic redshift, ⟨Δλ⟩ = 9.8 ± 3.1 Å, aligns well with the findings of Mediavilla et al. (2018), who reported a value of ⟨Δλ⟩ = 10.3 ± 5.9 Å.

To assess the microlensing-induced changes in the UV Fe III emission line blend, we shift the wavelength window λλ2039 − 2113 based on the estimated redshift applicable to that specific epoch and image. We then measure the magnitude difference by comparing the fits to the UV iron blends observed in images A and B within this shifted wavelength range. This analysis yields an average value of μ_Fe III = 0.8 ± 0.2 mag.

4. Microlensing simulations

4.1. Magnification maps

For simulating the microlensing effects of extended sources, we used the Fast Multipole Method – Inverse Polygon Mapping (FMM–IPM) algorithm², as detailed in Jiménez-Vicente & Mediavilla (2022). This innovative approach integrates the FMM algorithm from Greengard & Rokhlin (1987) for computing ray deflections with the IPM algorithm from Mediavilla et al. (2006); Mediavilla et al. (2011b) for generating magnification maps. Our simulations are based on 3000 × 3000 pixel² maps, covering an area of approximately $204\sqrt{M/0.3M_{\odot }}\times 204\sqrt{M/0.3M_{\odot }}$ light-days² on the source plane, where M represents the microlensing mass. For SDSS J1339+1310, the value of the Einstein radius is $R_{\mathrm{E}}=2.64\times {10}^{16}\sqrt{M/0.3M_{\odot }}\mathrm{cm}=10.2\sqrt{M/0.3M_{\odot }}$ light-days at the source plane. Our magnification maps for images A and B feature a high resolution of $0.068\sqrt{M/0.3M_{\odot }}$ light-days per pixel, which effectively samples the quasar’s optical accretion disk. These maps are defined by two key parameters: the local shear, γ, and the local convergence, κ, where κ is directly proportional to the surface mass density. The local convergence comprises two distinct components: κ = κ_c + κ_⋆, where κ_c denotes the convergence attributed to continuously distributed matter, such as dark matter, while κ_⋆ corresponds to the convergence caused by stellar-mass point lenses, like microlens stars within the galaxy. The microlenses adhere to a power-law mass function, characterized by a slope similar to that of the Salpeter initial mass function (Salpeter 1955). The value of α ≡ κ_⋆/κ represents the fraction of mass in compact objects and quantifies the contribution of stars to the total mass in the lens galaxy. The parameters κ, γ, and α for images A and B were taken from Shalyapin et al. (2021) and are listed in Table 5 for reference. We note that Shalyapin et al. (2021) provided ten lens models that match observational constraints with a χ² ∼ 0. From these, we selected the model that most accurately represents the observed time delay between the images.

4.2. Source profile

To model the structure of the unresolved quasar, we adopt circular Gaussian profiles ($I(R)\propto exp(-R^2/2r_{\mathrm{s}}^2)$) to simulate the luminosity distribution of the emitting regions (see, e.g., Fian et al. 2018b, 2021b). To determine the magnifications for a source of size r_s, we convolve the magnification maps with Gaussian profiles having a sigma of r_s. The consensus in the field is that the specific shape of the source’s emission profile has a minimal impact on microlensing flux variability studies, as demonstrated by Mortonson et al. (2005) and Muñoz et al. (2016). The results are mainly controlled by the half-light radius rather than the detailed intensity profile. Several studies have utilized Gaussian intensity profiles to model the BLR in microlensing studies of quasars (see Motta et al. 2017; Guerras et al. 2013a; Sluse et al. 2011; Wayth et al. 2005). These profiles are favored due to their simplicity and effectiveness when convolving them with magnification maps. Consequently, we perform convolutions of the maps with Gaussians of varying sizes of r_s, spanning an interval between 0.1 and 18 light-days for a mean stellar mass ⟨M⟩ = 0.3 M_⊙. All linear sizes can be appropriately rescaled for a different mean stellar mass, following the relation $r_{\mathrm{s}}\propto \sqrt{⟨M⟩}$. The characteristic size r_s is related to the half-light radius by R_1/2 = 1.18 r_s for Gaussian profiles. After the convolution process, we normalize each magnification map by dividing it by its average value. The resulting histograms of these normalized maps depict the expected microlensing variability. Subsequently, we create the microlensing difference histograms B–A for various r_s values. These simulated microlensing difference histograms are then set for comparison with the observational data, as elaborated in Section 4.3.

4.3. Bayesian source size estimation

Based on the differential microlensing observed in the wings and adjacent continua of various emission lines between lensed images, we can deduce the size of the emission regions. For this purpose, a statistical approach was employed, treating each microlensing measurement as an individual epoch event. Utilizing all the available observational epochs, we calculate the joint microlensing probability, P(r_s), to derive an average size estimate, following the procedures described in Guerras et al. (2013a) and Fian et al. (2018b, 2021b). In the case of the microlensing difference light curve, which provides numerous microlensing measurements, we adopt a more advanced approach. This method entails comparing the microlensing histograms derived from the model, which correspond to convolutions with various source sizes r_s, with the histogram of the observed data. For a more detailed explanation of this technique, we refer to previous works such as Fian et al. (2016, 2018a, 2021a). In this study, to accommodate the uncertainties associated with the selected baseline for no microlensing, we permitted the observed data histogram to shift within these uncertainties (see Table 2). We have also examined the effects of transverse velocity on the correspondence between the pixels in the simulated magnification map and the observed time span. The light curves utilized cover a span of 10.1 years and comprise 241 data points, with the respective magnification maps displaying a resolution of $0.068\sqrt{M/0.3M_{\odot }}$ light-days per pixel. Assuming an effective transverse velocity of 1000 km s⁻¹ and uniform sampling across the light curves, we calculated that approximately 1.3 data points are associated with a single pixel. This does not cause significant alterations to the distribution of the observed histogram.

5. Results and discussion

Given the sensitivity of microlensing to the size of the source, we use our measurements of microlensing magnification amplitudes to estimate the size of the continuum-emitting region across various wavelengths, as well as the dimensions of the regions emitting the blue and red wings of the BELs and the UV Fe III blend in the lensed quasar SDSS J1339+1310.

5.1. Continuum-emitting region structure

We studied the structure of the lensed quasar’s continuum-emitting region through the statistics of microlensing magnifications, as inferred from rest-frame UV spectral data and photometric monitoring analyses. Utilizing Bayesian methods, as detailed in Section 4.3, we estimated the probability of r_s based on the observed differential microlensing between images A and B. Using spectroscopic data, this approach yielded half-light radii for the continuum source ranging from 0.6 to 1.2 light-days across the wavelength range of 1150 Å to 2850 Å (see Table 3). Analysis of the pronounced microlensing variations in the 2009–2019 r-band light curves allowed us to constrain the half-light radius of the UV continuum-emitting region at 1930 Å. By comparing the microlensing difference histogram derived from the data with the simulated microlensing histograms for sources with varying r_s values (see Figure 5), we determined a half-light radius of R_1/2 = 2.2 ± 0.3 light-days for the region emitting the r-band continuum. We observed that shifting the observational histogram by −0.07 mag yielded the optimal alignment with the model histograms. Our result is in agreement within uncertainties with the recent findings of Shalyapin et al. (2021), who reported a size of $R_{1/2}=1.0_{-0.6}^{+7.3}$ light-days.

Fig. 5. Microlensing frequency distributions derived from the observed light curves (illustrated with a blue-colored histogram) compared to those obtained from simulated microlensing magnification maps (depicted by gray and black lines). The polygonal lines represent model histograms for various values of R1/2 ranging from 2.0 to 3.0 light-days. The thick line highlights the model’s best fit to the observed data, with R1/2 = 2.2 ± 0.3 light-days. We note that the observational histogram was shifted by −0.07 mag to achieve the best alignment with the model histograms.

Figure 6 displays the microlensing-based UV continuum-emitting region sizes as a function of wavelength. We fit our size spectrum using a disk model r_th(λ)∝λ^β, where the power-law index β – a marker of the disk’s temperature gradient – is adjustable. From Figure 6, it is evident that both the estimated sizes and the physical model roughly align with the expectations of an optically thick geometrically thin accretion disk (see, e.g., Shakura & Sunyaev 1973), where β = 4/3. However, our findings suggest a slightly flatter slope than this theoretical model predicts, with the best fit achieved when β = 0.9 ± 0.2. This observation of a flatter size-wavelength relation is consistent with several microlensing campaigns, such as those reported by Morgan et al. (2010), Blackburne et al. (2011), Jiménez-Vicente et al. (2014), and Muñoz et al. (2016). We note that the r-band emitting region, marked by a blue point in Figure 6, was excluded from the size spectrum fitting. This decision was based on the fact that this region encompasses a portion of the broad C III] emission line, resulting in a bias towards larger sizes. Additionally, it should be noted that the results regarding the continuum adjacent to the C III] should be interpreted with caution as the presence of (highly variable) Fe II and Fe III lines might influence the continuum blueward of this line.

Fig. 6. Microlensing-based continuum-emitting sizes as a function of wavelength. The dashed black line shows the best fit to the data, with a power-law index of β ∼ 0.9. The red dotted line is a fit with a fixed theoretical power-law index of β = 4/3, as expected for an optically thick and geometrically thin accretion disk. The blue-colored point corresponds to the r-band continuum, and the shaded blue area illustrates the wavelength range covered by the r-band.

5.2. BLR structure

We applied the methodologies outlined in Section 5.1 to estimate the dimensions of the continuum-emitting regions, this time focusing on the wings of the BELs and the Fe III emission blend. In contrast to the nuclear continuum emission region, the wings of the BELs originate from a more expanded area located at a greater distance from the central SMBH. We note that due to the limitations of the available spectroscopic observations, we did not estimate the size of the emission regions of the BELs separately. Instead, we provide a joint analysis with rough size estimates for the average emission region sizes of the blue and red line wings. In Figure 7, we display the probability density functions (PDFs) associated with the regions emitting the blue and red wings of the BELs, as well as the region responsible for the UV Fe III emission. Intriguingly, our analysis reveals a notably small average size for the region emitting the red wing of the BELs (R_1/2 = 2.9 ± 0.6 light-days), approximately 30% larger than the size of the r-band emitting region. Additionally, the blue BEL wings in this system also originate from a relatively compact area (R_1/2 = 11.5 ± 1.7 light-days), smaller than the average BLR size estimates reported by Guerras et al. (2013a) and Fian et al. (2018b, 2021b) for a sample of lensed quasars (a few tens of light-days). The observed difference between the red and blue wings may be attributed to the velocity field of the red wing residing in a region of higher magnification compared to the blue wing. In single-epoch observations, regions of high magnification yield smaller size estimates, whereas regions of low magnification result in larger size estimates.

Fig. 7. PDFs of the half-light radius R1/2 emitting the blue and red BEL wings, and the Fe III λλ2039−2113 blend. The vertical dashed lines indicate the expected size of the emission region, while the shaded regions represent the one-sigma intervals.

Based on the differential microlensing estimates for the UV Fe III blend between images A and B, we deduce that the size of its emission region is R_1/2 = 3.1 ± 1.2 light-days across. This indicates that the Fe III feature is formed close to the accretion disk, aligning with the conclusions drawn in studies by Guerras et al. (2013b), as well as in our previous works (Fian et al. 2021b, 2022).

5.3. SMBH mass

Assuming that the redshift of the Fe III λλ2039−2113 lines is gravitational in origin, we can use the previously determined redshift of the UV Fe III blend together with the inferred microlensing-based size of its emitting region to derive the mass of the central SMBH. This methodology has been successfully applied to diverse quasar samples yielding mass estimates that are statistically compatible with virial masses. Specifically, it was deployed for a sample of ten lensed quasars (Mediavilla et al. 2018), extended to another set of ten non-lensed quasars (Mediavilla et al. 2019), and employed in a detailed study of Q 0957+561 (Fian et al. 2022). Although this redshift method is generally unaffected by geometry and non-gravitational forces, it is nonetheless reliant on a size estimator expressed as:

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=\frac{2c^2}{3G}{\left(\frac{\mathrm{\Delta }\lambda }{\lambda }\right)}_{\mathrm{Fe}\mathrm{III}}{R}_{\mathrm{Fe}\mathrm{III}},\end{array}$$(1)

where Δλ/λ represents the redshift of the Fe III λλ2039−2113 emission lines, and R_Fe III denotes the size of the region emitting this blend. Inserting the average redshift, ⟨Δλ⟩, and the half-light radius of the UV Fe III blend, R_1/2 = R_Fe III, in Eq. (1), we obtain a SMBH mass of M_BH = (1.7 ± 0.5)×10⁸ M_⊙. This result is in agreement within uncertainties with the virial mass estimate provided by Shalyapin et al. (2021), $M_{\mathrm{BH}}=4.0_{-2.4}^{+6.0}\times {10}^8{M}_{\odot }$. The uncertainties in our SMBH mass estimate are proportional to the square root of the combined squared relative errors of the Fe III emitting region size determined via a maximum likelihood method, and the squared standard error of the mean redshift.

6. Conclusions

SDSS J1339+1339 is, to date, one of the most heavily microlensed quasars, and a substantial amount of photometric monitoring data and spectroscopic observations have been collected over the past decade. In this study, we utilized the newly published macro-model of the lens system by Shalyapin et al. (2021), alongside cutting-edge techniques to quantitatively model the microlensing statistics (see Jiménez-Vicente & Mediavilla 2022).

Firstly, it is worth highlighting that SDSS J1339+1339 represents an atypical system, characterized by variations in its light curves occurring over remarkably short timescales. These fluctuations could be attributed to a high transverse velocity and/or a low mass of the microlenses involved. Alternatively, the presence of hotspots or substructures within the accretion disk may also account for these observations. Such features could undergo (de-)magnification if image B crosses one or multiple caustic regions. By analyzing the rapid, microlensing-induced variations in the r-band light curves, along with the chromatic changes in the continua adjacent to the BELs, we gained insights into the structure of the continuum-emitting region. The microlensing signals indicate that image B is magnified ∼3 − 6 times relative to image A at wavelengths between ∼1000 and 3000 Å, exhibiting a trend of higher amplification at shorter wavelengths. This finding suggests that the sizes of the continuum-emitting regions increase with wavelength, lending support to the hypothesis of disk reprocessing. Additionally, it is noteworthy that the half-light radius we inferred for the region emitting the r-band continuum (R_1/2 = 2.2 ± 0.3 light-days) is consistent within uncertainties with the value reported by Shalyapin et al. (2021) ($R_{1/2}=1.0_{-0.6}^{+7.3}$ light-days). However, the size of the emission region in the r-band is overestimated by about three times when compared to the sizes obtained from spectroscopic data for the continuum-emitting regions. This overestimation is likely because the r-band filter includes part of the BEL C III], which is less affected by microlensing and therefore results in a bias towards a larger emission region size.

Through the assembly of a dataset comprising seven spectroscopic observations, we have conducted an in-depth analysis of various emission lines, including Lyα, Si IV, C IV, C III], and Mg II, as well as the Fe III λλ2039 − 2113 emission feature located redward of C III]. In our comparative study of the emission line shapes between images A and B, we have identified noticeable distortions attributable to microlensing. While the blue wings of these emission lines coincide well with each other, a pronounced enhancement is observed in the red wings of image B when compared to that of image A. The observation of microlensing-induced asymmetric deformations in all studied BELs provides strong evidence that the geometry of the BLR is not spherical. Microlensing of a spherically symmetric BLR would lead, in general, to symmetric variations in the emission lines. However, when the BLR is non-spherical, such as a disk-like structure, the microlensing effect is asymmetric, leading to asymmetric variations in the emission lines (Schneider & Wambsganss 1990; Lewis & Ibata 2004; Abajas et al. 2002; Sluse et al. 2012). Using Bayesian analysis techniques, we have estimated the average dimensions of the regions responsible for emitting the blue and red wings of the BELs. Our findings reveal that the minimal average size of the region emitting the blue wings is 11.5 ± 1.7 light-days. In contrast, the red wings are emitted from a region approximately four times smaller, measuring 2.9 ± 0.6 light-days. These sizes are significantly smaller compared to the BEL emitting regions typically found in quasars, which, as documented in the literature, are generally within the tens of light-days range (see, e.g., Guerras et al. 2013a; Fian et al. 2018b, 2021b). However, similar findings of small sizes for the BLR have also been observed in another system, SDSS J1004+4112, as reported in recent studies by Hutsemékers et al. (2023) and Fian et al. (2024).

We have also deduced that the spectral feature associated with the UV Fe III emission originates from a compact region of 3.1 ± 1.2 light-days across. By combining the microlensing-based size estimate with the measurement of the gravitational redshift for this blend, we were able to calculate the mass of the central SMBH, yielding a value of M_BH = (1.7 ± 0.5)×10⁸ M_⊙. Our inferred mass is consistent within uncertainties with the virial black hole mass estimate for this system ($M_{\mathrm{BH}}=4.0_{-2.4}^{+6.0}\times {10}^8{M}_{\odot }$; see Shalyapin et al. 2021), as well as with the average black hole mass determined for a sample of 27 lensed quasars ($M_{\mathrm{BH}}=2.4_{-0.4}^{+1.5}\times {10}^8{M}_{\odot }$; see Fian et al. 2021b). This consistency reinforces the reliability of the (gravitational redshift plus microlensing-based size) mass measurement, making this technique a compelling alternative to other existing methods (e.g., the viral plus reverberation-based size) for measuring black hole masses.

Acknowledgments

We thank the referee James Chan for the valuable comments and suggestions. This research was supported by the grants PID2020-118687GB-C31, PID2020-118687GB-C32, and PID2020-118687GB-C33, financed by the Spanish Ministerio de Ciencia e Innovación. J.A.M. is also supported by the Generalitat Valenciana with the project of excellence Prometeo/2020/085. J.J.V. is also financed by the project FQM-108, financed by Junta de Andalucía. D.C. and S.K. are financially supported by the DFG grant HA3555-14/1 to Tel Aviv University and University of Haifa, and by the Israeli Science Foundation grant no. 2398/19 and 1650/23.

References

All Tables

All Figures

	Fig. 1. Emission line profiles of images A (black) and B (blue) at seven different epochs in the rest-frame after the continuum has been subtracted and the line cores have been normalized. The observations reveal significant deformations in the red wings of image B for all emission lines.
In the text

Fig. 2. Average Lyα, Si IV, C IV, C III], and Mg II emission line profiles (from left to right) of the images A (black) and B (blue). Blue and red shaded regions indicate the selected windows for magnitude difference calculations, spanning 3000 − 8000 km s−1. We note that the blue wings of Lyα and Si IV are affected by absorption features and have been excluded from our analysis. The observations reveal significant microlensing-induced differences between the images in the red wings of all emission lines. The x-axis is represented on a velocity scale.

In the text

Fig. 3. Light curves and differential microlensing variability of the quasar SDSS J1339+1310. Upper panel: Light curves of the lensed images A (black) and B (light blue), as obtained by the GLENDAMA project (refer to Gil-Merino et al. 2018). The light curve of image B is shifted by its time delay (−48 days) relative to image A. Lower panel: Differential microlensing variability in light curve B, compared against a polynomial fit to light curve A. The differential microlensing variability has been corrected for the magnitude difference between the images, using the Mg II line core flux ratio as a baseline for no microlensing. The differential microlensing measurements in the continuum at rest-frame 2000 Å, chosen for its close wavelength proximity to the r-band at 1930 Å, are depicted as blue triangles for comparison. We note that one of the five available spectroscopic epochs has been discarded due to excessive noise in the data. The horizontal axes display Julian (bottom) and Gregorian (top) dates.

In the text

Fig. 4. Average spectra for images A (in black) and B (in blue), focusing on the regions around the C III] emission line and the Fe III λλ2039−2113 emission feature. The continuum has been subtracted, and the spectra normalized to the core of the C III] line. Vertical dotted blue lines mark the wavelengths of the Fe III lines as per the Vestergaard & Wilkes (2001) template. Dashed red and orange curves represent fits to the Fe III λλ2039−2113 blend in images A and B, respectively, shifted by approximately ∼8 Å and ∼12 Å to match the template rest-frame. The ordinate is in arbitrary flux units.

In the text

Fig. 5. Microlensing frequency distributions derived from the observed light curves (illustrated with a blue-colored histogram) compared to those obtained from simulated microlensing magnification maps (depicted by gray and black lines). The polygonal lines represent model histograms for various values of R1/2 ranging from 2.0 to 3.0 light-days. The thick line highlights the model’s best fit to the observed data, with R1/2 = 2.2 ± 0.3 light-days. We note that the observational histogram was shifted by −0.07 mag to achieve the best alignment with the model histograms.

In the text

Fig. 6. Microlensing-based continuum-emitting sizes as a function of wavelength. The dashed black line shows the best fit to the data, with a power-law index of β ∼ 0.9. The red dotted line is a fit with a fixed theoretical power-law index of β = 4/3, as expected for an optically thick and geometrically thin accretion disk. The blue-colored point corresponds to the r-band continuum, and the shaded blue area illustrates the wavelength range covered by the r-band.

In the text

	Fig. 7. PDFs of the half-light radius R1/2 emitting the blue and red BEL wings, and the Fe III λλ2039−2113 blend. The vertical dashed lines indicate the expected size of the emission region, while the shaded regions represent the one-sigma intervals.
In the text