Variation of Ionizing Continuum: The Main Driver of Broad Absorption Line Variability

Quasar outflows have been considered as a primary feedback source to explain the co-evolution of the central supermassive black holes (SMBHs) and their host galaxies. Theoretical studies and simulations have suggested that kinetic power of the order of only  1% of the Eddington luminosity is sufficient for significant feedback effects on the host galaxy (e.g., Scannapieco & Oh 2004; Di Matteo et al. 2005; Hopkins et al. 2006; Moll et al. 2007; Hopkins & Elvis 2010). Furthermore, in the past 10 years, UV and X-ray spectroscopic studies have revealed powerful high-velocity outflows carrying kinetic power of at least a few percent of bolometric luminosity in some bright active galactic nuclei (AGNs) and quasars (e.g., Tombesi et al. 2010, 2011; Borguet et al. 2012). Finally, massive molecular outflows on the galactic scale were detected in infrared luminous quasars, which are thought to occur in the early phase of nuclear activity, with velocities correlated with the AGN luminosity (e.g., Sturm et al. 2011).

Outflows may manifest themselves as blueshifted absorption lines when they intercept the line of sight to the continuum source. Broad absorption lines (BALs) up to velocities of 0.2 c from the correspondent emission lines are seen in about 10%–40% of quasar populations (e.g., Gibson et al. 2009; Allen et al. 2011). Technically, they are defined as absorption troughs with velocity widths >2000 $\mathrm{km}\,{{\rm{s}}}^{-1}$ at depths >10% below the continuum (Weymann et al. 1991). Most BAL quasars display absorption lines of only high-ionization species, such as N v λ1240, C iv λ1549, and Si iv λ1399, known as HiBAL (Weymann et al. 1991; Gibson et al. 2008; Filiz Ak et al. 2013), and a small percent also feature low ionization species such as Mg ii λ2798 and Al iii λ1857, known as LoBAL (Weymann et al. 1991; Zhang et al. 2010). However, the mass outflow rate and the kinetic power associated with BAL outflows are poorly determined due to the lack of density diagnostics or the distance to the central black hole.

BAL troughs are variable on timescales from several days to years, and the fraction of variable BALs increases with increasing observing intervals (Gibson et al. 2008; Capellupo et al. 2011; Wildy et al. 2014). Two mechanisms are invoked to explain the variability of BAL troughs: (1) changes in the ionization of gas; and (2) absorbing gas moving in and out of the line of sight. In the first case, the alteration of gas ionization can be caused by a varying incident ionizing continuum (case 1a) or the evolving gas density (case 1b) due to, e.g., thermal instability (Waters et al. 2017). In either case, the variability of BAL troughs can either provide important clues regarding the origins or crucial diagnostics for the physical conditions of outflows (Lundgren et al. 2007; Gibson et al. 2008, 2010; Capellupo et al. 2012, 2013; Filiz Ak et al. 2012, 2013; He et al. 2014, 2015; Welling et al. 2014; Wang et al. 2015). For example, in case (1a), the variability timescale will set an upper limit on the recombination timescale, which in turn gives a lower limit on the gas density. In the second case, assuming a Keplerian dominated transverse motion for the absorber, we can constrain the distance of the absorber to the central black holes using the variability timescale of the BAL troughs (Moe et al. 2009; Capellupo et al. 2011). In order to use these constraints, the first step is to establish which mechanism is responsible for the BAL variability.

