SDSS J1138+3517: a quasar showing remarkably variable broad absorption lines

Abstract

1 INTRODUCTION

Blueshifted quasar absorption lines, which indicate the presence of large amounts of outflowing material from the central engine, have become an increasingly important topic over the past two decades. This followed the realization of the potential influence these flows can exert on the host galaxy by restricting black hole growth (Springel, Di Matteo & Hernquist 2005; King 2010) and contributing to the evolution of galactic bulges (Magorrian et al. 1998; Silk & Rees 1998). This feedback process is also supported by evidence that outflows could extend hundreds of parsec into the surrounding interstellar medium (Barlow 1994), depositing significant amounts of kinetic energy in the process (Arav et al. 2013).

Absorption features observed in the ultraviolet (UV) portion of the quasar's rest-frame spectrum can be blueshifted by as much as 0.1c, where c is the speed of light in a vacuum, from the corresponding emission line centre. These absorption lines are divided into categories according to their velocity width, the widest being broad absorption lines (BALs) extending over at least 2000 km s⁻¹ where flux is less than 90 per cent of the continuum (Weymann et al. 1991). Quasars hosting at least one BAL are known as broad absorption line quasars (BALQSOs). Absorption can also be seen as narrow absorption lines (NALs), which extend over a few hundred km s⁻¹ and mini-BALs, which are of intermediate width between NALs and BALs (Hamann & Sabra 2004).

The population of BALQSOs, of which SDSS J113831.4+351725.2 (hereafter SDSS J1138+3517) is a member, account for approximately 10 to 20 per cent of the total quasar population (Reichard et al. 2003a; Knigge et al. 2008; Scaringi et al. 2009). However SDSS target selection effects may make this value as high as 41 per cent (Allen et al. 2011). If BAL observability is governed by line-of-sight orientation, then large-scale outflows may occur in the vast majority of quasars, even where no BALs are visible (Schmidt & Hines 1999). Common BAL subtypes include High-Ionization BALQSOs (HiBALs) (∼85 per cent) whose spectra only show broad absorption due to high ionization transitions such as C iv λλ1548,1550 Si iv λλ1394,1403 and N v λλ1239,1243 (Sprayberry & Foltz 1992; Reichard et al. 2003b). Low Ionization BALQSOs (LoBALs), in addition to high-ionization BALs, also exhibit broad absorption resulting from low ionization species such as Mg ii, Al iii and on rare occasions Fe ii (known as FeLoBALs).

Ionizing continuum changes are known to drive broad emission line (BEL) variability (Peterson et al. 1998; Vanden Berk, Wilhite & Kron 2004; Wilhite et al. 2006); however there is as yet no consensus behind the dominant mechanism involved in absorption line variability. Several large statistical studies of BAL variability have been published in recent years (Lundgren et al. 2007; Gibson et al. 2008, 2010; Capellupo et al. 2011, 2012, 2013; Filiz Ak et al. 2013; Wildy, Goad & Allen 2014). These studies refer to coordinated variability across widely separated absorption troughs as evidence favouring changes in input ionizing flux as the main mechanism behind variability (Filiz Ak et al. 2013), while variability over narrow velocity ranges within troughs is attributed to absorbing gas moving across the line of sight (Gibson et al. 2008; Capellupo et al. 2012).

As in Wildy et al. (2014), we use a novel method to reconstruct the un-absorbed spectral profile (Allen et al. 2011) against which the absorption can be measured. Absorption is identified based on this emission profile as it contains the total input spectral flux (continuum+emission line) entering the absorbing gas at each wavelength bin. We use this method to estimate the column density at each epoch for the Si iv and C iv ions. This has advantages over purely continuum-based absorber strength measurements as it allows absorption depth to be accurately identified in spectral regions where the absorption is overlapping with an emission line profile, allowing absorption features to be studied to arbitrarily low (less negative) velocities.

1.1 The unusual BAL behaviour of SDSS J1138+3517

The quasar examined in this paper, SDSS J1138+3517, has a redshift of z = 2.122 (Hewett & Wild 2010) and was investigated as part of the multi-BALQSO variability study undertaken in the rest-frame UV by Wildy et al. (2014). Using a sample of 50 quasars, Si iv λ1400 and C iv λ1549 BALs were identified and statistics gathered on their equivalent width (EW) changes over two epochs separated by quasar rest-frame time-scales ranging from approximately 1 to 4 years. Of the 50 BALQSOs, SDSS J1138+3517 showed the largest absolute change in EW in both Si iv and C iv absorption as well as a large fractional EW change (|ΔEW/〈EW〉|), occurring over a rest-frame time interval of ∼1 year. A third epoch indicated another dramatic variation in the BALs, while a fourth epoch showed relatively little change from the third except for significant variability in one C iv absorber, with these subsequent observations also being made in successive 1-year rest-frame time intervals.

In Filiz Ak et al. (2013), the variability of 428 C iv BALs in 291 quasars was found to produce a ΔEW dependence on rest-frame time interval which is well predicted by a random-walk model (see fig. 29 of that paper). Using the same BAL identification scheme as that study, the variation across the two observational epochs examined in Wildy et al. (2014) of the most variable C iv BAL region in SDSS J1138+3517 (spanning −23 400 km s⁻¹ to −7500 km s⁻¹) is ΔEW = −13.8 ± 2.9 Å. This value lies well outside the random-walk model prediction at the epoch separation time-interval of 360 d in the quasar rest frame, suggesting that this object may exhibit a variability mechanism which is rare in the BALQSO population.

In Filiz Ak et al. (2014), which use the same BAL identification scheme as above, the variability of C iv BALs was investigated over time-scales ranging from ∼1 to a few years, including those in quasars additionally containing Si iv BALs of overlapping velocity. Again examining the most variable C iv BAL in SDSS J1138+3517 identified using their definition, the value of ΔEW, along with a value of ΔEW/〈EW〉 = −1.3 ± 0.3, would place it among the top few most variable C iv BALs in their sample of 454 (includes C iv BALs both unaccompanied by Si iv absorption and overlapping in velocity with a Si iv BAL), highlighting the extreme variability of this source. In addition, an investigation by Gibson et al. (2010) included observations of two quasars over a similar rest-frame time interval of between 0.5 and 2 yr, neither of which showed variability as dramatic as in SDSS J1138+3517.

