The Sloan Digital Sky Survey Reverberation Mapping Project: How Broad Emission Line Widths Change When Luminosity Changes

To perform spectral fits on the multi-epoch spectra of our sample to measure broad-line properties, we adopt the publicly available quasar spectral fitting package QSOFIT (Shen et al., 2019). QSOFIT deploys a multi-component functional fitting approach similar to earlier work (e.g., Shen et al., 2008, 2011). Its continuum model consists of a power-law, an optical Fe ii template (Boroson & Green, 1992), an UV Fe ii template¹¹1The Fe ii template is the same as used in Shen & Liu (2012). See their §3.1 for more details. (Vestergaard & Wilkes, 2001; Tsuzuki et al., 2006; Salviander et al., 2007) and a 5th-order polynomial to account for any possible reddening. The emission line flux from the continuum and Fe ii subtracted spectrum is then fitted with multiple Gaussians. The fits were performed in the rest-frame of the quasar using systemic redshifts from Shen et al. (2019).

QSOFIT allows the user to specify the fitting range and switch on/off individual continuum and line components. To mitigate the difficulty of fitting the global continuum, we fit each spectrum locally in several line complexes by specifying the fitting range parameter. Table 3.1 summarizes the fitting ranges we use for each line complex, the continuum and line components included in the fitting, as well as the number of broad-line and narrow-line Gaussians used for each line component.

We now describe the multi-epoch spectral fitting procedure in detail. For each object, we first measure the median S/N over the spectrum for all spectroscopic epochs. Then the mean and standard deviation $\sigma$ of individual median S/N are computed in each observing season, and epochs that are 2$\sigma$ below the mean S/N of each season are rejected. This step excludes 2$\sim$6 epochs in the 2014 season, and 0$\sim$3 epochs in each of other three seasons. The remaining good S/N epochs are used to generate the mean spectrum using the SDSS-III spectroscopic pipeline idlspec2d. Next, we fit the mean spectra using QSOFIT. The high S/N mean spectrum is used to obtain robust estimates of the parameters describing the narrow line emission and the Fe ii emission, which are held fixed in the multi-epoch spectra fitting. The velocity division between narrow-line and broad-line Gaussian components is set to be 450 km/s in Gaussian $\sigma$. We tie the widths of the narrow lines within the same line complex (e.g., H$\beta$) but allow different narrow-line widths across different line complexes (see Table 3.1). Fixing the Fe ii emission is a simplification, since broad Fe ii does vary, albeit with smaller amplitude compared to other major broad lines (e.g., Barth et al., 2013; Hu et al., 2015). However, the individual epochs of spectra do not have sufficient S/N to allow reliable Fe ii subtraction, and therefore we use the average Fe ii emission constrained from the mean spectrum as a template to subtract this component.

As a demonstration, Figure 1 shows the multi-component functional fitting results of different line complexes from the mean spectra. Figure 2 shows the 2-D residual map to illustrate the quality of the multi-epoch fitting.

To quantify the line width, we calculate the Full-Width-at-Half-Maximum (FWHM) and the line dispersion $\sigma_{\rm line}$ as defined in Peterson et al. (2004), from the model sum of all broad-line Gaussian components. Following the method described in Wang et al. (2019), a [$-2.5\times$MAD, 2.5$\times$MAD] window is used for the calculation of $\sigma_{\rm line}$ to mitigate the adverse effects of noise and blending in the line wings. Here MAD is the Median Absolute Deviation, calculated by Equation (1) below:

where MED is the median wavelength, defined as the location dividing the line in equal flux on both sides. As discussed in Wang et al. (2019), this latter approach produces $\sigma_{\rm line}$ values that are reproducible, less affected by blending issues and noisy line wings, and is adaptive to the actual line width.

We employ a Monte Carlo approach to estimate uncertainties in the line width measurements (e.g., Shen et al., 2008, 2011). We perturb the original spectrum at each pixel by a random value drawn from a zero-mean Gaussian distribution whose $\sigma$ is set to the flux density uncertainty at that pixel and obtain a mock spectrum. We apply the same fitting procedure to the mock spectrum to obtain spectral measurements. We generate 30 mock spectra per epoch of each object, and estimate the measurement uncertainty of a spectral quantity as the semi-amplitude of the range enclosing the 16th and 84th percentiles of the distribution from the 30 trials.

