key: cord-0290829-fwashk96
authors: Williams, Peter R.; Treu, Tommaso; Dahle, Haakon; Valenti, Stefano; Abramson, Louis; Barth, Aaron J.; Brewer, Brendon J.; Dyrland, Karianne; Gladders, Michael; Horne, Keith; Sharon, Keren
title: Dynamical Modeling of the CIV Broad Line Region of the $z=2.805$ Multiply Imaged Quasar SDSS J2222+2745
date: 2021-03-19
journal: nan
DOI: 10.3847/2041-8213/ac081b
sha: c7a1310402b79e0ebc38a0935680082dfd828332
doc_id: 290829
cord_uid: fwashk96

We present CIV BLR modeling results for the multiply imaged $z=2.805$ quasar SDSS J2222+2745. Using data covering a 5.3 year baseline after accounting for gravitational time delays, we find models that can reproduce the observed emission-line spectra and integrated CIV fluctuations. The models suggest a thick disk BLR that is inclined by $sim$40 degrees to the observer's line of sight and with a emissivity weighted median radius of $r_{rm median} = 33.0^{+2.4}_{-2.1}$ light days. The kinematics are dominated by near-circular Keplerian motion with the remainder inflowing. The rest-frame lag one would measure from the models is $tau_{rm median} = 36.4^{+1.8}_{-1.8}$ days, which is consistent with measurements based on cross-correlation. We show a possible geometry and transfer function based on the model fits and find that the model-produced velocity-resolved lags are consistent with those from cross-correlation. We measure a black hole mass of $log_{10}(M_{rm BH}/M_odot) = 8.31^{+0.07}_{-0.06}$, which requires a scale factor of $log_{10}(f_{{rm mean},sigma}) = 0.20^{+0.09}_{-0.07}$.

Precise measurements of supermassive black hole (BH) masses across cosmic time are a necessary ingredient for understanding BH formation and growth and the connection between BHs and their host galaxies (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Ding et al. 2020) . The most successful technique for measuring M BH outside the local universe is reverberation mapping (RM, Blandford & McKee 1982; Peterson 1993 Peterson , 2014 Ferrarese & Ford 2005) , which measures the response of the broad emission-line region (BLR) to changes in the continuum. By combining emission-line widths with the time lag between continuum and emission-line fluctuations, one can obtain a virial estimate of the black hole's

The scale factor, f , accounts for the geometry, kinematics, and orientation of the BLR, which are intrinsic to each individual region and, in general, unknown. This is the largest source of uncertainty in RM M BH measurements, estimated to be ∼0.4 dex (Park et al. 2012) . Direct modeling of the BLR (Pancoast et al. , 2012 Brewer et al. 2011 ) avoids this issue entirely by including M BH as a free parameter. Additionally, models inform us of the structure and kinematics of the BLR and can provide f values for individual BLRs. This opens the possibility of finding correlations between f and other observables that could be used to improve M BH measurements for all active galactic nuclei (AGNs), not just those with data suitable for modeling (Williams et al. 2018 , Villafana et al. 2021 .

Until now, the dynamical modeling approach has been applied only to nearby AGNs with z < 0.1, and only one of these analyses (NGC 5548, Williams et al. 2020a) examined the UV-emitting BLR. This is due to the complexities of high-z reverberation mapping campaigns paired with the high signal-to-noise ratio (S/N ) spectra required for modeling. However, studies of BH growth require precise M BH measurements across all stages of the Universe, so an understanding of the UV BLR in the early Universe and in quasar-like environments is necessary.

The extraordinary data set of the multiply imaged quasar SDSS J2222+2745 (discovered by Dahle et al. 2013 ) described by Williams et al. (2020b, hereafter Paper I) is the first high-z monitoring campaign with data quality good enough for the modeling approach. In this paper, we model the C iv-emitting BLR in SDSS J2222+2745 using the data presented in Paper I. In Section 2, we briefly describe the data set and the BLR modeling approach used in the analysis. In Section 3, we present the model fits to the data and describe the inferred BLR geometry and kinematics. We also compute the scale factor f for the SDSS J2222+2745 C iv BLR and compare it to values for the Hβ BLRs of other AGNs determined using the same modeling approach. We conclude in Section 4. When necessary, we adopt a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

The data used in this analysis are the same data presented in Paper I, and the modeling approach is described in detail by Pancoast et al. (2014) . Here, we briefly summarize both components, but direct the reader to the respective papers for detailed explanations.

Beginning June 2016, we obtained monthly spectra of the three brightest images of SDSS J2222+2745 with the Multi-Object Spectrograph at Gemini Observatory North (GMOS-N; Hook et al. 2004 ). The spectra covered ∼5000 to 8200Å after dithering and flux calibration, covering the C iv and C iii] broad emission lines. In addition, we obtained g-band photometry with roughly twice-per-month cadence, beginning in September 2011 with the Alhambra Faint Object Spectrograph and Camera (ALFOSC) at the 2.56m Nordic Optical Telescope (NOT).