Answering the question of what drives the strong BAL variability observed in this quasar is the principal objective of this investigation. Observations of the source are detailed in Section 2. In Section 3, we quantify Si iv and C iv ionic column densities using different assumptions regarding the line-of-sight geometry, including models where there is incomplete coverage of the emission region by a homogeneous absorber. This is aided by the use of Gaussian profiles to fit the Si iv absorption components, which allows calculations to be performed based on the changes in these components even where they are not fully un-blended. In Section 4 we compute a grid of photoionization models, allowing examination of the theoretical variation of ionic column density over a range of ionization parameter at given hydrogen number densities and column densities. These grid models are compared to the best results for the ionic column densities from Section 3, allowing us to test the hypothesis of changes in ionizing flux being the driver of BAL variability. In Section 5 we discuss the results from Sections 3 and 4 and their implications for the variability mechanism. We also attempt to constrain the mass outflow rate and hence kinetic luminosity.

This paper assumes a flat Λcold dark matter cosmology with H₀ = 67 km s⁻¹ Mpc⁻¹, Ω_M = 0.315 and Ω_Λ = 0.685, as reported from the latest results of the Planck satellite (Planck Collaboration XVI 2014).

2 OBSERVATIONS

2.1 The Sloan Digital Sky Survey

The epoch 1 observation was obtained from DR6 of the Sloan Digital Sky Survey (SDSS). Beginning in 2000, the SDSS had imaged approximately 10 000 deg² of the sky as of 2010 June (Schneider et al. 2010), using the 2.5-m Telescope at the Apache Point Observatory, New Mexico, USA (York 2000). Imaging is carried out using a CCD camera that operates in drift-scanning mode (Gunn et al. 1998). Five broad-band filters, u,g,r,i and z, normalized to the AB system, are used, covering a wavelength range of 3900 to 9100 Å with an approximate spectral resolution of λ/Δλ = 2000 at 5000 Å (Stoughton et al. 2002). Further information on the SDSS data reduction process is given in Lupton et al. (2001) and Stoughton et al. (2002).

2.2 The William Herschel Telescope

The William Herschel Telescope (WHT) ISIS (Intermediate dispersion Spectrograph and Imaging System) double-armed (red and blue) spectrograph was used for long-slit spectroscopy of the target during semester A 2008, 2011 and 2014 corresponding to epochs 2, 3 and 4, respectively, with epochs 3 and 4 being observed as part of the ING service programme and epoch 2 taken from the sample used by Cottis et al. (2010). The observations provide spectral coverage from 3000 to 10 000 Å. For epoch 2 the 5700 dichroic was used, providing spectra from the two arms with overlap in the range 5400 to 5700 Å. The R158B and R158R gratings, giving a nominal spectral resolution of 1.6 Å per pixel in the blue and 1.8 Å per pixel in the red, were also utilized for epoch 2. Epoch 3 and 4 observations employed the 5300 dichroic, giving an overlap range of 5200–5500 Å, while the R300B and R316R gratings were used for the blue and red observations, providing spectral resolutions of 0.86 Å per pixel and 0.93 Å per pixel, respectively. In the dichroic overlap regions, an error weighted mean was used to combine the spectra from the two arms for each observation. For all WHT observations, a slit width of 1.5 arcsec was chosen to match typical seeing conditions and give reasonable throughput without compromising the spectral resolution. Multiple exposures were taken (in order to remove cosmic ray hits) bracketed by arc-lamp exposures for wavelength calibration. Standard star spectra were taken on the observing nights to allow flux calibration of the targets. Correction of CCD images using bias and flat-field exposures was performed within the irafccdproc routine. Spectra were extracted, along with sky background removal, using the apall task within the iraf longslit package. This used an extraction slit width of 10 pixels and optimal weighting, with 2 pixels per wavelength bin. Wavelength and flux calibration were applied to the extracted spectra within this same package. The spectra were placed on an absolute flux scale using spectra of a photometric standard observed at similar airmass. No attempt was made, such as through the use of grey shifts, to adjust the flux scale on the calibrated spectra towards levels which more closely matched across the epochs. Table 1 summarizes the SDSS J1138+3517 observations.

3 ANALYSIS

3.1 Spectral properties

The observed spectra at each epoch correspond to rest-frame UV wavelengths which include the region containing Si iv and C iv emission and the corresponding blueshifted absorption lines resulting from those same atomic transitions in outflowing gas. Since the spectral resolution differs between spectra, in order to accurately identify changes in absorption between each epoch the SDSS and WHT epoch 3 and epoch 4 spectra were convolved with Gaussians of appropriate full-width at half maxima to approximate the resolution of the WHT epoch 2 observation (using skyline widths as a guideline), followed by resampling of all four spectra on to a wavelength grid of 3.7 Å. This closely matches the bin width of the epoch 2 observation. Small differences in wavelength calibration were removed by aligning narrow emission and absorption features in the WHT spectra with the same lines in the SDSS spectrum. The SDSS and both WHT spectra were corrected for the effects of Milky Way Galactic Extinction using the method of Cardelli, Clayton & Mathis (1989) with A_V values taken from Schlegel, Finkbeiner & Davis (1998).

Using the reconstruction of the unabsorbed SDSS spectrum from Wildy et al. (2014) originally developed in Allen et al. (2011), appropriate reconstructions were also generated for the three WHT observations. This was achieved as follows: (i) a power-law continuum was fitted to each spectrum using line-free continuum bands as outlined in Wildy et al. (2014); (ii) the SDSS continuum was subtracted from the SDSS reconstruction; (iii) the WHT spectra had their respective continua subtracted; (iv) in the continuum-subtracted SDSS reconstruction, the Si iv λ1400, C iv λ1549, O i λ1304 and C ii λ1334 emission line boundaries were identified; (v) the emission lines from part (iv) were scaled to match the red side (unabsorbed) of the emission lines in each of the WHT observations, and (vi) the continuum was added to the scaled emission lines for each of the WHT observations, creating appropriate reconstructions for epochs 2, 3 and 4.