These mechanisms predict very different relations between continuum and absorption line variabilities. In case (1a), the variability of continuum and absorption lines are highly correlated; in the other two cases, they are independent (see Waters et al. 2017 for the thermal instability case). Furthermore, along the line sight, different parts of outflows see the same ionizing continuum, and thus different absorption components also vary coordinately in the case (1a), while in the other two cases, we do not expect such coherent variations. Previous studies showed that the same BAL troughs of C iv and Si iv BALs or different BAL components of C iv in one object vary in a coordinated manner (strengthening or weakening; e.g., Filiz Ak et al. 2012, 2013; Wildy et al. 2014), indicating that ionization change due to a varying ionizing continuum plays an important role (Hamann et al. 2011; Wang et al. 2015). It is inconclusive whether variations of BAL equivalent width (EW) correlate with continuum flux in individual quasars (as found in, e.g., Trevese et al. 2013; Capellupo et al. 2014). Wang et al. (2015, hereafter Paper I) presented a statistical study of ensemble correlation between BALs and the continuum for a large sample of variable BAL quasars with a multi-epoch spectroscopic data set from SDSS DR10. They found that the EWs of the BALs decrease (increase) statistically when the continuum brightens (dims) for about 70% of BALs, and concluded that a varying ionizing continuum drove BAL variability in most of these quasars. They also showed that the emergence or disappearance of BALs is also connected to the continuum variability.

However, a number of issues were not addressed in that paper. First, coordinated variations of emission and continuum are also expected in the scenario of clouds moving in and out of the picture, if BAL outflows carry dust, as suggested by some studies (He et al. 2014; Leighly et al. 2015; Zhang et al. 2015). Second, the relation between BAL and continuum variability is not monotonic according to photoionization models (Figure 3 in Section 3). Consequently, it lowers the ensemble correlation between BAL and continuum variability. A quantitative analysis should take this into account.

In this paper, we will examine the role that variations of the continuum play in the BAL variabilities and give a quantitative estimation for the fraction of BAL variability driven by the photoionization. The paper is organized as follows. We describe the sample and data analysis in Section 2. An indicator of the ionization state is described in Section 3. The results are presented in Section 4 and Section 5. The conclusions are given in Section 6.

We merged the BAL quasar in the catalog of SDSS data release 7 (DR7, Shen et al. 2011)³ with that of DR12 by Pâris et al. (2016).⁴ Duplicated entries were removed. We compared this catalog with the SDSS spectroscopic catalog and selected quasars with multi-spectroscopic observations. To investigate the variability of C iv and Si iv, we adopt a redshift cut 1.9 < z < 4.7. To ensure detection of major absorption lines, we only keep quasars with at least one spectrum with a signal to noise (S/N) > 10. After these cuts, we obtain a sample of 2005 BAL quasars with two spectra or more. For one quasar with m $(m\geqslant 2)$ spectra, there are ${C}_{m}^{2}=m(m-1)/2$ spectra pairs. As a result, there are total 9918 spectrum pairs in the sample of 2005 BAL quasars. In the following subsection we describe the method of construction of a variable absorption line quasar sample.

2.1. Spectra Flux Calibration

2.2. Using Unabsorbed Quasar Templates to Fit the Spectra

In order to determine the wavelength span of the C iv and Si iv BAL trough, we fit the spectra using unabsorbed template quasar spectra, following Paper I, that are derived from 38,377 non-BAL quasars with $1.5\lt z\lt 4.0$ and ${\rm{S}}/{{\rm{N}}}_{1350}\gt 10$ in SDSS Data Release 7 (DR7). Following Paper I, we fit these templates to the spectra using a double power-law function as a scale factor,

$$\begin{eqnarray}&&{S}_{\lambda }=A[1]{\left(\displaystyle \frac{\lambda }{2000\mathring{\rm A} }\right)}^{A[2]}+A[3]{\left(\displaystyle \frac{\lambda }{2000\mathring{\rm A} }\right)}^{A[4]},\end{eqnarray} \tag{ 1 }$$

where coefficients (A[1], A[3]) and exponents(A[2], A[4]) are determined by minimizing ${\chi }^{2}$. This fitting spans a wavelength range of 1150 to 2850 Å and we add a additional Gaussian component at the Lyα, N v, Si iv, and C iv locations to the unabsorbed templates to improve the fits therein. We iteratively mask spectral pixels lower than the model with a significance of $\geqslant 3\sigma$ on the blue side of Lyα, N v, Si iv, and C iv to exclude possible absorption lines. We select the best template according to the criterion of the maximum number of pixels within 1σ error as the final fitting. The fitting procedure is demonstrated in Figure 1.