Since we are modeling only the emission of the C iv BLR, we first need to isolate the C iv broad emission line from the other spectral components. As described in Paper I, we model C iv as a fourth-order Gauss-Hermite polynomial, and we use those fits in this analysis. We combine the time-series of spectra for each image by first multiplying the fluxes by the corresponding image magnification and then shifting the times by the measured time delays (∆τ AB = −42.44 days and ∆τ AC = 696.65 days; Dyrland 2019), setting image A as the reference.

The 2020 spectra for the leading image C suffer from low S/N as a result of the small image magnification paired with relaxed observing condition constraints in 2020 due to the COVID-19 pandemic. For this reason, we remove these data from the analysis.

To speed up the modeling code, we also combine spectra that fall within 7 days of each other (observed frame), after accounting for gravitational time delays. This corresponds to ∼2 days in the rest frame, which is over an order of magnitude smaller than the expected BLR size. Since our smooth BLR model would be unable to resolve variations on this timescale, we lose no constraining power by combining these spectra. After making these changes, our final data set consists of 63 spectra covering a 1944 day (5.3 year) baseline, after accounting for gravitational time delays.

Finally, we bin the individual spectra by a factor of 8 in wavelength to further decrease computation time, giving 93 wavelength bins from 5503.5 to 6239.5Å. The smooth BLR model is unable to resolve emission-line structure on the 1Å/pix (50 km/s/pix) scale of the raw data, and we are using smooth Gauss-Hermite fits to C iv, so attempting to fit all pixels would introduce unnecessary computational burden while providing no additional constraining power.

We model the BLR emission using a collection of point particles surrounding the central black hole and ionizing source. The ionizing light propagates outwards and as it reaches the particles, they instantaneously re-process the light and emit it towards the observer in the form of emission lines.

The positions and velocities of the particles are determined by a number of free parameters, described in detail by Pancoast et al. (2014) . To summarize, the particles are distributed in a thick disk with half-opening angle θ o and inclined relative to the observer's line of sight by θ i (θ i = 0 deg is face-on). An additional parameter, γ, determines if the particles are distributed uniformly throughout the thick disk (γ = 1) or if they are concentrated at the faces of the disk (γ = 5). The distances of the particles from the origin are drawn from a Gamma distribution with shape parameter β, and mean µ, that has been shifted from the origin by a minimum radius r min .

A parameter κ determines the relative brightness of each particle by controlling if particles re-emit preferen-tially back towards the origin (κ = −0.5), isotropically (κ = 0), or away from the origin (κ = 0.5). The disk midplane can range from being fully opaque (ξ = 0) to fully transparent (ξ = 1).

A fraction of the particles, f ellip , are assigned to have near-circular orbits with radial and tangential velocities drawn from a distribution centered on the circular velocity. The remaining particles all have either inflowing or outflowing trajectories based on the binary parameter f flow (< 0.5 inflow, > 0.5 outflow). The radial and tangential velocities are drawn from a distribution centered on the radial escape velocity that is rotated by an angle θ e in the v r − v φ plane. This allows for particles with purely radial motion (θ e ∼ 0 deg) or on highly elliptical, bound orbits (θ e ∼ 45 deg).

