Variation of ionizing continuum: the main driver of Broad Absorption Line Variability

TL;DR

This study analyzes quasar broad absorption line variability using SDSS data, showing that ionizing continuum changes primarily drive BAL variability, with some outflows responding on very short timescales of days.

Contribution

It provides statistical evidence that ionizing continuum variations are the main driver of BAL variability, supported by a large multi-epoch quasar dataset and photoionization model predictions.

Findings

Over 80% of BAL variability is driven by ionizing continuum changes.

BAL outflows can respond within days, indicating short recombination times.

The ratio of Si iv to C iv BALs correlates with continuum variations as predicted by models.

Abstract

We present a statistical analysis of the variability of broad absorption lines (BALs) in quasars using the large multi-epoch spectroscopic dataset of the Sloan Digital Sky Survey Data Release 12 (SDSS DR12). We divide the sample into two groups according to the pattern of the variation of C iv BAL with respect to that of continuum: the equivalent widths (EW) of the BAL decreases (increases) when the continuum brightens (dims) as group T1; and the variation of EW and continuum in the opposite relation as group T2. We find that T2 has significantly (P_T<10-6 , Students T Test) higher EW ratios (R) of Si iv to C iv BAL than T1. Our result agrees with the prediction of photoionization models that C +3 column density increases (decreases) if there is a (or no) C +3 ionization front while R decreases with the incident continuum. We show that BAL variabilities in at least 80% quasars are driven by the variation of ionizing continuum while other models that predict uncorrelated BAL and continuum variability contribute less than 20%. Considering large uncertainty in the continuum flux calibration, the latter fraction may be much smaller. When the sample is binned into different time interval between the two observations, we find significant difference in the distribution of R between T1 and T2 in all time-bins down to a deltaT < 6 days, suggesting that BAL outflow in a fraction of quasars has a recombination time scale of only a few days.

Equations

$$S_{\lambda }=A[1]{\left(\frac{\lambda }{2000\mathring{A}}\right)}^{A[2]}+A[3]{\left(\frac{\lambda }{2000\mathring{A}}\right)}^{A[4]}$$

$${\mathrm{R}}_{\mathrm{EW}}=\frac{\mathrm{E}{\mathrm{W}}_{\mathrm{Si}\text{sc}\mathrm{iv}}}{\mathrm{E}{\mathrm{W}}_{\mathrm{C}\text{sc}\mathrm{iv}}},$$

$$N(r)=a_1\exp \left[\frac{-(r-a_2{)}^2}{2a_3^2}\right]\left(1+erf[a_4(r-a_2)]\right)$$

Full text

1. Introduction

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 suggested that kinetic power of order of only  1% of the Eddington luminosity is deemed sufficient for significant feedback effects on the host galaxy (e.g., Scannapieco & Oh 2004; Di Matteo , Springel & Hernquist 2005; Hopkins et al. 2006; Moll et al. 2007; Hopkins & Elvis 2010). Furthermore, in the past ten 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., Borguet et al. 2012; Tombesi 2010, 2011). Finally, massive molecular outflows on the galactic scale were detected in infrared luminous quasars, which is thought in the early phase of nuclear activity, with velocities correlated with the AGN luminosity (e.g., Sturm 2011).

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

BAL troughs are variable on time scales 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). So far, 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. 2016). In either case, the variability of BAL troughs can provide either important clues to origin or crucial diagnostics for the physical conditions of outflows (Lundgren et al. 2007; Gibson et al. 2008, 2010; Filiz Ak et al. 2012, 2013; Capellupo et al. 2012, 2013; Welling et al. 2014; He et al. 2014, 2015; Wang et al. 2015). For example, in case 1a), the variability time scale will set an upper limit on recombination time scale, 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 time scale of the BAL troughs (Moe et al. 2009; Capellupo et al. 2011). In order to use these constraints, a first step is to establish which mechanism is responsible for the BAL variability.

These mechanisms predict very different relations between continuum and absortpion line variabilities. In case 1a), variability of continuum and absorption lines are highly correlated; while in other two cases, they are independent (see Waters et al. 2016 for thermal instability case). Furthermore, along the line sight, different part of outflows see the same ionizing continuum, thus different absorption components vary also coordinately in the case 1a), while in 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 varying ionizing continuum plays important role (Hamann et al. 2011; Wang et al. 2015). It is inconclusive whether variations of BAL EW correlates with continuum flux in individual quasars (as found in e.g, Trevese et al. 2013; Capellupo, Hamann & Barlow 2014). Wang et al. (2015, hereafter paper I) presented a statistical study of ensemble correlation between BALs and continuum for a large sample of variable BAL quasars with multi-epoch spectroscopic dataset of the SDSS DR10. They found that the EWs of the BALs decrease (increase) statistically when the continuum brightens (dims) for about 70% BALs, and concluded that varying ionizing continuum drove BAL variability in most of these quasars. They also showed that emergence or disappearance of BALs are also connected to the continuum variability.

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