For the SDSS reconstruction used in the subsequent analysis, the procedure was carried out to scale the Si iv emission line, as the original reconstruction underestimates the Si iv emission [the reconstructions were optimized for C iv as described in Allen et al. (2011)]. The region used for calculating the appropriate C iv line scaling is restricted to the interval between line centre and an observed wavelength of 4866 Å. This is to avoid the He ii+O iii absorption and emission complex immediately longward of C iv. The Si iv emission line was shifted by two wavelength bins blueward in the epoch 2 reconstruction relative to the SDSS reconstruction to match the shifted peak of the emission line in the epoch 2 spectrum. Based on the synthetic BAL-based error estimation method described in section 3 of Wildy et al. (2014), the fractional error on the reconstruction for a quasar of this r-band signal-to-noise ratio and redshift is estimated to be 5 per cent of the flux value at each resolution element for all calculations described in this paper. The observed spectra for all four epochs (solid black lines), as well as their final reconstructions (solid red lines), are illustrated in Fig. 1.

Figure 1.

Open in new tabDownload slide

Spectra for epochs 1 to 4 of SDSS J1138+3517, with top panel to bottom panel in order of observation date starting with the earliest. The observed spectrum is in black while the reconstruction of the unabsorbed spectrum is in red. Vertical dotted lines indicate the laboratory wavelengths of the Si iv and C iv emission lines. The vertical dashed line indicates the maximum wavelength used in calculating the C iv emission scaling, emission longward of this point is not accurately reconstructed.

The dramatic variability in both ions’ absorbers across the first three epochs is clear when the spectra are normalized to the reconstructions. The transitions of interest are both doublets, with rest-frame laboratory wavelengths of 1393.76 and 1402.77 Å contributing to Si iv λ1400 and correspondingly 1548.20 and 1550.77 Å for C iv λ1549. Relative to the laboratory-frame rest-wavelength of the red component of each doublet, absorption is significant between 0 and ∼−13 000 km s⁻¹ for Si iv and between 0 and ∼−20 000 km s⁻¹ for C iv, where a negative value indicates outflowing material (blueshifted). As can be seen in Fig. 2, the C iv troughs are deeper than the Si iv troughs at similar velocities. The deepest troughs show the least variability and are located at the lowest (least negative) velocities, as has been reported in previous BAL studies (Lundgren et al. 2007; Capellupo et al. 2011).

Figure 2.

Open in new tabDownload slide

Si iv (upper panel) and C iv (lower panel) absorption regions in red component velocity space. Epoch 1 (SDSS) spectra are in black, epoch 2 (first WHT) spectra are in red, epoch 3 (second WHT) spectra are in green and epoch 4 (third WHT) spectra are in blue. Spectra are normalized to the un-absorbed reconstruction. The Si iv Gaussian components (described in Section 3.3) are numerically labelled in order of increasing outflow velocity.

Due to its relative lack of variation from the previous observation, especially in the Si iv absorption lines where identifiable doublets are present, epoch 4 is left out of the analysis in subsequent parts of Sections 3 and 4. Epoch 4 is instead discussed further in Section 5.3.

3.2 Lower limits for outflow column densities using direct integration

Since the C iv doublet is unresolved, this method is the only means of estimating a lower limit for the column density of this ion. For Si iv, methods involving modelling the components need to be checked for consistency with the lower limit derived here.

3.3 Gaussian components of the Si iv outflow

The Si iv doublet velocity separation is ∼1920 km s⁻¹, allowing some individual components to be identified. Assuming narrow components of quasar absorption lines follow an approximately Gaussian profile, e.g. Hamann et al. (2011), model profiles can be constructed for the doublet lines. Spectral fitting is performed using the specfit package within iraf, requiring three free parameters for each Gaussian component, namely the wavelength at line centre, line full width at half-maximum (FWHM) and line EW. Appropriate restrictions are applied to these parameters, i.e. the EW ratio between the blue and red components of the doublet must be between 2:1 and 1:1, the difference in velocities at line centre must not differ from the expected doublet separation by more than one velocity bin, and the widths must be the same. By using an initial ‘guess’ for these parameters’ values, spectral fitting is performed using chi-square minimization. In total, five Si iv doublets are identified across the three epochs, spanning the range −10 000 ≤ v ≤ 0 km s⁻¹. These are shown in Fig. 3.

Figure 3.

Open in new tabDownload slide

Epoch 1 (top panel), epoch 2 (middle panel) and epoch 3 (lower panel) models of Si iv doublet absorption lines in red component velocity space (related doublet components are shown in the same colour). The black dotted line describes the total model absorption profile, while the red dotted line near the bottom of each panel indicates the difference between the model profile and the observed profile.

Comparing the Gaussian fit to the almost un-blended red component of the lowest velocity doublet in epoch 2 (shown in red in Fig. 3) to the data points at these velocities gives a fractional error in the fit of 0.033. To be conservative a larger error of 5 per cent on each point of all Gaussian profiles is assumed for calculations involving these model profiles. While the EW and occasionally the FWHM change from one epoch to the next, there is no significant change in the line centre velocity and hence no evidence for acceleration of the outflowing material. This lack of acceleration has been noted in previous studies (Hamann et al. 2008; Rodríguez Hidalgo, Hamann & Hall 2011). The best-fitting parameters obtained by specfit for each doublet component are provided in Table 3.

3.4 A pure-partial coverage model