When interpreting the model parameters, it is important to keep in mind that we are modeling the BLR emission rather than the underlying gas. We do not include in our model the complex photoionization process, which would require additional assumptions about the BLR environment-such as the gas density, temperature, and metallicity distributions-and the relation between the observed g-band continuum and the ionizing spectrum.

Using the observed continuum light curve, we can compute the time-series of emission-line spectra that a given model produces by summing the contributions of each BLR particle, taking into account the positioninduced time lag and velocity-induced wavelength shift. Our goal is to explore the model parameter space such that the model-produced spectra best fit the observed spectra. We assume that the ionizing continuum follows the observed g-band continuum and use Gaussian processes as a way to flexibly interpolate between observed data points.

We use a Gaussian likelihood function to compare the emission-line spectra with the model spectra. In postprocessing, we also introduce a statistical temperature, T , that softens the likelihood function, effectively increasing the spectra uncertainties. This accounts for under-estimates of uncertainties as well as the challenge of fitting a complex BLR with a simple model. To explore the parameter space of both the BLR model and continuum model, we use the diffusive nested sampling code DNest4 (Brewer & Foreman-Mackey 2016) . Diffusive nested sampling is a modification of the nested sampling technique that is particularly efficient at exploring complex, high-dimensional probability spaces. From the code output, we produce a posterior sample from which we can infer the model parameter values. 2 : One possible model fit to the observed spectra.

3 : Normalized residual ([Data − Model]/Data uncertainty). 4 : Observed spectrum from one of the epochs with uncertainties multiplied by √ T (black points) and the spectrum produced by the model shown in panel 2 (red). 5 : Integrated C iv emission-line flux of the data (black points) and model (red). 6 : Observed g-band continuum light curve (black) and the model fit to the light curve (red).

The modeling code is able to find regions of the BLR model parameter space that reproduce the observed emission-line shape and fluctuations (Figure 1) . To avoid over-fitting, we soften the likelihood function with a temperature T = 180 when constructing the posterior sample from the DNest4 output. This is equivalent to multiplying all uncertainties on the spectra by a factor of √ T = 13.4, which is reflected in the error bars in Panels 4 and 5. The model displayed in Figure 1 

Examining the posterior probability density functions (PDFs) for the model parameters, we find two clusters of solutions: those with θ i and θ o around 45 degrees and those with θ i and θ o around 85 degrees. Due to known degeneracies resulting from the flexible parametrization of the model (see, e.g., Grier et al. 2017) , different combinations of parameters can produce identical particle distributions and velocities. For θ i ∼ θ o → 90 deg, the particles are arranged in a sphere, but as γ → 5, they are increasingly concentrated along the axis of the sphere which is close to perpendicular to the observer's line-of-sight. Since RM data cannot resolve rotations in the plane of the sky, these solutions are equivalent to near-face-on thick disk models and are thus degenerate with the first family of solutions. The 2D posterior PDFs reveal that this is the case for the models in our sample with θ i , θ o ∼ 85. Unfortunately, the other model parameters are difficult to interpret in this arrangement, so we choose to exclude solutions with θ i > 65 deg.

Taking the median and 68% confidence intervals for each parameter, we find a thick disk BLR with θ o = 52.6 +5.5 −7.8 deg that is inclined by θ i = 40.8 +5.6 −6.5 deg. The parameter determining if emission is concentrated on the faces of the disk is not constrained for these models (γ = 2.6 +1.4 −1.2 ). The radial distribution of BLR emission drops off with radius faster than exponentially (β = 1.32 +0.11 −0.09 ), is shifted from the origin by a minimum radius r min = 8.2 +1.3 −1.1 ld, and has a median radius r median = 33.0 +2.4 −2.1 ld. This is visible in Figure 2 with the shell-like concentration of points around the BLR center and a rapid drop-off in point density at larger radii. The particles preferentially emit back towards the ionizing source (κ < −0.47), which agrees with photoionization predictions, and the disk midplane is found to be mostly transparent (ξ > 0.77). The modeling approach is more reliable since it does not depend on Hβ-to-C iv conversions for f , and is more precise since it avoids the additional uncertainty introduced by the intrinsic scatter in f .