In this paper, we will examine the role of variation of continuum playing 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 §2. An indicator of the ionization state is described in §3. The results are presented in §4 and §5. The conclusions are given in §6.

2. Variable Absorption Line Quasar Sample

We merged the BAL quasar in the catalog of SDSS data release 7 (DR7, Shen et al. 2011) 111http://das.sdss.org/va/qso_properties_dr7/dr7.htm with that of DR12 by Pâris et al. 2016 222http://dr12.sdss3.org/datamodel/files/BOSS_QSO/DR12Q/ DR12Q_BAL.html. Duplicated entries are 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 signal to noise SNR $>$10. After these cuts, we obtain a sample of 2005 BAL quasars with two spectra or more. For one quasar with m $(m\geq 2)$ spectra, there are $C^{2}_{m}=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

One of the primary scientific goals of Baryon Oscillation Spectroscopic Survey (BOSS) is to search for the baryon acoustic oscillation (BAO) signal in the Lyman $\alpha$ forest. As mentioned in Paper I, the fiber positions of BOSS quasar targets were purposefully offset in order to optimize the throughput of light at 4000Å, while the standard stars used for flux calibration are positioned for 5400Å  (Dawson et al. 2013). As a result, there exists the systematic offest between the spectrophotometric flux calibrations of SDSS and BOSS.

The quasar locus for high redshift quasars overlaps with the stellar locus in SDSS color-color space. As a result, many stars are targeted as quasars until the spectra were taken. These stars were observed in the same manner as the observed quasars (Harris et al. 2016). Thus, they have the same systematic offset with the quasars in BOSS. We give an average correction to the BOSS quasar spectra. The correction is small at the bluest wavelength and up to about 17% at the reddest wavelength (Fig. 13 of Harris et al. 2016).

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 as in Paper I, which are drived from 38,377 non-BAL quasars with $1.5<z<4.0$ and $SNR_{\rm 1350}>10$ in SDSS Data Release 7 (DR7). Like Paper I, we fit these templates to the spectra using a double power-law function as a scale factor,

[TABLE]

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 Lyman $\alpha$, 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 $\geq 3\sigma$ on the blue side of Lyman $\alpha$, N v, Si iv, C iv to exclude possible absorption lines. We select the best template according to the criterion of the maximum number of pixels within 1$\sigma$ error as the final fitting. The fitting procedure is demonstrated in Fig. 1.

Like Paper I, we search contiguous deficient pixels for intrinsic absorption lines region of C iv in the normalized spectrum. We screen the moderate to broad absorption lines, which are more likely intrinsic, deficiency over a width of $\Delta\ln\lambda\geqslant 10^{-3}$, or 300 km s*-1* in velocity, and statistically significant than 5$\sigma$. 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 $SNR$ spectrum of the pair of a quasar as a reference, and then rescale it using the double power-law function (Eq.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 the 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 $\pm$500 km/s of the line center at the source rest frame. The add/subtract Gaussian will be used as an indicator of the increase or decrease of 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 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$\sigma$. Adjacent marked pixels are then connected to form a variable region. Then we expand such regime into neighboring pixels that have the same sign in the difference spectrum but at the less than 3$\sigma$ 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 $\lambda 1549$) and with the variations of continuum (1500Å) beyond 5% is identified (Fig. 2).

3. an indicator of the ionization state

We run a series of photoionization simulations to see how ionic column densities respond to the variations of ionizing continuum using version c13.03 of Cloudy, described by Ferland et al. (2013). We consider a typical gas density of 10^{6}\hbox{cm{}^{-3}}, 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 Fig. 3, the C iv and Si iv column densities increase first, then reach a peak, and decrease after 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:

[TABLE]

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 the R decreases with incident continuum. The R values should have the significant statistical differences between the high and low ionization states.

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.

4. The distributions of R for T1 and T2