Evidence for pure-partial covering can be obtained by comparing the residual flux velocity dependence of a (almost) saturated absorption line to the profile of 1 − C. Similar profiles would suggest the shape of the absorption profile is determined by the velocity-dependent covering fraction of an opaque absorber rather than intrinsic differences in optical depth. An additional test is to compare the profile of e^−τ to the red component of an absorption doublet. If the profiles do not match, the absorption line profile cannot be due to an absorber which is homogeneous and completely covers the emission region. In order to carry out these tests, a doublet with unblended components is desirable. The closest example available in the data set is absorber 1 in the epoch 2 spectrum. From the values of the EWs of the red and blue components (see Table 3), the blue:red ratio is less than 2:1, indicating that the blue component is probably saturated. After placing the blue component on to the red component velocity grid by interpolation, six velocity bins are found which follow the general shape of the corresponding Gaussian models for both the red and blue components. Five of these points satisfy I_r ≥ I_b ≥ |$I_{\rm r}^{2}$| and are illustrated in Fig. 4.

Figure 4.

Open in new tabDownload slide

The absorber 1 profiles of the red and blue components interpolated on to the red component's velocity grid. The solid red and blue lines represent the Gaussian models, the dotted lines represent the observed fluxes, the green short-dashed line indicates the 1 − C values at these points and the magenta long-dashed line represents the e^−τ values of the red component. Appropriately coloured triangles and squares represent the observed points in velocity space. While the 1 − C profile is a good match for the blue absorber, the e^−τ profile does not match the red profile.

From the profile traced by 1 − C (Fig. 4, dashed green line), it is clear that the saturated blue component could have its profile determined by the uncovered fraction. It is also obvious that the residual intensity of the red component is a poor match to e^−τ. Together these provide strong evidence that partial coverage is a significant effect in this absorber and possibly others in this quasar.

Given that this is the only doublet where meaningful results can be obtained on an individual velocity bin basis, and even in this case it can only be done close to the peak of the absorption components, it is necessary to use the Gaussian models to achieve estimates of column density across a greater proportion of the Si iv absorption for the three epochs. Individual velocity bins in the Gaussian models have no physical realization, so unlike for the previous example it is impossible to construct 1 − C profiles using equation (3). However if for each doublet component model fluxes are averaged over all velocities spanned by the absorption, individual components can be treated as individual velocity bins. This method can subsequently be used to estimate average values of covering fraction and true optical depth for each doublet. Using this true optical depth in equation (1) gives the doublet's column density along the line of sight when N_ion is multiplied by the covering fraction. These values are listed in Table 4. Due to the weakness of absorbers 3 to 5 in epochs 2 and 3 this can only be performed across all absorbers in epoch 1, since at small absorber depths optical depth and covering fraction measurements are unreliable due to their sensitive dependence on doublet component ratio.

3.5 An inhomogeneous absorber model

An alternative to the PPC model is a situation whereby the optical depth varies with spatial location across the line of sight. This was originally developed using a power-law-dependent optical depth model in de Kool, Korista & Arav (2002) and was subsequently investigated by Arav et al. (2005), who found power-law and Gaussian-shaped optical depth dependencies on a one-dimensional spatial location to be inadequate in describing O vi absorbers in Mrk 279. Here, we attempt a similar investigation using a power-law model of the form τ(x, λ) = τ_max(λ)x^a, where τ(x, λ) is the optical depth at position x on the surface of the absorber and at wavelength λ, τ_max(λ) is the maximum optical depth of the absorber and a is the power-law index.

Following Arav et al. (2005), we make two simplifying assumptions to aid calculations without any loss of generality. First, we assume that the surface brightness is uniform across the face of the emission source such that S(x, y, λ) = 1, where x and y are spatial positions in the plane of the line of sight and S is arbitrarily normalized. We also simplify the optical depth spatial dependence by using a one-dimensional model, so that τ(x, y, λ) = τ(x, λ). The x values are fixed to span an interval 0 ≤ x ≤ 1 for simplicity. The observed residual flux at a given velocity bin will then be found using I_res = |$\int _{0}^{1}\,e^{-\tau (x)}\,\mathrm{d}x$|⁠. To assess a broad range of power-law indices, values from 0 to 10 with step size of 0.1 are tested. In order to integrate over the x range, it is effectively divided into 1000 bins, giving dx = 0.001.

Similar to the procedure in Section 3.4, we first examine the plausibility of a power-law model describing the real behaviour of the absorption features by attempting to predict the profile of the blue component of absorber 1 in epoch 2 after interpolation of the blue absorption line profile into the velocity space of the red component. Given a value of a to be tested, the process is carried out as follows: (i) find the value of τ_max which gives the value of the residual red component flux I_r = |$\int _{0}^{1} e^{-\tau _{\rm max}x^{a}}\,\mathrm{d}x$| at each velocity bin; (ii) predict each I_b value using |$I_{\rm b} = \int _{0}^{1}e^{-2\tau _{\rm max}x^{a}}\mathrm{d}x$|⁠, and (iii) compare predicted I_b values to the observed values and calculate the corresponding χ² value over the entire velocity range. This procedure is repeated for all values of a, the value adopted being the one which minimizes the χ² calculation in step (iii). This provides a best value of a = 3.3|$^{+0.4}_{-0.1}$|⁠; the corresponding predicted blue profile is shown in Fig. 5. Finding N_ion requires the knowledge of the average optical depth |$\bar{\tau }$|⁠, where |$\bar{\tau } = \int _{0}^{1}\tau (x)\;\mathrm{d}x$|⁠. As in figs 2 and 3 of Arav et al. (2005), the dependence of the doublet component residual intensities on |$\bar{\tau }$| as well as comparisons with complete homogeneous coverage and PPC models is shown in Fig. 6.

Figure 5.

Open in new tabDownload slide

Inhomogeneous absorber model (red dotted line) of the blue doublet component for a power-law index a = 3.3. The solid blue line is the blue component Gaussian model and the dotted blue line is the observed spectrum in velocity range of the blue component, both of which have been interpolated on to the red component velocity grid. Power-law model error bars are based on the errors of a. Blue triangles and red squares indicate the points on the observed spectrum and inhomogeneous absorber model, respectively.

Figure 6.

Open in new tabDownload slide