The rest-frame time lag one would measure based on these models is τ median = 36.4 +1.8 −1.8 days, which agrees very well with the cross-correlation measurement of τ cen = 36.5 +2.9 −3.9 by Paper I. For the BLR model shown in Figure 3 , we also compute the mean emission-line lag in the five wavelength windows used to measure the velocity-resolved lags in Paper I. We should note that the lag for the central wavelength bin in Paper I was based on a narrow-band filter transmission curve, rather than the top-hat function used here. Since the transmission is higher near the center of the bin where the lags are longest, this will bias the filter-based measurements towards longer lags. Regardless, all of our model-based measurements are consistent with the velocity-resolved lags based on cross correlation.

In Figure 4 , we show the 2D posteriors for log 10 (M BH /M ), θ i , and θ o . There is an anti-correlation between the black hole mass and inclination angle, which is expected-as the BLR is tilted closer to face on, a larger black hole is required to reproduce the observed line width. We also find two streaks of solutions that can be separated by the line θ i = θ o . Those with θ i > θ o tend to have slightly more transparent midplanes and a higher fraction of particles on near-circular orbits, but all other parameter distributions remain indistinguishable. Importantly, both cases give the same black hole mass.

Using the M BH posterior PDF from our BLR model along with the C iv line widths and lags measured in Paper I, we can compute the appropriate scale factor, f , for the SDSS J2222+2745 BLR. We use the rest-frame lag, τ cen = 36.5 +2.9 −3.9 days, and the four emission-line widths: ∆V mean,FWHM = 7734 ± 59 km s −1 , ∆V mean,σ = 4261 ± 49 km s −1 , ∆V rms,FWHM = 9219 ± 458 km s −1 , and ∆V rms,σ = 5907 ± 148 km s −1 .

We propagate all uncertainties utilizing the full M BH posterior as follows: For each sample in the posterior, we draw an emission-line width from a normal distribution, N (∆V, σ 2 ∆V ), where ∆V is the median value reported above and σ ∆V is the corresponding uncertainty. We then draw a time lag using the same approach, assuming an uncertainty that is equal to the average of the upper and lower uncertainties, 3.4 days. From this, we construct a sample of log 10 f values and compute the median and 68% confidence intervals. We find log 10 (f mean,σ ) = 0.20 +0.09 log 10 (f mean,σ ) Note that the uncertainties on σ and τcen are not included in this conversion from log 10 (MBH/M ) to log 10 f , so the true posterior distribution for log 10 f is slightly broader.

log 10 (f mean,FWHM ) = −0.32 +0.09 −0.07 , log 10 (f rms,σ ) = −0.08 +0.09 −0.07 , and log 10 (f rms,FWHM ) = −0.46 +0.10 −0.09 . The uncertainties on log 10 f are larger than those for log 10 (M BH /M ) since they include the uncertainties on ∆V and τ .

In Figure 5 , we show the scale factor plotted against various BLR model and AGN parameters, along with the data and fits for the Hβ BLR by Williams et al. (2018) and the C iv and Hβ BLRs for NGC 5548 (Williams et al. 2020a) . To determine the bolometric luminosity, we use log 10 (L 1350 /erg s −1 ) = 44.66 ± 0.18 from Paper I and assume a bolometric correction of BC 1350 = 3.81, computed by Shen et al. (2011) using the composite quasar spectral energy distributions of Richards et al. (2006) .

Paper I also computed conversion factors between f Hβ and f CIV by comparing AGNs with line width and lag measurements available for both Hβ and C iv. They found that for a relation of the form log 10 (f CIV ) = Figure 5 . Correlations between the scale factor f and BLR model parameters and AGN properties. Each point and contour pair shows the median and 68% confidence region of the 2D posterior PDFs. The blue point is SDSS J2222+2745, orange (black) is the C iv (Hβ) BLR of NGC 5548 (Williams et al. 2020a) , and the grey points and fits are those for the Hβ BLR from Williams et al. (2018) . For the two C iv measurements, we also show a shifted contour (dashed) based on the f Hβ to fCIV conversions calculated in Paper I.