Zoom In Zoom Out Reset image size

Like Paper I, we search contiguous deficient pixels for intrinsic absorption line regions of C iv in the normalized spectrum. We screen the moderate to BALs, which are more likely intrinsic; and we screen the deficiency over a width of ${\rm{\Delta }}\mathrm{ln}\lambda \geqslant {10}^{-3}$, or 300 km s⁻¹ in velocity, statistically more significant than 5σ. Then we check by eye to exclude the false ones, caused by an improper fit in most cases. The absorption line regions of Si iv are defined according to the corresponding velocity of C iv.

2.3. Identification of Variable Absorption Lines

As described in Paper I, in order to identify the variation region of C iv BAL, we select the higher S/N spectrum of the pair of quasars as a reference, and then rescale it using the double power-law function (Equation (1)) to account for potential variations of the continuum shape. As shown below, with this recipe we can fit the observations very well. Examples of the fit are displayed in Figure 1. To account for variations of emission line EW, we add/subtract a Gaussian to/from the rescaled spectrum. We restrict the center of the Gaussian to lie within ±500 km s⁻¹ of the line center at the source rest-frame. The added/subtracted Gaussian will be used as an indicator of the increase or decrease of the C iv emission line EW. With this fit, we found that the scatter in the difference spectrum over emission line regions is similar to that of continuum regions. Compared with unabsorbed quasar template matching, the rescaled reference matching usually produces a better fit outside the absorption line region. As a result, we will measure the absorption line variability from the difference spectrum rather than by comparing the EWs obtained in the template fitting.

To identify variable absorption line components in the difference spectrum, we searched for contiguous negative and positive bins to determine the region of variable components (Paper I). We mark all pixels where the difference is larger than 5% of the average value and more than 3σ. Adjacent marked pixels are then connected to form a variable region. Then we expand such regions into neighboring pixels that have the same sign in the difference spectrum but lie at the less than 3σ significant level. After that, we merge the neighboring regions with the same variable sign and with a separation of less than four pixels. Then it will generate the variable region of BAL.

After the visual check, a sample containing 2324 spectrum pairs from 843 quasars, with detected variable absorption lines (C iv λ1549) and variations of continuum (1500 Å) beyond 5%, is identified (Figure 2).

Zoom In Zoom Out Reset image size

We run a series of photoionization simulations to see how ionic column densities respond to variations of ionizing continuum using version c13.03 of Cloudy, described by Ferland et al. (2013). We consider a typical gas density of 10⁶ cm⁻³, since gas ionization is not sensitive to density at a given ionization parameter ($U={Q}_{{\rm{H}}}/[4\pi {r}^{2}{n}_{{\rm{H}}}c]$, where ${Q}_{{\rm{H}}}$ is the source emission rate of hydrogen-ionizing photons, r is the distance to the absorber from the source, c is the speed of light, and ${n}_{{\rm{H}}}$ is the hydrogen number density). As demonstrated in Figure 3, the C iv and Si iv column densities increase first, then reach a peak, and decrease as the ionization parameter increases. Due to the lower ionization potential, Si iv will reach the peak earlier than C iv as the ionization parameter increases. Assuming a Gaussian distribution of optical depth, the EWs of Si iv and C iv BALs can be deduced. The EW ratio of Si iv to C iv BALs,

$$\begin{eqnarray}&&{R}_{\mathrm{EW}}=\displaystyle \frac{{\mathrm{EW}}_{\mathrm{Si}{\rm{}}\mathrm{iv}}}{{\mathrm{EW}}_{{\rm{C}}{\rm{}}\mathrm{iv}}},\end{eqnarray} \tag{ 2 }$$