In addition to the broad-line widths, we also measure the broad-line flux from our spectral fitting. We compare our flux measurements with those reported in G17, H20, and G19, which are based on PrepSpec outputs, and found reasonable agreement (Figure 3).

In our following analysis we use the merged $g$-band light curves as our baseline continuum flux, taken from G17²²2The units of the light curves in G17 are erg s^-1 cm^-2 Å^-1 and the flux scales are slightly adjusted in the merging process (see G17 for details). We first calculate the $g$-band flux of each spectroscopic epoch by convolving the PrepSpec-processed spectrum with the SDSS $g$-band filter curve (Doi et al., 2010) using the equation $\int{f_{\lambda}Rd\lambda}/\int{Rd\lambda}$, where $R$ is the response curve and $f_{\lambda}$ is the flux density. Next, we scale the G17 $g$-band synthetic light curves to have the same mean and standard deviation as ours. We found good matches between the two sets of light curves. Finally we perform the same re-scaling procedure on the Bok and CFHT $g$-band light curves., H20 and G19. These continuum fluxes are not yet host subtracted. Park et al. (2012) found that in Mrk 40 (Arp 151), it can result in much steeper slopes than the virial expectation if the host contribution is not removed. The reason is that the extra contribution from host in the continuum flux will suppress the dynamic range in luminosity variations, making the slope of the line width – continuum variability relation in log-log space steeper than the intrinsic value.

To remove the contribution of host contamination in the continuum light curves, we use the host fraction measured in Yue et al. (2018) (data obtained by private communication)³³3We have tested alternative estimation of the host fraction from spectral decomposition performed in Shen et al. (2015b). We found that only a few individual objects show noticeable differences, and the median values of the measured breathing effects remain more or less the same. Therefore we conclude our statistical results are insensitive to the detailed scheme of host correction in the continuum light curves.. Yue et al. (2018) decomposed the host and nuclear light using CFHT $g$ and $i$ band coadded images and measured the host fraction for 103 SDSS-RM quasars at $z<0.8$. The host fraction is calculated within the 2″ diameter BOSS fiber and could be directly applied to spectroscopy. As listed in Table LABEL:tab:sampleinformation, the $g$-band host fraction $f_{{\rm host},g}$ is moderate or low for most of our objects, but could reach as high as 0.9 (RM772 at $z=0.249$) for low redshift objects. For objects at $z<0.8$, we measured the $g$-band host flux from the mean spectrum. We convolve the mean spectrum with the $g$-band response curve (Doi et al., 2010) to derive the total flux, which is multiplied by the $g$-band host fraction to derive the host flux. The obtained $g$-band host flux is subtracted from the light curves as a constant offset. For objects at $z>0.8$, we assume the host fraction is negligible in $g$-band.

We define a line S/N to restrict the breathing analysis to the subset of quasars with reliable spectral measurements. The line S/N is calculated from the continuum model subtracted residual spectrum and defined as the ratio between the broad-line model flux and the flux uncertainty. The broad-line model flux is calculated within the velocity range of [$-2.5\times$ MAD, $2.5\times$MAD] (same as the range for $\sigma_{\rm line}$ calculation). This is a better line S/N criterion than using the continuum-unsubtracted spectral S/N since we could lose some objects with large equivalent widths (EWs) but low continuum levels.

We first fit all epochs of spectra for all objects in the parent sample, and calculate the median line S/N from all epochs. We select our final sample by requiring the median line S/N$>$2 for the G17 sample, and $>$3 for G19 and H20 samples. The more stringent criterion for the G19 and H20 samples is because the multi-year observations in these two papers have a wider range of S/N variation than that for the single-season observations in G17.

In addition to the restriction on the median line S/N, we also exclude individual epochs with line S/N less than 2 in our following analysis. Table 2 summarizes the statistics of different lines covered by our final sample, and Figure 4 displays their redshift and luminosity distribution.

We now examine the three major broad lines accessible in optical spectroscopy for low-redshift quasars, H$\alpha$, H$\beta$, and Mg ii. For the normal breathing effect we are seeking anti-correlations between variations of the continuum flux and line width. If the broad line flux responds to continuum flux variations, we can use the broad-line flux as a surrogate to track luminosity variations. The advantage of using the broad-line flux is that there is no time delay in the anti-correlation between line flux and width variations. However, the line flux light curve is much less sampled than the continuum light curve and generally has worse S/N. Therefore, we chose to use the better sampled continuum light curve as the reference light curve, and take into account the time lag between the continuum variability and broad-line response.