log 10 (f Hβ ) + α, the best fits are α mean,σ = 0.087 ± 0.007, α rms,σ = −0.021 ± 0.028, and α mean,FWHM = 0.694 ± 0.008. In Figure 5 , we also show the two C iv-based contours shifted by these conversion factors. An anti-correlation between f and θ i is expected to exist, regardless of the emission line or AGN properties. The parameter θ i simply describes the orientation of the BLR in the sky, so for a BLR model with disk-like kinematics, as θ i increases, f must decrease to account for the larger line-of-sight velocity component contributing to the emission-line width. The SDSS J2222+2745 measurement happens to fall on the same relation as the Hβ BLR measurements, but it is possible that C iv BLRs follow a different relation and this is a chance alignment. A larger sample of C iv BLR models is necessary to distinguish between the two options.

While θ i itself is not an observable, it may be deduced in AGNs with radio jets if the BLR rotation axis and jet are aligned. BLR modeling of 3C 120 (Grier et al. 2017) and the spatially resolved BLR of 3C 273 from interferometry (Gravity Collaboration et al. 2018 ) support that this is the case, although a larger sample size is still necessary to determine if this relation holds for every AGN.

We also find that the SDSS J2222+2745 measurement lies below the Hβ relations with M BH , and the offset is exaggerated by the f CIV to f Hβ correction factor. The Hβ relation between f mean,FWHM and M BH is currently detected only at the 2σ level (Williams et al. 2018) , so a larger sample will be necessary to solidify the relationship (Villafana et al. 2021, in prep) . Assuming that the relationship is real for the Hβ BLR, the offset of SDSS J2222+2745 could mean that the correction factors do not hold for all AGNs, or that there is something inherently different about the SDSS J2222+2745 BLR. One possibility is that C iv BLRs follow a different relationship than Hβ BLRs, although the NGC 5548 measurement is better aligned with the Hβ-based measurements. Another possibility is that BLR environments are different at z = 2.805 than at z < 0.1. Determining the source of the offset will require additional modeling of BLRs using multiple emission lines and across a wider range of redshifts.

The results presented here are the first of their kind for an AGN at z > 0.1 and at the peak of AGN activity. The tight constraints on M BH and other BLR properties demonstrate that gravitational lensing is a powerful tool that can provide data good enough for detailed BLR studies. While opportunities like that of SDSS J2222+2745 are rare, they serve an important role in expanding our understanding of the BLR outside of the local Universe.

The main results of our analysis can be summarized as follows:

1. The C iv-emitting BLR for SDSS J2222+2745 is a thick disk with a size r median = 33.0 +2.4 −2.1 light days. The kinematics are dominated by nearcircular Keplerian motion, with the remaining 10-20% of emission indicating infalling trajectories.

2. The median rest-frame lag produced by the BLR models is 36.4 +1.8 −1.8 days which agrees closely with the τ cen = 36.5 +2.9 −1.9 day measurement from cross correlation. The velocity-binned mean lags from the model are also consistent with the velocityresolved lags from cross correlation.

3. The black hole mass for SDSS J2222+2745 is log 10 (M BH /M ) = 8.31 +0.07 −0.06 . This corresponds to scale factors of log 10 (f mean,σ ) = 0.20 +0.09 −0.07 , log 10 (f mean,FWHM ) = −0.32 +0.09 −0.07 , log 10 (f rms,σ ) = −0.08 +0.09 −0.07 , and log 10 (f rms,FWHM ) = −0.46 +0.10 −0.09 .

Master's thesis

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

Astronomical Society of the Pacific Conference Series

The data presented here are based in part on observations obtained at the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOTSA. Partly based on observations made with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea. We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance.This research made use of Astropy, 1 a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013 PW and TT gratefully acknowledge support by the National Science Foundation through grant AST-1907208 (Tody 1986 (Tody , 1993 , Scipy (Virtanen et al. 2020)