Behaviour of doublet components with respect to average optical depth for the power-law model adopted (a = 3.3). The residual intensity of the blue component I_b (blue line) is found by doubling the optical depth of the simulated red component I_r (red line). The |$I_{\rm r}^{2}$| value (magenta line) indicates what the blue component would have been for a homogeneous absorber completely covering the emission source. The green line is (1 − C) for a PPC model given the values of I_r and I_b at a specific |$\bar{\tau }$|⁠. The ratio of |$\frac{\bar{\tau }}{\tau _{\rm {PPC}}}$| indicated by the black line, where τ_PPC is the optical depth of the PPC model predicting the (1 − C) curve, shows the divergence of the power-law model from the PPC model.

From Fig. 5 it is clear that, like the PPC model, an inhomogeneous (power-law) absorber model is a good predictor of the behaviour of the blue component in this particular doublet. As in Section 3.4, this is extended to the Gaussian models with doublet components treated as single points with flux averaged over all velocity bins showing absorption. Interestingly, this leaves the upper bound on a unconstrained for all absorption doublets except absorber 4 at epoch 1 (see Table 5).

Given that the best-fitting power-law index of the Gaussian model for absorber 1 at epoch 2 is not consistent with the measured value of a from the six velocity bins in the observed spectrum, this analysis is considered inconclusive. The constraints on the strongest doublet (absorber 1) suggest that this feature has a comparatively large value of a. At large values of a the dependence of optical depth on position, x, approximates a step-function, and thus cannot easily be distinguished from a PPC model.

4 PHOTOIONIZATION SIMULATIONS

4.1 cloudy setup

Modelling the column densities and ionization fractions of ions in quasar outflows by specifying the shape of the ionizing continuum and various gas (column) densities can provide insight into the physical location of the outflowing plasma and the changing nature of the relationship between the ionizing photon flux and ionization state, especially if ionization changes are the principal driver of absorption line variability. Photoionization models in this investigation are performed using cloudy (Ferland et al. 1998), version c13.02.

A grid of models containing column density data for Si iv and C iv is generated using the input continuum, with grid-points at hydrogen column densities log(N_H/cm⁻²) = 21, 22 and 23, hydrogen number densities log(n_H/cm⁻³) = 5, 7 and 9 and log U values spanning −5.0 to 3.0 in intervals of 0.2 dex. These hydrogen number densities span typical narrow line region to broad line region (BLR) densities, the maximum density being close to the critical density of C iii]. The density is unlikely to be higher than this as this would suggest the absorber originates from within the BLR, conflicting with the fact that individual absorption components are relatively narrow in comparison to the BELs. This is also apparent from the need for the deepest troughs to be absorbing both the BEL and the continuum. The ionization parameter range encompasses all reasonable values given the existence of the absorbers, while the hydrogen column density range was chosen to span a range above the minimum ionic column density for C iv, |$N_{\rm C\,\small {iv}}$|⁠, and below a density at which Thomson scattering becomes significant. This rules out column densities of log(N_H/cm⁻²) ≤ 20 as these generate |$N_{\rm C\,\small {iv}}$| values below the minimum value calculated from direct integration (Fig. 7).

Figure 7.

Open in new tabDownload slide

Predicted C iv column densities as a function of hydrogen number density for given hydrogen column densities. Solid squares indicate the cloudy input hydrogen densities, which are linked by solid lines of constant hydrogen column density. The horizontal dashed line gives the minimum C iv column density predicted by the SDSS observation, which shows the strongest absorption of all epochs. It is clear that log(N_H/cm⁻²) ≤ 20 is ruled out by this limit.

4.2 Application of cloudy models to estimated column densities

The cloudy models can be applied to results from the previous sections in this paper to provide parameters governing the properties of the quasar continuum source which affect the outflowing gas. Initially it is necessary to find the correct normalization of the continuum SED described in Section 4.1, which can be used to find the total ionizing luminosity L_ion and ionizing photon emission rate Q(H) and hence the distance from the continuum source to the inner edge of the outflowing gas (facing the source) and the mass outflow rate. The correct ionizing SED is found by scaling the output to the rest-frame-corrected flux of the observed SDSS spectrum. Output parameters from the simulation include the total ionizing photon flux ϕ(H) which can be used to find the true values of the ionizing luminosity, L_ion, the rate of emission of Lyman continuum photons, Q(H) and the ionizing source distance to the absorbing gas, R.

The calculated value of R will vary depending on the grid parameters (U, N_H and n_H); however since the same SED is used in each case, values of ionizing luminosity and ionizing photon production rate are found which are applied universally; they are L_ion = 4.17 × 10⁴⁶ erg s⁻¹and Q(H) = 9.80 × 10⁵⁶ s⁻¹. The elemental abundances used by cloudy are a solar composition derived from Grevesse & Sauval (1998), Holweger (2001), and Allende Prieto, Lambert & Asplund (2001, 2002).

In order to apply the grid of models calculated in Section 4.1, limits for the column densities of Si iv and C iv must be found. Since there are no resolvable doublets available for C iv, only the lower limits found from direct integration (see Section 3.2) can be used for this ion. For Si iv, the strength of the doublets at epoch 1 means total column densities can be estimated using the PPC model for all five absorbers as listed in Table 4. For epochs 2 and 3, PPC can be used for absorbers 1 and 2 which are the predominant contributors to the total line-of-sight column density at these epochs, plus a contribution from absorbers 3 to 5 using the peak optical depth (POD) method outlined in Section 3.3 modified to take account of partial coverage. Here, we assume that the covering fractions of absorbers 3 to 5 at epochs 2 and 3 are the same as those at epoch 1. The true optical depth can be found once the covering fraction is known. Multiplying the original POD-derived |$N_{\rm Si\,\small {iv}}$| by the covering fraction and the ratio of the integrated real optical depth over the integrated apparent optical depth spanning the velocity range of the red component of the doublet gives an estimate of the true column density. The |$N_{\rm Si\,\small {iv}}$| and their upper and lower limits are shown in Table 6. The total column density lower limits are assumed to be the totals derived from the POD method (listed in Section 3.3). Upper limits on the value of |$N_{\rm Si\,\small {iv}}$| at epoch 2 are comparatively large due to the large uncertainty in the covering fraction of absorber 2; however the best value lies towards the lower end of this interval.