As demonstrated in Fig. 3, the ionic column density of a specific species may respond to a continuum variation, negatively (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 continuum: the EW of the BAL decreases (increases) when the continuum brightens (dims) as group T1; and the variation of EW and continuum in the opposite relation as 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. Further more, the mean value of R distribution for T1 should 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$\pm$2.2% of EWs of the C iv emission lines decrease (increase) when the continuum brightens (dims), which is 10.2$\sigma$ from no preference (50%). The intrinsic Baldwin effect (e.g., Kinney et al. 1990) is detected at a higher conficence level (23.2$\sigma$) in this paper. So, we firstly 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 Fig. 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_{\rm 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.

4.2. Sub-sample of SDSS and BOSS

As shown in Fig. 2, there are total 67 spectrum pairs from the SDSS/only (MJD $\leq$ 55176) and 1177 pairs from the BOSS/only (MJD $>$ 55176). Like Paper I, we use the variability amplitude $>$5% of continuum flux as a threshold. It is shown in Fig. 5 that the mean value of R distribution for T2 is significantly larger than that for T1 ($\Delta R$ = 0.164, $P_{\rm T}=1.2E-02$ for SDSS; $\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_{\rm KS}=2.5E-03$) and BOSS ($P_{\rm KS}=1.1E-12$). It is also the same as predicted by the photoionization model. The $\Delta R$ of the SDSS is larger than that of the BOSS. The reason maybe that the calibration uncertainties of BOSS is larger than that of the SDSS.

4.3. The whole sample

There are total 2324 spectrum pairs from 843 quasars with detected variable absorption lines and the variations of continuum beyond 5%. We gradually increase the amplitude of variations of continuum to compare the distributions of R for T1 and T2. As shown in Fig. 6, the distributions of R for T1 and T2 are significantly ($P_{\rm 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 the $\Delta R$ increases with the amplitude of variations of continuum.

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. It agrees with the predication of the photoionization model. These evidences support the idea that variation of the ionizing continuum is the main driver for BAL variability.

4.4. Estimate of the fraction of BAL variability driven by the photoionization

On the one hand, the uncertainty of flux calibration doesn’t affect the measurement of BAL EW. On the other hand, the variation of covering factor of the BAL gas doesn’t 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 which originally belong to T1 (T2), may be assigned to T2 (T1) at the same probability ($P_{mis}$ in Fig. 7) for the two types. As shown in top panle of Fig. 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_{mis}=n/(A-m+n)=m/(1-A-n+m)$. Then, the fraction of BAL variability driven by the photoionization is $F_{pho}=1-2(n+m)$.

We use the skew gaussian function,

[TABLE]

to describe the intrinsic distributions of T1 and T2 of our sample, where $erf(r)=\int_{-r}^{r}e^{-t^{2}}dt$ is the error function and the parameter $a_{4}$ 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 sub-samples at different thresholds of continuum variability. During the fitting, we tie all the parameters together except $a_{1}$ which varies independently for sub-samples 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.

It shows that $F_{pho}$ increases with the 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 thresholds of continuum variability is beyond 30%. As shown in Fig. 8, the R isn’t thoroughly monotonic with ionization parameter and there exists a plateau when the $N_{\rm H}$ is large. The plateau elongates as the $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.

5. The recombination time scale

As mentioned in Paper I, the high concordance between continuum and absorption line variations requires a recombination time $t_{rec}=(\alpha n_{e})^{-1}$ shorter than both the time scale of typical continuum variations and the interval between the two observations (e.g., Barlow et al. 1992), where $\alpha$ is the recombination rate and $n_{e}$ is the electron density . As a result, if the time scale 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 Fig. 9, it still has a significant difference ($P_{\rm KS}=0.031$, $P_{\rm T}=0.032$) between the distributions of T1 and T2 when $\Delta T<6~{}days$ and the amplitude of variations of continuum is larger than 15%. This result suggests that there exists the BAL outflows whose recombination time scales are only about several days. Here note that, $\Delta R$ clearly increases (0.019 to 0.134) with the amplitude of variations of continuum (5% to 25%) when $\Delta T<5~{}days$. It implies that the upper limit of recombination time scale may be shorter than $5~{}days$. It agrees with the result of Grier et al. 2015 which they found the shortest BAL variability time scale (1.2days) yet reported. They concluded that it’s most likely a rapid response to changes in the ionizing continuum.

For low-redshift AGN, the most direct method to determine the outflow radius is using spatially resolved spectroscopy, especially integral field unit (IFU) spectroscopy, e.g. (Barbosa et al. 2009; Riffel & Storchi-Bergmann 2011; Rupke & Veilleux 2013; Liu et al. 2013a,b, 2014, 2015), but for high-redshift ($z>2$) BAL quasars, the realistic approach is deriving 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; de Kool et al. 2001, 2002a, b; Hamann et al. 2001; Arav et al. 2008; Moe et al. 2009; Bautista et al. 2010; Dunn et al. 2010; Aoki et al. 2011; Borguet et al. 2012; Arav et al. 2013, 2015; Chamberlain et al. 2015). Due to line blending, the method can only apply to relatively narrow line. Complementary to these conventional approaches, we can monitor some quasars continuously whose BAL variabilities and recombination time scales are only several days. Then we can constrain on the electron density from the time lags of BAL to continuum variability.

6. Conclusions

In this paper, 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 the variations of continuum beyond 5% from the 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_{\rm KS}=5.4E-19$) different, while the mean value of R of T1 is significantly ($P_{\rm T}=7.5E-07$) lower than the 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 above two aspects. We give an estimate for the fraction of BAL variability driven by the variation of ionizing continuum. The fraction of BALs 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 is beyond 30%.

(3) It still has a significant difference ($P_{\rm KS}=0.031$, $P_{\rm T}=0.032$) between the distributions of T1 and T2 when the $\Delta T<$ 6 days and the amplitude of variations of continuum is larger than 15%. This result suggests that there exists the BAL outflows whose recombination time scales are only about several days. To determine the electron density, we can monitor some quasars continuously whose BAL variabilities and recombination time scales are short. Once the density is known, the size of the BAL outflow can be measured.

7. ACKNOWLEDGMENTS

Authors thank the anonymous referee for constructive comments and suggestions for improving 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), National Basic Research Program of China (grant No. 2015CB857005) and 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).