decreases as the ionization parameter increases at a given ${N}_{{\rm{H}}}$. Generally speaking, ${C}^{+3}$ column density increases (decreases) if there is a (or no) ${C}^{+3}$ ionization front, while R decreases with incident continuum. The R values should have significant statistical differences between the high and low ionization states.

Zoom In Zoom Out Reset image size

For each of the spectrum pairs, we take the mean value of the two observations $R=({R}_{1}+{R}_{2})/2$ as the indicator of the ionization state.

As demonstrated in Figure 3, the ionic column density of a specific species may correspond to a continuum variation, negatively (the EW of the BAL decreases when the continuum brightens) or positively, depending on the ionization of absorbing gas. We divide the sample into two groups according to the pattern of the variation of BAL with respect to that of the continuum: the EW of the BAL decreases (increases) when the continuum brightens (dims) as group T1; and the variation of the EW and continuum in the opposite relation of group T2. As a result, the distributions of R for T1 and T2 should be different as long as the BAL variabilities are driven by changes of the ionizing continuum. Furthermore, the mean value of R distribution for T1 should be smaller than that for T2 in statistics.

4.1. Using the Intrinsic Baldwin Effect

In Paper I, we used the variation of emission line EW as an independent check for the continuum variation. 72.5 ± 2.2% of EWs of the C iv emission lines decrease (increase) when the continuum brightens (dims), which is 10.2σ from no preference (50%). The intrinsic Baldwin effect (e.g., Kinney et al. 1990) is detected at a higher confidence level (23.2σ) in this paper. Thus, we first use the variation of emission lines instead of the continuum variation to check the distributions of R for T1 and T2. We check both the mean value (Student's T Test, T) and distributions of R (Kolmogorov–Smirnov test, KS). As demonstrated in Figure 4, the mean value of R distribution for T2 is significantly (${P}_{{\rm{T}}}=1.7E-02$) larger than that for T1. And the distributions of R for T1 and T2 are significantly different (${P}_{\mathrm{KS}}=3.9E-06$). It is the same as the photoionization model predictions. Like Paper I, about 75% of the spectrum pairs are at the high ionization state.

Zoom In Zoom Out Reset image size

4.2. Sub-sample of SDSS and BOSS

As shown in Figure 2, there are a total of 67 spectrum pairs from SDSS/only (MJD ≤ 55176) and 1177 pairs from BOSS/only (MJD > 55176). Following Paper I, we use a variability amplitude >5% of continuum flux as a threshold. It is shown in Figure 5 that the mean value of R distribution for T2 is significantly larger than that for T1 (${\rm{\Delta }}R$ = 0.164, ${P}_{{\rm{T}}}\,=1.2E-02$ for SDSS; ${\rm{\Delta }}R$ = 0.033, ${P}_{{\rm{T}}}=2.6E-04$ for BOSS). The distributions of R for T1 and T2 are significantly different for both SDSS (${P}_{\mathrm{KS}}=2.5E-03$) and BOSS (${P}_{\mathrm{KS}}=1.1E-12$). It is also the same as that predicted by the photoionization model. The ${\rm{\Delta }}R$ of the SDSS is larger than that of BOSS. The reason may be that the calibration uncertainties of BOSS are larger than those of SDSS.

Zoom In Zoom Out Reset image size

4.3. The Whole Sample

There are a total of 2324 spectrum pairs from 843 quasars with detected variable absorption lines and variations of continuum beyond 5%. We gradually increase the amplitude of variations of the continuum to compare the distributions of R for T1 and T2. As shown in Figure 6, the distributions of R for T1 and T2 are significantly (${P}_{\mathrm{KS}}=5.4E-19$) different. The mean value of R distribution for T2 is significantly (${P}_{{\rm{T}}}=7.5E-07$) larger than that for T1 and ${\rm{\Delta }}R$ increases with the amplitude of variations of continuum.

Zoom In Zoom Out Reset image size

The above distributions of R for T1 and T2 are all significantly different. And the mean value of R distribution for T2 is significantly larger than that for T1. This agrees with the prediction of the photoionization model and this evidence supports the idea that variation of the ionizing continuum is the main driver of BAL variability.

4.4. Estimate of the Fraction of BAL Variability Driven by the Photoionization

On one hand, the uncertainty of flux calibration does not affect the measurement of BAL EW. On the other hand, the variation of the covering factor of the BAL gas does not change the UV continuum. As a result, considering the uncertainty of flux calibration or the variation of covering factor of the outflow, part of the BALs will be randomly assigned to the two types (T1 and T2). The BALs that originally belonged to T1 (T2) may be assigned to T2 (T1) at the same probability (P_mis in Figure 7) for the two types. As shown in the top panel of Figure 7, the fraction of observed T1 is A. The fraction of misassignment from T1 to T2 is n, while the fraction from T2 to T1 is m. Then the misassignment possibility is ${P}_{\mathrm{mis}}=n/(A-m+n)=m/(1-A-n+m)$. Then, the fraction of BAL variability driven by the photoionization is ${F}_{\mathrm{pho}}=1-2(n+m)$.

Zoom In Zoom Out Reset image size

We use the skew Gaussian function,

$$\begin{eqnarray}&&N(r)={a}_{1}\exp \left[\displaystyle \frac{-{(r-{a}_{2})}^{2}}{2{a}_{3}^{2}}\right](1+\mathrm{erf}[{a}_{4}(r-{a}_{2})]),\end{eqnarray} \tag{ 3 }$$

to describe the intrinsic distributions of T1 and T2 of our sample, where $\mathrm{erf}(r)={\int }_{-r}^{r}{e}^{-{t}^{2}}{dt}$ is the error function and the parameter a₄ is the coefficient of skewness. The observed T1 consists of two components, the blue dotted line and the red shadow region, while T2 consists of the red dotted line and the blue shadow region. The blue (red) shadow region is the fraction of misassignment from T1 (T2) to T2 (T1). We fit the distributions combining subsamples at different thresholds of continuum variability. During the fitting, we tie all the parameters together except a₁, which varies independently for subsamples at different thresholds of continuum variability, with the minimizing ${\chi }^{2}=2.29$. The error of each bin of the distributions is estimated using the bootstrap method.

This shows that F_pho increases with thresholds of continuum variability. This is expected because the fraction of misassignment caused by the uncertainty of flux calibration gradually decreases with the amplitude of continuum variability. F_pho is about 80% when the threshold of continuum variability is beyond 30%. As shown in Figure 8,  R is not thoroughly monotonic with ionization parameter and there exists a plateau when ${N}_{{\rm{H}}}$ is large. The plateau elongates as ${N}_{{\rm{H}}}$ increases. The intrinsic distribution of T1 is broader than T2 from the fitting results. This may imply that many T2 quasars are at the plateau.

Zoom In Zoom Out Reset image size

As mentioned in Paper I, the high concordance between continuum and absorption line variations requires a recombination time ${t}_{\mathrm{rec}}={(\alpha {n}_{e})}^{-1}$ shorter than both the timescale of typical continuum variations and the interval between the two observations (e.g., Barlow et al. 1992), where α is the recombination rate and n_e is the electron density. As a result, if the timescale of typical continuum variations or the interval between the two epochs of the spectrum pair is shorter than the recombination time, there will be no significant difference between the distributions of R for T1 and T2. We gradually change the time interval between the two epochs and the amplitude of variations of continuum. As shown in Figure 9, it still has a significant difference (${P}_{\mathrm{KS}}=0.031$, ${P}_{{\rm{T}}}=0.032$) between the distributions of T1 and T2 when ${\rm{\Delta }}T\lt 6\,\mathrm{days}$ and the amplitude of variations of continuum is larger than 15%. This result suggests that there exist BAL outflows whose recombination timescales are only about a few days. Note that ${\rm{\Delta }}R$ clearly increases (0.019–0.134) with the amplitude of variations of continuum (5%–25%) when ${\rm{\Delta }}T\lt 5\,\mathrm{days}$. This implies that the upper limit of the recombination timescale may be shorter than $5\,\mathrm{days}$. This agrees with the results of Grier et al. (2015), who found the shortest BAL variability timescale (1.2 days) yet reported. They concluded that this result is most likely a rapid response to changes in the ionizing continuum.

Zoom In Zoom Out Reset image size

For low-redshift AGNs, the most direct method to determine the outflow radius is to use spatially resolved spectroscopy, especially integral field unit (IFU) spectroscopy, (e.g., Barbosa et al. 2009; Riffel & Storchi-Bergmann 2011; Liu et al. 2013a, 2013b, 2014, 2015; Rupke & Veilleux 2013), but for high-redshift ($z\gt 2$) BAL quasars, the realistic approach is to derive the galactocentric distance of the outflow by the hydrogen number density ${n}_{{\rm{H}}}$, which can be determined from the absorption lines of the excited states of ions (e.g., Fe ii*, Si ii*, S iv*). During the last decade or so, the outflow radii have been measured for a number of individual quasars using density-sensitive absorption lines from excited levels (Arav et al. 1999, 2008, 2013, 2015; de Kool et al. 2001, 2002a, 2002b; Hamann et al. 2001; Moe et al. 2009; Bautista et al. 2010; Dunn et al. 2010; Aoki et al. 2011; Borguet et al. 2012; Chamberlain et al. 2015). Due to line blending, the method can only apply to relatively narrow lines. Complementary to these conventional approaches, we can monitor some quasars continuously whose BAL variabilities and recombination timescales are only a few days. Then we can constrain the electron density from the time lags of BAL to continuum variability.

We present a statistical analysis of the variability of BAL for a sample of 2324 spectrum pairs in 843 quasars (1.9 < z < 4.7) with detected variable absorption lines and variations of continuum beyond 5% from SDSS DR12. The main results are as follows.

• (1)  
  Our study shows that the distributions of R for T1 (high ionization state) and T2 (low ionization state) are significantly (${P}_{\mathrm{KS}}=5.4E-19$) different, while the mean value of R of T1 is significantly (${P}_{{\rm{T}}}=7.5E-07$) lower than that of T2. These results suggest that the variabilities of BALs are dominated by the variation of ionizing continuum.
• (2)  
  Considering the uncertainty of flux calibration and the variation of covering factor of the outflow, part of the BALs will be randomly assigned to the two types. As a result, the significance of difference of R for T1 and T2 will be weakened by the above two factors. We give an estimate for the fraction of BAL variability driven by the variation of ionizing continuum. The fraction of BAL variability driven by the variation of ionizing continuum increases with the thresholds of continuum variability. It is about 80% when the thresholds of continuum variability are beyond 30%.
• (3)  
  There is still a significant difference (${P}_{\mathrm{KS}}=0.031$, ${P}_{{\rm{T}}}=0.032$) between the distributions of T1 and T2 when ${\rm{\Delta }}T\lt 6$ days and the amplitude of variations of continuum is larger than 15%. This result suggests that there exist BAL outflows whose recombination timescales are only about a few days. To determine the electron density, we can monitor some quasars continuously whose BAL variabilities and recombination timescales are short. Once the density is known, the size of the BAL outflow can be measured.

The authors thank the anonymous referee for the constructive comments and suggestions, which improved the clarity of the manuscript. We acknowledge the financial support by the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences (XDB09000000), NSFC (NSFC-11233002, NSFC-11421303, U1431229), the National Basic Research Program of China (grant No. 2015CB857005), and the National Science Foundation of China (nos. 11373024 & 11233003).

Guilin Liu is supported by the National Thousand Young Talents Program of China, and acknowledges the grant from the National Natural Science Foundation of China (No. 11673020 and No. 11421303) and the Ministry of Science and Technology of China (National Key Program for Science and Technology Research and Development, No. 2016YFA0400700).