Using the cloudy output at each combination of log N_H and log n_H, plots can be made of predicted ionic column densities across the entire range of ionization parameter U. The span of log U covers the peak of both N_ion for all combinations, resulting in cases where there are two possible regions in which the ionic column density values can be located. Given that only lower limits of the C iv column density can be found, the column densities for this ion are estimated by locating the positions in log U space of allowed |$N_{\rm Si\,\small {iv}}$| and using the equivalent values for C iv. If it is assumed that ionizing continuum changes are the main factor driving variability, it is possible to identify the range of log U which satisfies the ionic column density limits at each epoch. Although there are nine possible combinations of log N_H and log n_H, only two are shown (Figs 8 and 9 ). It is apparent from these figures that the lower region of allowed log U values (where the gradient of log N_ion versus log U is positive) is invalid as the allowed range of |$N_{\rm Si\,\small {iv}}$| for epoch 3 does not permit any |$N_{\rm C\,\small {iv}}$| values above the lower limit. This was found to be true for all nine density/column density combinations, giving one continuous span of allowed log U values in each case.

Figure 8.

Open in new tabDownload slide

Predicted Si iv (solid red lines) and C iv (solid blue lines) column densities as a function of ionization parameter at log(N_H/cm⁻²) = 22 and log(n_H/cm⁻³) = 7. Black, red and green dashed lines with arrows indicate the lower limits of the C iv column density for epochs 1, 2 and 3, respectively. Using the same colour scheme, the solid triangles represent the best values of the Si iv column densities, while the dash–dotted lines show the confidence intervals on these values with the area enclosed being shaded by vertical dotted lines. The thick solid lines at the top show the corresponding log U limits.

Figure 9.

Open in new tabDownload slide

Predicted Si iv (solid red lines) and C iv (solid blue lines) column densities as a function of ionization parameter at log(N_H/cm⁻²) = 23 and log(n_H/cm⁻³) = 9. Black, red and green dashed lines with arrows indicate the lower limits of the C iv column density for epochs 1, 2 and 3, respectively. Using the same colour scheme, the solid triangles represent the best values of the Si iv column densities, while the dash–dotted lines show the confidence intervals on these values with the area enclosed being shaded by vertical dotted lines. The thick solid lines at the top show the corresponding log U limits.

4.3 Outflow properties derived from models

In Table 7 information on mass outflow rate, distance to the inner face of the outflowing gas from the ionizing continuum source and C iv column density is provided at the values of U allowed by the limits on |$N_{\rm Si\,\small {iv}}$| and lower limits on |$N_{\rm C\,\small {iv}}$| at each grid-point. This gives an overview of the allowed range in physical state of the outflow, as well as possible effects a changing ionization parameter could exert. Linear interpolation between the 0.2 dex interval-separated log U values is used to estimate values at intervals of 0.02 dex.

As expected, the range of log U allowed is strongly dependent on N_H rather than n_H, leaving the value of the hydrogen number density uncertain, with the caveat that it probably does not exceed the BLR density. If ionizing continuum flux variations are driving the variability, then given the values of log U and their possible ranges indicated by the third and fourth columns of Table 7 it is possible to estimate how the ionization parameter is changing across the time-scales sampled. Between epochs 1 and 2, log U will be changing by increments of between 0.06 and 0.76 dex with best estimates of ∼0.5 to ∼0.6 dex. For the epoch 2 to epoch 3 interval, log U changes by increments of between 0 and 0.38 dex, with ∼−0.2 dex being the best estimate. This would require a similar change in the ionizing flux reaching the outflowing gas. Changes of similar scale to these best estimates for the change in log U have been observed in the extreme ultraviolet (EUV) of at least one other quasar of similar redshift over a comparable rest-frame time-scale e.g. in Reimers et al. (2005), where the flux at a rest-frame wavelength of 335 Å changed by a factor of ∼3 in a rest-frame time-scale of ∼0.65 yr. Since the EUV is the main contributor to outflow ionization, this indicates our results are plausible.

As indicated in Column 7 of Table 7, the log R values do not change significantly between epochs for given combinations of logN_H and logn_H. The dependence of logR on logN_H is indicated in Fig. 10. It is appropriate to consider how the predicted absorber distances compare to the size of the BLR. Using the method of Kaspi et al. (2005), the BLR radius R_BLR can be estimated using R_BLR = AL^B × 10 lt-days, where L is λL_λ(1450 Å)/10⁴⁴ (erg s⁻¹) and A and B are constants, with a set of values found using the FITEXY method and the mean Balmer line's time lag in Kaspi et al. (2005) to be A = 2.12 and B = 0.496. This set of values is used as the index B is consistent with subsequent studies, e.g. Bentz et al. (2006) which found it to be ∼0.5 This is consistent with the fact that U ∝ Q_H/R², leading to R ∝ L^0.5. The procedure provides a value of R_BLR ∼ 7.6 × 10¹⁷ cm (a few hundred light-days). Comparing this with the values found for R in Table 7 suggests that the outflowing gas is outside the BLR in eight out of nine cases, the exception being the grid-point log(N_H/cm⁻²) = 23, log(n_H/cm⁻³) = 9.

Figure 10.

Open in new tabDownload slide

Dependence of log R on log N_H for constant mass outflow rates. Values are from Table 7.