Figure 5 presents two examples of time variations of continuum flux (host subtracted, if applicable), line flux, and broad-line widths of H$\alpha$, H$\beta$, and Mg ii. There is clear evidence of anti-correlated variations between the continuum flux $f_{\rm con}$ and the line widths (both FWHM and $\sigma_{\rm line}$) for all three lines in these two quasars, i.e., the normal breathing effect.

Next, following Shen (2013), we study the changes in line width as a function of changes in continuum flux in log-log space using the following equation:

To synchronize each continuum light curve and the time series of line width variations, we interpolate the continuum light curve using JAVELIN (Zu et al., 2011). JAVELIN models the continuum variability of quasars assuming a Damped Random Walk (DRW) model (e.g., Kelly et al., 2009; Kozłowski et al., 2010), which provides empirical and more accurate interpolations of the light curves than using simple linear interpolations. The interpolated continuum light curves are plotted in the top panel of Figure 5. For the uncertainties of the interpolated flux, we take into account both the uncertainties reported in the JAVELIN interpolated light curves and the uncertainties of the lag measurements. For the interpolated flux uncertainty from lag measurements, we assume the true lag is uniformly distributed between $[\tau-\tau_{nerr},\tau+\tau_{perr}]$, where $\tau$, $\tau_{nerr}$ and $\tau_{perr}$ are the measured lag, lower and upper uncertainties of the lag, respectively. We generate 50 simulated lags within the range above and obtain the interpolated flux and flux uncertainty for each simulated lag. The final interpolated flux and uncertainty are calculated as follows:

where $w_{i}$ is the weight calculated by $\frac{1}{f_{{\rm err},i}^{2}}$. If an interpolated flux point is located in a seasonal gap, it will have larger uncertainties as showed in Figure 5.

We then pair the synchronized continuum flux and line widths at two different epochs. For each pair of epochs with synchronized line width and continuum flux, we calculate two symmetric sets of differences: $[\log f_{\rm con,1}-\log f_{\rm con,2}$, $\log W_{1}-\log W_{2}]$ and $[\log f_{\rm con,2}-\log f_{\rm con,1}$, $\log W_{2}-\log W_{1}]$, where subscripts 1 and 2 denote the two epochs. For an object with $N$ pairs, we have $N(N-1)$ sets of $[\Delta(\log f_{\rm con})$, $\Delta(\log W)]$. This redundancy of data points enforces a symmetric distribution of the flux and line width changes and hence a zero intercept $\beta$. Therefore, the slope $\alpha$ measures the breathing effect⁴⁴4If we do not use duplicated points in the regression, the slope $\alpha$ does not change. The slope uncertainty increases by $\sim 40\%$.. The uncertainty of $\Delta(\log f_{\rm con})$ and $\Delta(\log W)$ are derived using error propagation from the uncertainties in $f_{\rm con}$ and $W$, respectively. As in Shen (2013), most points are located around zero and these small changes in continuum flux and line widths are consistent with measurement uncertainties. We remove those points with absolute value of $\Delta(\log f_{\rm con})$ less than 1.5 times their uncertainties $\Delta(\log f_{\rm con,err})$. Our selection criterion of SNR${}_{\rm Var,con}>$2 almost guarantees that there are enough points left for slope measurements. There is only one object, RM392 in the H20 sample, that does not have enough points to meaningfully measure the slope after light curve interpolation. We thus exclude this object from further analysis.

To measure the slope $\alpha$ between $\Delta(\log W)$ and $\Delta(\log f_{\rm con})$ we use the Bayesian linear regression package linmix from Kelly (2007). We also measures the Pearson correlation coefficient $r$ between the variations in line width and continuum flux. The final results are summarized in Table 3.

We present the $\Delta(\log W)-\Delta(\log f_{\rm con})$ correlation of the three low-ionization lines in Figure 6 and 7 for the two examples shown in Figure 5. In RM840, the slopes (based on both FWHM and $\sigma_{\rm line}$) of H$\beta$ are consistent with the expected value of $-0.25$ from the virial relation within 3$\sigma$ uncertainties, but the slopes of H$\alpha$ are shallower. In RM303, the slope based on Mg ii FWHM is closer to the expected value of $-0.25$ than $\sigma_{\rm line}$.