The actual accretion rate can be calculated from the bolometric luminosity, which can be found using the luminosity at 1450 Å by applying L_bol ∼ 4.36λL_λ (at 1450 Å) (Warner, Hamann & Dietrich 2004). For a bolometric luminosity of L_bol ∼ 9.5 × 10⁴⁶ erg s⁻¹, the resultant value of the accretion rate in terms of that at Eddington is ∼0.15|$\dot{M}_{\rm Edd}$| (∼17.0 M_⊙ yr⁻¹), meaning that this object's accretion rate is a substantial portion of the value at the Eddington limit. This allows the ruling out of two combinations of hydrogen column density and hydrogen number density [log(N_H/cm⁻²) = 22, log(n_H/cm⁻³) = 5 and log(N_H/cm⁻²) = 23, log(n_H/cm⁻³) = 5] due to their excessively high outflow rates compared to the accretion rate, as indicated in Table 7. Comparing the remaining allowed values suggests the mass outflow rate is approximately between 1 and 200 per cent of the mass accretion rate, corresponding to a kinetic luminosity range of ∼3.3 × 10⁴² to 6.6 × 10⁴⁴ erg s⁻¹. The upper limit of this range gives a kinetic energy to luminosity ratio of ∼0.7 per cent, high enough to have a significant impact on the evolution of the surrounding galaxy (Hopkins & Elvis 2010); however definite conclusions regarding the outflow's importance in this regard are impossible given the weakness of the constraint.

5 DISCUSSION

5.1 Are changes in the outflow geometry responsible for variability?

Changes in covering fraction within the context of the PPC model are one possibility for the dramatic variability seen in this quasar's absorption lines. From Section 3.4 (Fig. 4) it is clear that the deepest parts of the blue absorber 1 trough in epoch 2 are accurately described by a saturated absorber whose depth is determined by the fraction covering the emission source. However, Section 3.5 (Fig. 5) indicates that an inhomogeneous absorber modelled over the same velocity range is also an excellent predictor of the same trough. Further insight can be gained from the Gaussian model profiles, specifically absorbers 1 and 2 which are well-defined and hence can be studied using both PPC and inhomogeneous models. From Table 5 it appears that these features are not well described by an inhomogeneous absorber, since, given the maximum power-law index of 10, many of the calculated values are unconstrained. Large power indices produce line-of-sight geometries which resemble a step-function, which is equivalent to a PPC absorber. The difference between the power-law index which matches the deepest part of Absorber 1 in epoch 2 (a = 3.3|$^{+0.4}_{-0.1}$|⁠) and the value in Table 5 is also much greater than the uncertainty. This leads to the conclusion that an inhomogeneous absorber is probably not the predominant influence on absorption line profile shape.

If it is assumed that absorption takes place within the context of PPC, it is possible to determine how the covering fraction has changed between epochs using the Gaussian profile models. From Table 4 it is clear that there is no significant change in C for absorber 1 between epochs (changes being within 1σ for σ ≤ 0.042). There is also no discernible change in the covering fraction of absorber 2. Conclusions regarding covering fraction variations in absorber 2 are weak for two reasons. First, the epoch 2 doublet is entirely unsaturated and therefore the covering fraction estimation is completely unconstrained. Secondly, when comparing the epoch 1 value (C = 0.187 ± 0.017) with the epoch 3 value (C = 0.275 ± 0.165), the large uncertainty in the latter allows a significant range which prevents any strong conclusion regarding the presence of covering fraction variations. Given the lack of significant variability in either the strongly constrained absorber 1 or the weakly constrained absorber 2, we conclude that there is no evidence for covering fraction being the dominant mechanism behind the variability.

5.2 Evidence for ionization changes driving variability

Ionization fraction changes within the outflowing gas can be generated either by intrinsic changes in the ionizing radiation emitted from the central ionizing continuum source or by shielding of this radiation by gas in the region between the accretion disc and the BAL gas. There are several indicators that suggest ionization fraction changes are implicated in the absorption line variability in this object, the strongest being coordinated changes seen across large velocity separations. A similar result was noted in Trevese et al. (2013) regarding the quasar APM 08279+5255, where they observed coordinated variability in two Gaussian-modelled components of a C iv BAL separated by ∼5000 km s⁻¹. Large velocity separations between absorption troughs imply that the absorption arises in physically distinct regions, hence making covering fraction changes an unlikely explanation due to the level of coordination required across such large separations (e.g. Filiz Ak et al. 2013).

This phenomenon is especially evident in the C iv absorption troughs (see Fig. 2). The visibly saturated absorption trough centred at ∼−4000 km s⁻¹ undergoes weakening between epochs 1 and 2 before strengthening again in epoch 3, which is a pattern also seen in the features at ∼−9500 km s⁻¹ and ∼−18 500 km s⁻¹. Since there is a lot of blending of the Si iv components, the same effect is not as obvious in this ion. However, it is instructive to look at the behaviour of the modelled Gaussian Si iv doublets in Fig. 3, where the same effect is apparent.

No part of the power-law continuum contributing to this spectrum is observable shortward of Ly-α; however it is reasonable to expect that it extends into the EUV portion of the spectrum (Kriss et al. 1999), which is largely responsible for the flux which controls the ionization fractions of the plasma in the accretion disc's environment. Therefore, variability in the observed continuum could indicate variability in the ionizing radiation from the accretion disc. There is evidence for an excess in the blue part of the spectrum during epoch 2 (first WHT observation) when compared to the prior and subsequent observations, as indicated in Figs 11 and 12, which could indicate an increase in the ionizing continuum flux at this epoch. At the third epoch, the blue excess appears to disappear, with the flux in this range appearing to follow the continuum from epoch 1 once again.

Figure 11.

Open in new tabDownload slide

Upper panel: epoch 1 observation (black line) and epoch 2 observation (red line) with power-law continuum (dashed black line) fitted to the epoch 1 spectrum. Vertical dotted lines indicate the boundaries of where the blue excess is apparent. The thick black horizontal line indicates the dichroic overlap region in the WHT observation. Lower panel: same as upper panel but with red line indicating the epoch 3 observation.

Figure 12.

Open in new tabDownload slide

Epoch 1 observation (black), epoch 2 observation (red) and epoch 3 observation (green), normalized to the epoch 1 continuum shown in Fig. 11 over the blue excess range, with dotted lines representing the errors. The epoch 2 spectrum appears to lie almost everywhere above the epoch 1 continuum.