We report the values of $\alpha_{F}$ and $\alpha_{S}$ for H$\beta$, H$\alpha$, and Mg ii, measured for individual objects in our sample, in Table 3. Figure 8 displays the distributions of the two slopes for the three lines. Median (MED) slopes and their uncertainties are estimated by bootstrap resampling. We adopt the median of the bootstrap distribution and the semi-amplitude of the range enclosing the 16th and 84th percentiles as the measured MEDs and their uncertainties, respectively. In general, the Lag0 sample and the full sample (Lag0 plus Lag1 sample) show consistent MED slopes, suggesting that using the lag from an alternative line to synchronize the continuum light curve with line width changes is a reasonable approximation. We also find that the fractions of the three cases, i.e., breathing, non-detection of breathing, and anti-breathing, are similar in the Lag0 sample and the full sample for H$\beta$. For Mg ii and H$\alpha$, the fraction of non-detections in the full sample is larger than that in the Lag0 sample. This is because objects in the Lag1 sample generally have smaller flux variation, and H$\beta$ is usually more variable than H$\alpha$ and Mg ii. The statistics of the slope distribution are summarized in Table 3.4. The median (MED) values of $\alpha_{F}$ and $\alpha_{S}$ for H$\beta$ are closer to the expected value of $-0.25$ than for H$\alpha$ and Mg ii. On the other hand, for H$\beta$, the slope based on $\sigma_{\rm line}$ is closer to the expected value of $-0.25$ than the slope based on FWHM.

While the median slope is negative, which is indicative of normal breathing to some extent, the distribution of the slope in our sample is broad for all three low-ionization lines, and some objects show a positive slope indicative of very different breathing behaviors. To quantify the different breathing behaviors, we define three breathing modes: normal breathing, no breathing and anti-breathing:

where $\alpha$ is either $\alpha_{F}$ or $\alpha_{S}$, and $\alpha_{\rm err}$ is the slope uncertainty reported by linmix. We calculate the fractions of normal breathing, no breathing, and anti-breathing, $q_{-}$, $q_{0}$ and $q_{+}$, and summarize the results for all three lines in Table 3.4. H$\beta$ is the line that is most consistent with normal breathing. On the contrary, both Mg ii and H$\alpha$ show a large fraction of no-breathing, consistent with earlier studies (e.g., Shen, 2013). Anti-breathing is rare for all three low-ionization lines. We discuss the implications of these results further in §4.

In our fiducial approach we fixed the Fe ii parameters in the multi-epoch spectral decomposition. The weak Fe ii variability may have a non-negligible effect on the wings of the Mg ii profile, which $\sigma_{\rm line}$ is more sensitive to. To test this effect, we allow the Fe ii parameters to vary and refit objects included in the Mg ii Lag0 sample. The median slopes in this case are $-0.11\pm 0.08$ and $-0.02\pm 0.04$, using FWHM and $\sigma_{\rm line}$, respectively, which are fully consistent with those derived in the fiducial case. Therefore our results are insensitive to the details of treating the Fe ii emission around Mg ii.

Our sample includes a hypervariable quasar SDSS J141324${\rm+}$530527 (RM017). It is identified as a changing-look quasar in Wang et al. (2018) using spectroscopic data from the SDSS-RM project. Dexter et al. (2019) studies the response of different emission lines in terms of line flux and width to the dramatic continuum flux variation using monitoring data from 2009 to 2018. In our work, we only include the light curve in the 2014 observing season for this particular object. In 2014, RM017 showed moderate variations in the continuum flux and line width. The derived $\alpha_{S}$ for H$\alpha$, H$\beta$, and Mg ii in this work are $-0.16\pm 0.03$, $-0.22\pm 0.02$ and $-0.07\pm 0.01$, which are smaller but generally consistent with the values measured in Dexter et al. (2019): $-0.21\pm 0.02$, $-0.26\pm 0.03$ and $-0.12\pm 0.04$, respectively.

Finally, as mentioned in §2, using a more strict criterion for our sample, e.g., ${\rm SNR}_{\rm Var,con}>3$, the distributions of breathing slopes are generally consistent. We present the results derived with the samples defined by the new criterion in Figure 19. While the sample size of H$\alpha$, H$\beta$, and Mg ii Lag0 sample decreases to 8, 17, 9 (compared to 14, 25, 13 for ${\rm SNR}_{\rm Var,con}>2$), the median slopes are generally consistent within uncertainties with our fiducial results.