Evidence from the cloudy simulations also supports an ionizing continuum change driving the variability. Given the decrease in the Si iv column density between epochs 1 and 2, an increase in ionization parameter must occur between these observations (see Figs 8 and 9), with an estimated log U change of ∼0.5 to ∼0.6 dex. An increase in log U agrees with the observed steepening of the continuum at the blue end of the spectrum between epochs 1 and 2, assuming the continuum extends to ionizing energies. There is a large overlap in predicted Si iv column densities between epochs 2 and 3; however the epoch 3 range lies at the higher end of the epoch 2 range, with the best epoch 3 value being 0.3 dex above the best epoch 2 value. A log U change of −0.2 dex is needed to drive this increase in Si iv column density between epochs 2 and 3. A decrease in ionization parameter between these two epochs is supported by the apparent disappearance of the aforementioned blue excess in the continuum.

5.3 Fourth observation

The quasar SDSS J1138+3517 was observed again by the WHT in 2014 May and is included in Fig. 2. In contrast to the earlier WHT observations, the Si iv components are largely unchanged from the previous epoch's spectrum. The deep C iv features centred at ∼−4500 and ∼−9500 km s⁻¹ do not show changes comparable to previous variability when going from epoch 3 to epoch 4. However, the weaker feature at ∼−18 500 km s⁻¹ shows a significant strengthening, with its profile becoming comparable to that of the epoch 1 observation. The absorption line centred at ∼−13 000 km s⁻¹, which is very conspicuous at epoch 1 but almost unidentifiable at epochs 2 and 3, shows a slight increase in strength; however the significance of this change is only 1.08σ.

Although these changes from epoch 3 to epoch 4 are in the weakest C iv components and therefore do not significantly change the minimum column density estimate made in epoch 3, the localization in velocity space of this strengthening to a range of ∼6000 km s⁻¹ is interesting. If the increase in absorption is only considered to be significant in the stronger of the two features, then the region over which the change occurs is reduced to ∼4600 km s⁻¹. A single varying trough allows the possibility of a change in the covering fraction of the absorber to explain the phenomenon. This trough is also interesting in that it does not show as great a change between epoch 2 and 3 as the feature centred at −9500 km s⁻¹. This could be due to either (a) the line being the result of a moderately saturated outflow whose peak trough depth describes the peak covering fraction, or (b) the outflow is located in a region where it is relatively shielded from the ionizing flux variations (the large velocity separation from the other troughs suggests it is in a physically distinct location.).

6 SUMMARY

This following provides a summary of the conclusions that can be drawn from this investigation.

• When comparing to large-sample statistical studies of BALQSOs (Filiz Ak et al. 2013, 2014; Wildy et al. 2014), BALs in SDSS J1138+3517 are exceptionally variable, especially in the 360 d rest-frame time interval between epoch 1 and epoch 2. When measuring BAL variability using the method of Filiz Ak et al. (2013, 2014), the most variable C iv BAL in this object undergoes a change of ΔEW = −13.8 Å. When compared to fig. 29 of Filiz Ak et al. (2013), this change in ΔEW is far outside the range predicted by their random-walk model at an equivalent time interval. This could indicate that SDSS J1138+3517 belongs to a rare class of highly variable BALQSOs whose variability mechanism is not the same as that driving the random-walk behaviour observed in Filiz Ak et al. (2013).

• When applying Gaussian fitting to model the Si iv outflows in Section 3.3, it is found that the components of individual doublets rarely exhibit the expected 2:1 blue to red optical depth ratio and also do not reach zero intensity at line centre. This suggests saturation of the blue component combined with spatially dependent line-of-sight optical depth (either inhomogeneous coverage or PPC).

• Although both inhomogeneous power-law and PPC models successfully predict the central velocity range of the deepest saturated blue trough in epoch 2, modelling the absorption as Gaussians over the entirety of the velocity range suggests a PPC model is more successful. This is due to the fact that the best values of the power-law index suggest a very steep slope, therefore approximating a step-function which is the mathematical equivalent of PPC. The fact that both models can be successfully applied over a limited range in absorption line velocity is something that should be noted in future investigations which attempt to distinguish these models.

• There is no evidence for a change of covering fraction across the epochs when calculating C from the (Gaussian modelled) doublets representing the two deepest absorbers, suggesting line-of-sight covering changes are not the dominant mechanism driving the absorption line variability. However, the weak C iv feature centred at ∼−18 500 km s⁻¹ may have undergone a covering fraction change between epochs 3 and 4.

• There is evidence that ionization changes play a significant role in the variability, most notably from the coordinated changes in both Si iv and C iv absorption troughs separated widely in velocity. There is also a hint that the ionizing flux may have increased in epoch 2, given the apparent steepening of the continuum in the blue part of the spectrum. This conclusion would be supported if AGN have steeper α_ox at higher luminosities as indicated in Vignali, Brandt & Schneider (2003). The photoionization simulations are consistent with this scenario, since the column densities of both ions decrease with increasing ionization parameter in the log U range spanned by the outflow. The best values for the changes in ionization parameter are within the scope of quasar EUV variability.

• The location of the inner face of the outflowing gas is not well constrained, being located between approximately 0.1 and 200 parsec from the continuum source. This range places it somewhere in a region spanning the BLR radius to narrow line region distances. However, the lowest velocity (and deepest) C iv trough must be absorbing both BLR and continuum emission, meaning that this outflow is located outside the BLR. Given the width of the absorption components, which are substantially less than BLR widths, and the value of R_BLR, it seems likely that the inner face is at least ∼1 parsec from the continuum source.

• The magnitude of the mass outflow rate could be a substantial fraction of or even exceed the mass accretion rate. However due to the uncertainties on R, strong constraints on the true value of the mass outflow rate are not achieved. The derived kinetic luminosity range of ∼3.3 × 10⁴² to 6.6 × 10⁴⁴ erg s⁻¹ indicates that feedback effects such as quenching of star formation may be possible if the true value is near the upper limit of this range.

This work is supported at the University of Leicester by the Science and Technology Facilities Council (STFC) and is based on observations made with the WHT/ISIS operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias and spectroscopic observations from Data release 6 of the SDSS. We wish to thank our anonymous referee for their useful comments and thorough review of the draft.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

REFERENCES