Investigating the Origin of the Absorption-line Variability in the Narrow-line Seyfert 1 Galaxy WPVS 007

WPVS 007 has been observed with UV COS spectroscopy with HST five times. The details of the observations including instrument, grating, date, wavelength range, spectral resolution, exposure time, and signal-to-noise ratio (S/N) at 1450 Å are provided in Table 1. UV spectral observations with HST FOS in 1996 and FUSE in 2003 exist in the archive but were not included because they did not have comparable S/N (in the case of HST FOS) or wavelength range (in the case of FUSE) for a comparative analysis. Interestingly, there is no evidence for the UV BALs in 1996, though the narrow, low-velocity mini-BAL is still evident in the spectrum (Leighly et al. 2009).

The final data products for all epochs were downloaded from the archive in 2020 January for analysis. No recalibration was necessary. Spectra were adjusted to rest-frame wavelengths using a redshift of 0.02861 adopting the NASA/IPAC Extragalactic Database redshift value.

The objective was to accurately model the shape of the available UV absorption lines in order to determine the physical properties of the outflow that cause them. Modeling required that features have sufficient S/N. We used the spectra in the 1065–1600 Å range to include the P V λ λ1118, 1128 and the C iii*λ1175 multiplet BALs at the blue end of the wavelength range and the C IV λ λ1548, 1551 line at the red end of the wavelength range. The Ciii] complex at 1900 Å has low S/N, and so we did not include this in our analysis.

Strong geocoronal lines are present from 1173–1189 Å and from 1257–1273 Å. We masked these sections from the analysis. We identified possible Galactic absorption lines from C II and Si II. Candidate Galactic lines were identified as narrow absorption lines that showed no change in optical depth between epochs. The region around the C II λ1334 line was masked in all epochs. We identified Galactic lines of Si II overlapping the C iii* trough, most obvious in the observation from 2015 when the C iii* BAL was weak. These lines are likely Si II λ1190 and λ1193 that have been shifted to 1157 Å and 1161 Å, respectively. Because these latter features overlapped with BALs, we included them in the modeling rather than masking them out. There is a detector gap in the 2010 epoch from 1139 Å–1230 Å in the rest frame.

Swift has observed WPVS 007 since 2005 October typically with a cadence of once per week (Grupe et al. 2013). Observations prior to 2013 April are listed in Grupe et al. (2008, 2013). The Swift X-ray telescope (XRT; Burrows et al. 2005) was operating in photon-counting mode (Hill et al. 2004), and the data were reduced by the task xrtpipeline version 0.12.6., which is included in the HEASOFT package. Source counts were selected in a circle with a radius of 24.8'' and background counts in a nearby circular region with a radius of 247.5''. Due to the extremely low X-ray count rate of WPVS 007 of about 1 × 10⁻⁴, it cannot be detected in a single observation. An X-ray detection of WPVS 007 with Swift requires stacking the data of several years together. The UV-optical telescope (UVOT; Roming et al. 2005) data of each segment were coadded in each filter with the UVOT task uvotimsum. Source counts in all six UVOT filters were selected in a circle with a radius of 5'' and background counts in a nearby source-free region with a radius of 20''. UVOT magnitudes and fluxes were measured with the task uvotsource based on the most recent UVOT calibration as described in Poole et al. (2008) and Breeveld et al. (2010). The UVOT data were corrected for Galactic reddening (E_B−V = 0.012; Schlegel et al. 1998). The correction factor in each filter was calculated with Equation (2) in Roming et al. (2009), who used the standard reddening correction curves by Cardelli et al. (1989). Due to new UVOT calibrations files, it was necessary to reanalyze all previously published data.

We are providing updated photometry from continued Swift monitoring since 2013 April (Grupe et al. 2013). This observing campaign is the only regular monitoring of such an unusual low-luminosity source with a fully developed BAL. In Figure 1, we show the results of the Swift monitoring campaign carried out over a period of 16 yr. The Swift UVOT fluxes were corrected for Galactic extinction. As shown in Leighly et al. (2015), the 2015 occultation event resulting in the near disappearance of the BAL shows a corresponding minimum in UV flux. The 2010 BAL was found when WPVS 007 was in a high state at the brightest magnitude.

Zoom In Zoom Out Reset image size

Figure 2 shows the SED using observations taken near the time of HST COS spectral observations, and is an updated version of the left panel of Figure 2 from Leighly et al. (2015). Although the SED for 2010 (where the BAL was at maximum absorption) shows higher flux and 2015 (during the occultation event) shows reduced flux, the SEDs for 2013 and 2017 December with variable BALs are nearly identical. These observed changes are in contrast to the 2013 June and December epochs that show variation in the SED with almost no change in the BALs.

Zoom In Zoom Out Reset image size

Figure 3 shows the evolution of WPVS 007 in UV color–magnitude space, an updated plot from Figure 1 in Leighly et al. (2015). There is a clear correlation between color and magnitude. The 2010 and 2015 observations were taken when WPVS 007 occupied two extremes in the color–magnitude diagram with the NLS1 appearing bluest in 2010 and reddest in 2015. The 2013 and 2017 colors are intermediate. In Figure 1 of Leighly et al. (2015), the authors showed that variable reddening is a better fit to the observed changes in color and brightness observed from the Swift photometry than intrinsic optical/UV spectral variability.

Zoom In Zoom Out Reset image size

SimBAL (Leighly et al. 2018) is a novel forward-modeling and spectral-synthesis software package that has been used to model quasars with broad overlapping absorption troughs, leading, for example, to the discovery of the most powerful quasar outflow to date (Choi et al. 2020). SimBAL uses grids of ion column density created using Cloudy (Ferland et al. 2017) to produce a synthetic spectrum from available atomic data and Markov Chain Monte Carlo (MCMC) python package emcee (Foreman-Mackey et al. 2013) to step through the parameter grids.

To generate the synthetic spectrum at each step, SimBAL uses the ionization parameter, density, covering factor, and hydrogen column density. The column density is parameterized as $\mathrm{log}({N}_{{\rm{H}}})-\mathrm{log}(U);$ here U is the unitless ionization parameter defined as $U=\tfrac{Q}{4\pi {R}^{2}{n}_{{\rm{H}}}c}$, where Q is the number of ionizing photons per second, R is the distance of the BAL from the SMBH, n_H is the hydrogen density, and c is the speed of light. Partial covering is now widely accepted as endemic for BAL absorption lines (e.g., Barlow & Sargent 1997; Arav et al.1999), and can manifest as non-black saturation of BALs. Often this absorption is modeled as a step-function in column density, where the absorber covers part of the continuum uniformly, but light from the remaining uncovered region comes through unobstructed. SimBAL uses a power-law partial-covering model (see de Kool et al. 2002; Arav et al. 2005; Sabra & Hamann 2005; Leighly et al. 2019) where $\tau ={\tau }_{\max }{x}^{a}$ with x as the projection of the 2D region onto one dimension with a value between 0 and 1 and $\mathrm{log}(a)$ is the parameter fit by SimBAL. Power-law partial covering assumes a power-law distribution of optical depth, where the power-law index is a fit parameter for SimBAL. When the index value is 0, this is equivalent to the homogeneous partial-covering case, e.g., a sharp-edged, solid absorber. SimBAL uses power-law partial covering because it is computationally tractable and explains the observed difference in covering fraction between high- and low-ionization lines (e.g., Hamann et al. 2001). See Leighly et al. (2018) for further details about SimBAL and the spectral-synthesis method of modeling BAL quasar absorption and Leighly et al. (2019) for further discussion of the covering fraction.

An advantage of SimBAL over traditional analysis methods is that the solution is self-consistently constrained by all of the lines simultaneously, as well as from the upper limits of the absorption lines that are not detected. Therefore, in a velocity-resolved analysis (such as the analysis presented here), the physical parameters of the gas are constrained by multiple lines. This allows for excellent constraints on physical conditions of the gas due to the diagnostic power of lines with different physical properties such as oscillator strengths and abundances. The first velocity-resolved analysis using SimBAL was the of low-ionization BAL source SDSS J085053.12 + 445122.5 (Leighly et al. 2018, 2019). Applied to multiepoch data, this methodology offers the capacity to study the change in the physical outflow properties over time as a function of velocity.

In Table 2, we list key lines modeled in absorption, including their ionization potential and primary diagnostic power for determining physical conditions of the gas. Figure 4 shows the differences in absorption between epochs after continuum normalization. To find a generic absorption model that would work for all epochs, we first look at the overall absorption structure in all epochs and at which BALs and BAL subfeatures change over time. Subsequently, a starting model for absorption lines was chosen taking into consideration the morphology of the lines, total width of the absorption lines in all epochs, and variation within the structure of individual BAL lines between epochs.

Zoom In Zoom Out Reset image size

The ionizing continuum used in the Cloudy grids is a relatively soft quasar SED (Hamann et al. 2011), which we use because WPVS 007 has been consistently found to be X-ray weak, with a soft X-ray spectrum (Grupe et al. 1995, 2007, 2008, 2013). We also assumed an enhanced metallicity (Z = 3 Z_⊙) consistent with previous findings that BALQs have metallicities Z ≳ Z_⊙ as measured using emission-line ratios (Hamann et al. 2002, but see also Temple et al. 2021).

Of the physical outflow parameters listed in Table 2, density is the hardest to constrain with only the C iii* line. In order to constrain density, we need absorption from the collisionally excited state of an ion (e.g., Gabel et al. 2005; Arav et al. 2013). In this case, C iii* is the only ion within this wavelength range producing absorption from an excited state. Lines with high ionization potential, such as N V, are sensitive to changes in ionization parameter. As was discussed in Leighly et al. (2018), the P V line is very useful as a constraint on column density because it is rarely optically thick due to the low abundance of P V (Hamann 1998) and also constrains the ionization parameter as a high ionization is required to produce P V (see also Leighly et al. 2009). We observe the P V line in all COS spectra, which is an indication that the much more abundant C IV will likely be saturated (Borguet et al. 2012).

The kinematic differences and opacity differences (in the case of C iii*) required us to look for ways to tie or fix parameters in the model. As a first step, we look for ways to group structures within the BAL. Based on relative velocity widths of various lines, we divided the velocity profile into four main regions, a high-velocity gas component with all absorption between ∼−12,000 km s⁻¹ and ∼–8000 km s⁻¹, a medium-velocity component consisting of absorption between −8000 and −5000 km s⁻¹, and two low-velocity components with absorption between ∼ −5000 km s⁻¹ and ∼ −3000 km s⁻¹ and ∼ −3000 km s⁻¹ and ∼ −600 km s⁻¹, which we call low (1) and low (2), respectively. The C iii* line (that we rely on to constrain the density; see Table 2) absorbs only at low velocities (∼ −5000 km s⁻¹ to ∼ −600 km s⁻¹) in all epochs, but due to differences in absorption between low (1) and low (2) in the 2015 epoch compared to other epochs, we subdivided this velocity group further. The high-velocity gas is only present for some high-ionization lines such as C IV and N V in the 2010 and 2017 epochs, and not low-ionization lines such as C iii* or C II, and is not present in P V. The mini-BAL feature at velocities of ∼ −600 km s⁻¹ overlaps the C IV emission line (with strong absorption for other high-ionization lines of Si IV and N V). The mini-BAL absorption line is deeper than the BAL and does not show variation between epochs. The observed lack of variation implies a larger covering fraction for the mini-BAL than the BAL. For the remainder of the paper, we discuss absorption variability as a function of the velocity groups defined.

We used the tophat setting in SimBAL (see Leighly et al. 2018, 2019; Choi et al. 2020) to model the spectra. The tophat accordion model uses rectangular bins of equal velocity width to span the BAL. Both the maximum velocity and width per bin are fit as parameters in SimBAL, but the number of bins remains fixed. The starting bin size was chosen to be 950 km s⁻¹, roughly half the separation width of the Si IV doublet, or the velocity separation between the N V lines. We used C IV to select the number of bins necessary to model absorption across all epochs. The number of bins for each epoch was chosen so that variability in each subfeature of key BALs could be studied between epochs. Specifically, we wanted to ensure that the number of bins in each epoch resulted in the final velocities of each bin roughly aligned between epochs. We chose the smallest number of bins to achieve this objective for each epoch. Using the velocity groups defined above, we use 12 bins to span absorption from velocities of −13,000 km s⁻¹ to −1200 km s⁻¹ in 2010. For the 2013 epochs, we used seven bins to span velocities from −8000 km s⁻¹ to −1200 km s⁻¹, the 2015 epoch was fit with eight bins spanning from −9000 km s⁻¹ to −1200 km s⁻¹, and the 2017 epoch was fit with 11 bins spanning velocities from −12,000 km s⁻¹ to −1000 km s⁻¹. We use two tophat bins to fit the mini-BAL absorption. These tophats had widths of 300 km s⁻¹ spanning velocities from −600 km s⁻¹ to 0 km s⁻¹. Leighly et al. (2018) looked at the difference in final fit parameters when varying the total number of bins that span the BAL and found that the total number of bins used to model the absorption did not significantly affect the values of the final fit parameters.

To reduce the number of free parameters, we tied some absorption parameters between SimBAL bins within the same velocity group. Following the analysis of Leighly et al. (2018), each bin has a unique value of $\mathrm{log}({N}_{{\rm{H}}})-\mathrm{log}(U)$ and partial-covering parameter $\mathrm{log}(a)$. Ionization parameter is tied between bins in a given velocity group (with unique ionization parameters assigned to low (1) and low (2)) to reduce the number of fit parameters (and therefore computation time). Ionization parameter is more difficult to constrain with the lines available in the spectrum, and we further found that freeing ionization parameter between bins did not improve the fit. With C iii* as the only density-sensitive line present in all spectra, we fit only one density value for each epoch. The low-velocity bins overlap with the C iii* line and are constraining the density value for this model. We tested whether freeing the density between velocity groups improves the fit, but found that this failed an F-test for significance.

As a first step, we fit each epoch independently, letting all outflow, continuum, and emission-line parameters vary between epochs. In Figures 5 and 6, we present the fits and best-fit parameters, respectively, for the BAL and mini-BAL for all epochs. We note that the continuum in the 2013 epochs shows an uneven blue wing for the C IV line. We tested the impact of this shape on our final results by fitting a separate model for the 2013 December epoch with additional constraints to force a smooth shape for the final C IV emission line. The resulting final parameters were consistent with the original fit as expected because the C IV line is saturated, and therefore is not strongly constraining on the best-fitting model parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The physical parameters as a function of velocity are shown for each HST COS epoch fit in the right panel of Figure 6. Density is not shown for the mini-BAL as it could not be constrained. The left panel of Figure 6 shows synthetic ion models produced by SimBAL using the final best-fit models for each epoch as a function of velocity. Dashed lines are used to show the corresponding BAL features for each velocity group. Velocity bins between main velocity groups were given unique ionization parameter values and not tied. Although the density was fit for the 2010 epoch (and shown in Figure 6), we do not present any parameters derived from the density value for 2010 as no density-sensitive lines were present in the bandpass. We note that the high-velocity component in 2017 has a lower ionization parameter than the low-velocity gas and yet shows no evidence of absorption from low-ionization line C iii*. The C iii* line is sensitive to more than ionization state alone. SimBAL is fitting the column density, ionization parameter, and covering fraction simultaneously for this component, all of which contribute to the opacity produced (or not produced in the case of the high-velocity gas component) in the synthetic spectrum.

We simultaneously fit two epochs at a time, tying all absorption parameters between epochs except one to see if a single parameter can explain the variability between epochs. These simultaneous fits of WPVS 007 are the first time that more than one epoch of observation was fit simultaneously using SimBAL. This analysis was done with the 2013, 2015, and 2017 December HST COS observations. The 2010 observation was excluded due to missing observations of C iii* and N V lines, and we chose the 2013 December observation over the 2013 June observation as it is chronologically closer to the 2015 and 2017 observations. The 2013 and 2017 epochs were chosen first due to the comparable signal-to-noise and little change in continuum or line parameters between the two epochs. In the end, simultaneous fits of 2013–2017 (not including 2015), 2013–2015, and 2015–2017 epochs were carried out.

The scenarios tested are as follows: (1) there is a change in ionization state of the gas, which will be reflected as a change in ionization parameter (U); (2) there is a change in column density; or (3) there is a change in partial covering ($\mathrm{log}(a)$). As SimBAL fits the column density with respect to the hydrogen ionization front ($\mathrm{log}({N}_{{\rm{H}}})-\mathrm{log}(U)$) instead of a total column density, we tied the column density for scenario (1) listed above.

In Figure 7, we show the results of the 2013 and 2017 simultaneous fitting for the three test cases outlined above. The models fail to fit the full depth of the P V and Si IV BALs in both of the following scenarios: (1) where we assume that the ionization parameter is causing the variability, and (2) where we assume that the column density alone is controlling the variability. In contrast, all absorption features of both epochs can be fit by varying the partial-covering parameter, $\mathrm{log}(a)$ only. The variation in absorption lines can be therefore explained as largely driven by a change in the amount of partial covering.

Zoom In Zoom Out Reset image size

Given that the $\mathrm{log}(a)$ parameter is the main driver for absorption-line variability between epochs, we further investigate if freeing either ionization parameter or column density in addition to partial covering improves the observed fit from the $\mathrm{log}(a)$ fits alone; these models fail the F-test for improvement to the model given the additional parameters.

The changing covering fraction best explains the changes between 2013 and 2017 and is less satisfactory for the ultra low state in 2015. We test whether these results are consistent with the 2015 data by fitting the 2015 spectrum with the best-fit values of density, ionization parameter, and column density from the combined 2013–2017 fit and allowing only the continuum parameters, velocity, and $\mathrm{log}(a)$ parameters to vary. Using an F-test, we find that the fit is improved statistically when all absorption parameters (not just $\mathrm{log}(a))$ are allowed to vary. The residuals show that the difference between the models (where only $\mathrm{log}(a)$ is free versus where all absorption parameters are free) is mostly overlapping the emission lines and mini-BAL area, indicating that there may be either some small change in the mini-BAL in 2015 or that there is some degeneracy between the continuum parameters and the mini-BAL absorption parameters. As far as differences between the BAL models, the difference comes in a slightly poorer fit for C iii* in the two lowest velocity bins where we see the BAL disappear in 2015. The change in covering fraction is not able to explain the low opacity in 2015 for these lowest velocity bins, which may indicate that there are slight variations in other parameters at low velocities between these epochs that are not being taken into account when only $\mathrm{log}(a)$ is allowed to vary.

We wanted to see if we can determine the cause of variability using the individual fits of each epoch to confirm the result we obtained with the simultaneous fits of the 2013 and 2017 epochs. To do this, we look for correlations between the change in depth of the BAL and the change in physical parameters from the final SimBAL fits presented in Figure 6. Throughout this analysis, we used a change in bin depth between epochs to measure the absorption strength of the BAL at a certain velocity. As our velocity bins are ∼ 1000 km s⁻¹, this change in bin depth is comparable to the statistic ΔA used by Capellupo et al. (2011, 2012, 2013) to quantify changes in BALs over time. Their ΔA parameter is the change in the fraction of the continuum-normalized flux removed over 1000–2000 km s⁻¹ velocity bins (Capellupo et al. 2011).

Using the final synthetic spectra from each single-epoch fit, we produce bin-by-bin synthetic spectra for the C IV, C iii*, N V, Si IV, and P V lines and measure the change in bin depth per bin, for each ion between the same bins in different epochs. We examined the change in bin depth between each pair of epochs always with the form (later epoch − earlier epoch) for each bin. Combinations of epochs are 2010–2013, 2010–2015, 2010–2017, 2013–2015, and 2013–2017. We also calculate the change in each of the physical parameters (ionization parameter, covering parameter, $\mathrm{log}({N}_{{\rm{H}}})-\mathrm{log}(U)$, and column density) per bin between pairs of epochs. There are fewer data points for ionization parameter ($\mathrm{log}(U)$) as all bins of a particular velocity group have the same $\mathrm{log}(U)$ value, and so we only include one data point per unique value of ionization parameter. The bin depth for those data points is the mean of the change in bin depth for all bins that share the same ionization parameter. We exclude one outlier data point with a change in $\mathrm{log}(U)$ of less than −2, as this data point came from the difference in ionization parameter from the high-velocity gas in 2010 and 2017, which may be poorly constrained due to a lack of absorption from C iii* and P V in both epochs. We additionally exclude three data points with a $\mathrm{log}(U)$ change greater than 2.1, as all three data points corresponded to the changes in the lowest-velocity component from the 2010 epoch where the 2010 data has almost no opacity. We use the scipy.stats package to calculate the Spearman-rank correlation coefficient and corresponding p-value in each case. The results are shown in Figure 8.

Zoom In Zoom Out Reset image size

The high-ionization lines (C IV, Si IV, N V, and P V) all show a negative correlation between the $\mathrm{log}(a)$ parameter and the change in bin depth (−0.88 for C IV, −0.65 for Si IV, −0.88 for N V, and −0.34 for P V) with very small p-values (1.8 × 10⁻¹⁵, 1.3 × 10⁻⁶, 2.9 × 10⁻¹⁰, and 0.02 for C IV, Si IV, N V, and P V, respectively) indicating that this correlation is strong, and a decrease in covering fraction manifests as a decrease in line depth across all of the high-ionization lines. The low-ionization line (C iii*), which is unsaturated, did not show the correlation with $\mathrm{log}(a)$, but showed a positive correlation with column density (0.67 with a p-value of 7.1 × 10⁻⁵) and $\mathrm{log}({{\rm{N}}}_{H})-\mathrm{log}({\rm{U}})$ (0.65 with a p-value of 1.6 × 10⁻⁴). These correlations may indicate that variability in different lines may be controlled by different parameters. This analysis supports the result obtained from Section 4.1 that covering fraction is the primary driver of the variability observed.

4.2.1. Regression

The physical wind parameters that are directly fit by SimBAL can be used to calculate derived physical properties of the outflow including total hydrogen column density, mass-loss rate, radius, and kinetic luminosity, given in Table 3. We discuss in further detail below how the properties were calculated and their implications. For column density and radius, these values are given both by velocity group and in total.

Table 3. Derived Parameters from the Best-fit SimBAL Model

Parameter	Velocity Group	2010	2013 June	2013 Dec.	2015	2017
$\mathrm{log}({{\rm{N}}}_{{\rm{H}}})({\mathrm{cm}}^{-2})$ a	High	${22.31}_{-0.29}^{+0.21}$	⋯	⋯	⋯	${19.35}_{-0.18}^{+0.32}$
	Medium	${20.23}_{-0.30}^{+0.26}$	${21.33}_{-0.19}^{+0.19}$	${21.11}_{-0.24}^{+0.23}$	${21.41}_{-0.42}^{+0.31}$	${21.96}_{-0.32}^{+0.18}$
	Low (1)	${21.34}_{-0.48}^{+0.22}$	${21.66}_{-0.08}^{+0.09}$	${21.47}_{-0.11}^{+0.11}$	${22.17}_{-0.22}^{+0.13}$	${22.39}_{-0.14}^{+0.14}$
	Low (2)	${18.61}_{-0.06}^{+0.06}$	${21.36}_{-0.07}^{+0.07}$	${21.16}_{-0.09}^{+0.09}$	${21.73}_{-0.49}^{+0.60}$	${22.47}_{-0.11}^{+0.12}$
	Total BAL	${22.76}_{-0.20}^{+0.17}$	${22.28}_{-0.06}^{+0.06}$	${22.11}_{-0.08}^{+0.09}$	${22.59}_{-0.32}^{+0.51}$	${23.24}_{-0.12}^{+0.12}$
	Mini-BAL	${19.87}_{-0.05}^{+0.08}$	${20.13}_{-0.05}^{+0.05}$	${20.21}_{-0.08}^{+0.08}$	${20.12}_{-0.09}^{+0.12}$	${20.66}_{-0.06}^{+0.07}$
$\mathrm{log}({\rm{R}})(\mathrm{pc})$	High	N/A b	⋯	⋯	⋯	${0.55}_{-0.13}^{+0.10}$
	Medium	N/A	$-{1.58}_{-0.15}^{+0.15}$	$-{1.60}_{-0.24}^{+0.21}$	$-{1.12}_{-0.28}^{+0.31}$	< −1.16
	Low (1)	N/A	$-{1.59}_{-0.13}^{+0.12}$	$-{1.54}_{-0.20}^{+0.17}$	$-{1.20}_{-0.24}^{+0.25}$	$-{0.99}_{-0.08}^{+0.07}$
	Low (2)	N/A	$-{1.29}_{-0.11}^{+0.11}$	$-{1.25}_{-0.19}^{+0.15}$	$-{0.93}_{-0.41}^{+0.38}$	$-{0.83}_{-0.04}^{+0.04}$
	Mini-BAL	N/A	N/A	N/A	N/A	N/A
$\mathrm{log}(\dot{{\rm{M}}})({{\rm{M}}}_{\odot }\,{\mathrm{yr}}^{-1})$	BAL	$-{0.64}_{-0.10}^{+0.14}$	$-{0.81}_{-0.14}^{+0.15}$	$-{0.97}_{-0.22}^{+0.20}$	${0.10}_{-0.30}^{+0.30}$	${0.73}_{-0.11}^{+0.11}$
$\mathrm{log}({{\rm{L}}}_{{KE}})(\mathrm{erg}\,{{\rm{s}}}^{-1})$	BAL	${42.92}_{-0.11}^{+0.14}$	${42.14}_{-0.14}^{+0.15}$	${41.97}_{-0.21}^{+0.19}$	${43.14}_{-0.29}^{+0.26}$	${43.80}_{-0.11}^{+0.11}$
$\mathrm{log}(\dot{{\rm{P}}})(\mathrm{dyne})$	BAL	${34.18}_{-0.11}^{+0.14}$	${33.69}_{-0.14}^{+0.15}$	${33.53}_{-0.21}^{+0.19}$	${34.64}_{-0.29}^{+0.28}$	${35.28}_{-0.11}^{+0.11}$

Note. In each velocity section, values are the mean of the bins in each velocity group where the parameters were well sampled. Upper and lower limits represent the 95% confidence interval sampled from the distributions of each parameter. Velocity ranges given are approximate for each epoch.

^a $\mathrm{log}({{\rm{N}}}_{H})+\mathrm{log}(\tfrac{1}{1+{\rm{a}}})$. ^b N/A is given for regions and parameters where the density could not be constrained due to a lack of C iii* absorption.

Download table as:  ASCIITypeset image

5.1.1. Radius

The radius of the BAL is calculated from the density and ionization parameter values from the SimBAL fits, and Q, the number of ionizing photons per second, using the ionization parameter equation defined in Section 3.1. The value of Q is found by scaling a standard quasar SED (the soft quasar SED used to create the SimBAL grids; Hamann et al. 2011) to best match the UV SED and photometry available for each epoch, and then integrating the photon flux for energies larger than 13.6 eV. Q is estimated to be $\mathrm{log}(Q)=53.8-54.1$ (phot s⁻¹); the range in values is due to the differences in the normalization for each epoch (see Figures 2, 3).

The radius of the medium- and low-velocity gas has lower uncertainties because the density is better constrained by the presence of C iii* absorption. To estimate the radius of the high-velocity gas with no independent C iii* constraints, we assumed the same density values as the lower-velocity gas for each epoch. The lowest (mini-BAL) velocity absorbers have no independent density constraints, but are likely to have larger radius values than the BAL gas because the mini-BAL absorbing gas completely covers both the continuum and emission-line regions.

The radius per velocity bin is presented in Figure 9. The density is measured per epoch and ionization parameter is measured per velocity group. The low- and medium-velocity BAL gas has an average value across all epochs of 0.07 pc. The much larger radius calculated for the high-velocity component (3.5 pc) may be an artifact of the assumption that both components share the same density when they have much different values of ionization parameter. If the high-velocity component has a higher density, the radius would be lower; a density value of $\mathrm{log}({{\rm{n}}}_{H})\sim 8$ would put the high-velocity gas at the same radius as the lower-velocity gas. No radius value is presented for the BAL in the 2010 epoch or for the mini-BAL in all epochs due to the lack of C iii* absorption (the line is not in the bandpass in 2010) to directly constrain the density.

Zoom In Zoom Out Reset image size

For reference, the dust sublimation radius is estimated to be 0.01 pc calculated using the equation ${R}_{\mathrm{sub}}=0.2{L}_{46}^{1/2}$ pc from Laor & Draine (1993) and a bolometric luminosity of 5.20 × 10⁴³erg s⁻¹ (Leighly et al. 2015). With a radius estimate of 0.07 pc from SimBAL models, this places the winds in the torus, consistent with Leighly et al. (2015) who found the radius to be in the vicinity of the torus with a distance of 0.17–1.47 pc based on photometric variability.

5.1.2. Column Density

5.1.3. Mass-loss Rate

5.1.4. Kinetic Luminosity

As the power-law partial-covering parameter is the most empirical and least intuitive of the SimBAL physical parameters, we calculated the effective covering fraction per bin from the partial-covering parameter. The effective covering fraction is the fraction of the emission region covered. Since power-law partial covering must, by definition, assume that no region is truly uncovered, to calculate the total effective covering fraction, we assume some cutoff in optical depth, below which we assume there is no contribution to total covering. Effective covering fraction is taken to be the fraction of the surface area covered by gas with an opacity greater than 0.5. The τ > 0.5 criteria is given in Arav et al. (2005) as a good cutoff to estimate the effective covering fraction. Optical depth as a function of velocity in the power-law partial-covering case is defined as $\tau ={\tau }_{{\max }}{x}^{a}$ where ${\tau }_{\max }=2.654\times {10}^{-15}f\lambda {N}_{{ion}}(v)$ (Savage & Sembach 1991) and f is the oscillator strength of the transition, λ is the wavelength (Å), and N_ion (v) is the ion column density (${\mathrm{cm}}^{-2}{\left(\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{-1}$). This calculation was done for the high-ionization lines of C IV, Si IV, N V, and P V, as well as the low-ionization lines C II and C iii* to examine changes in effective covering fraction as a function of velocity, epoch, and ion line as was done in Figure 16 of Leighly et al. (2019).

In Figure 10, we show the effective covering fraction per velocity bin as a function of transition and time from 2010–2017. These results show that the high-ionization lines cover more of the continuum than the low-ionization lines. The covering fraction was higher for the medium-velocity gas in 2010 and 2017 when the BALs were deepest and the high-velocity component was present. The effective covering fraction was lowest for the medium- and low-velocity gas in 2015 when the high-velocity BAL component was absent and the BALs were the shallowest. At the lowest velocities, the effective covering fraction does not change significantly from 2013–2017 and even appears to increase in 2015, although these bins overlap with emission lines so small changes between epochs could reflect changes in the emission lines and not a change in absorption. For the bins around 8000 km s⁻¹ there is a change in covering fraction of Si IV from <0.1 in 2015 to 0.4 in 2017, whereas we see almost no change in Si IV for the low-velocity bins between the two epochs. Unlike previous work with WPVS 007, we find absorption at low velocities in 2015 due to our finding that the line emission is uncovered. The lowest-velocity gas seemed to disappear from C IV and P V in 2010, but the Si IV component remained well covered. Overall, we observe small changes in effective covering align with small changes in BAL depth for a particular velocity bin and large changes in effective covering align with regions of the BAL where larger changes were observed.

Zoom In Zoom Out Reset image size

The mini-BAL parameters are the same within errors between epochs for $\mathrm{log}(U)$, $\mathrm{log}({N}_{{\rm{H}}})-\mathrm{log}(U)$, $\mathrm{log}({N}_{{\rm{H}}})$, and $\mathrm{log}(a)$, consistent with the observed lack of variation in mini-BAL depth between epochs. This gas shows lower ionization parameter ($\mathrm{log}(U)$ average of −0.79 across all epochs) compared to the BAL gas ($\mathrm{log}(U)$ average of 1.27 across all epochs and velocities), lower column density (average $\mathrm{log}({N}_{{\rm{H}}})$ of 20.6 cm⁻²) compared to the total column density of the BAL (average $\mathrm{log}({N}_{{\rm{H}}})$ of 22.6 cm⁻²), and higher covering fraction (average $\mathrm{log}(a)$ of 0.71 for the mini-BAL) than the BAL (average $\mathrm{log}(a)$ of 1.5 across all epochs and velocities). We determined that the BAL does not cover the line-emitting region, whereas the mini-BAL fully covers the emission-line region. The difference in covering of the line and continuum from the BAL gas, together with the large covering fraction, suggests that the mini-BAL is an outflow separate from the BAL gas, existing on larger scales and likely farther away (although this cannot be confirmed without the density constraint needed to calculate R). An estimate for the ratio of the radius of the mini-BAL outflow to the radius of the BAL outflow for different epochs, assuming that the mini-BAL has the same density of the BAL for a particular epoch, shows that the mini-BAL could be around 10 times the distance of the BAL outflow.

With our new multiepoch SimBAL analysis with simultaneous fits of the continuum and several absorption lines, we are in a position to better constrain the physical properties of the outflow, and so we compare this work to previous studies of WPVS 007. From the SimBAL fit results, we found that the ionization parameter of WPVS 007 is high, up to $\mathrm{log}(U)$ = 1.88 (as seen in 2015). This ionization parameter is significantly higher than previous constraints on the ionization parameter from studies of this object. For example, Li et al. (2019) gave a maximum value for $\mathrm{log}(U)$ of 0.0.

The SimBAL analysis gives us confidence in our results, as we simultaneously fit multiple lines as a function of velocity. The strong, high-ionization lines of P V and N V in particular are powerful diagnostics. The model fits presented here show strong N V and P V BALs, with a P V apparent optical depth almost matching that of C IV at the lowest velocities. These lines together imply a large ionization parameter and a large column density (e.g., Hamann 1998) because of the low cosmic abundance of P V (P/C = 0.00093; Grevesse et al. 2007). Specifically, we find that the ionization parameter varies from $\mathrm{log}(U)$ = −0.28 to $\mathrm{log}(U)$ = 1.88 for bins that cover the P V BAL. Our value is consistent with the previous analysis of the 2003 FUSE observations by Leighly et al. (2009) that gave the lower limit of $\mathrm{log}(U)$≳0.0 from the strong, almost saturated P V absorption line and saturated N V line. In addition, the assumption of Li et al. (2019) that the N V BAL is weak is not consistent with our spectral-fitting results. Assuming complete covering of the continuum and emission lines leads to a significant underestimate of the amount of absorption from the N V line. BALs not fully covering emission lines have been observed before (e.g., Arav et al. 1999; Borguet et al. 2012; Choi et al. 2022b). With a robust determination of the ionization parameter and the density constraints from the C iii* line, we used the definition of ionization parameter to calculate the radius of the outflow over the velocity bins where the C iii* feature is present. The average radius across velocity regions (assuming the density is uniform) and epochs (excluding the high-velocity gas in 2017 and the 2010 epoch) is 0.07 pc. The initial launching radius of the outflow can be estimated from ionization and kinematic arguments (based on assuming the Keplerian velocity at the launch radius and terminal velocity of the outflow are comparable), and gives much smaller values of ∼10⁻⁴ pc (Leighly et al. 2009; Li et al. 2019). Alternately, Leighly et al. (2015) gave an estimate for the outflow radius of 0.17–1.47 pc based on the variability, with a torus radius of R_{τ K} ∼ 0.036 pc. Our value for the outflow is thus also consistent with being found in the vicinity of the torus.

Understanding the physical reason behind BAL variability offers promise for understanding more about the geometry and structure in quasar outflows. With the use of SimBAL, we have the unique opportunity to study the difference between the physical properties of the gas over time. To date, many studies of variability use an observed change in equivalent width of the C IV line (or C IV and Si IV) to inform the understanding of the cause of the variability in the BAL (e.g., Capellupo et al. 2011, 2012; Filiz et al. 2012; Capellupo et al. 2013; Filiz et al. 2013, 2014). However, studies such as Lundgren et al. (2007) have emphasized the importance of having lines other than C IV (which is often saturated) when studying variability.

In the case of WPVS 007, we find this to be true; had we studied only the C IV line, we would not have been able to distinguish between a change in ionization parameter, column density, or covering fraction. In the multiepoch fits in Section 4.1, all tested scenarios were able to fit the C IV line for the 2013 and 2017 December epochs (see Figure 7). As mentioned above, the presence of P V in the spectrum, given the low abundance of phosphorus compared to carbon and its high ionization potential, is an indicator of high column density and high ionization parameter (e.g., Leighly et al. 2009; Capellupo et al. 2014; Leighly et al. 2018). In this parameter regime, the C IV absorption line will always be highly saturated (Borguet et al. 2012), and therefore large changes in some physical properties of the outflow may show little effect in the C IV EW. In particular, saturated lines—as C IV is whenever a P V BAL is present—are insensitive to changes in ionization parameter (McGraw et al. 2017; Vivek 2019). For WPVS 007, it was only with the combined simultaneous fit of several lines, including P V, N V, Si IV, and C IV, that a change in covering fraction was revealed to be the driving source of variability between epochs. From Figure 7, the models where the ionization parameter or column density were allowed to vary could not produce the full depth of the P V BAL. C iii* also was critical in the analysis as the only density-sensitive line. Like P V, C iii* is not saturated, but is also not present at all velocities where C IV is observed.

Previous studies of BAL variability have linked the cause of variability to one of two causes: (1) a change in ionization of the gas caused by a change in strength or shape of the ionizing continuum, or (2) a change in absorption caused by gas moving transverse across our line of sight with the quasar, the "cloud-crossing" scenario. The origin of the cause of variability is frequently used to determine the physical conditions of the gas. If a change in ionization is causing the observed variability, then the timescale between observations can be used to estimate the electron density and then the radius (e.g., Barlow 1993). If variability is caused by cloud crossing, the timescale and degree of change in the depth of the BAL can be used to estimate the crossing speed and therefore the radius (assuming Keplerian motions; e.g., Capellupo et al. 2013). Using SimBAL, the constraints on ionization parameter and density that allow us to calculate the radius are obtained directly from the MCMC fits.

Overall, it is unlikely that a change in continuum flux is driving the observed changes in the BALs. A changing ionization parameter causing the observed changes in absorption would indicate that the variability is likely being caused by a changing continuum. We found no support for a change in ionization parameter driving the variability in the BALs as discussed above through our use of simultaneous fitting as well as poor correlation between the change in ionization parameter and change in depth of the BALs. This interpretation of the variability differs from that proposed by Li et al. (2019). Li et al. (2019) examined patterns between the WISE infrared photometry, the Swift photometry, and the EWs of the C IV and Si IV lines to find that the cause of the variability is due to a change in the ionizing continuum instead of a change in a line-of-sight absorber. They were able to reproduce the observed Swift light curve using stochastic changes to the emission from the disk. They further calculate the concordance index (absorption-line EWs are assigned a +1 if the EW and continuum change in the same direction and −1 for the opposite case) to compare the stochastic changes in the disk to the variation in the UV Si IV and C IV BALs. In comparing the UV luminosity and BAL EW for the same five epochs of spectroscopy used in this analysis, Li et al. (2019) found a positive concordance index (the luminosity and BAL EW vary in the same direction) in all cases for C IV and most cases for Si IV. This method has been used to link a change in continuum with a change in absorber optical depth before (e.g., Wang et al. 2015). We observe the same stochastic changes in the photometry (see Figures 1 and 3) due to variation in the accretion disk, but disagree that these stochastic changes are the reason for the observed variation in BAL EWs. Using SimBAL, we are able to test these underlying assumptions that a changing EW is controlled by a change in the ionizing continuum. Given the saturation of the C IV, the EW of this line is not sensitive to changes in ionization state. Moreover, the high-ionization parameter required by the presence of P V absorption means that the concordance scenario proposed by Li et al. (2019) predicts a decrease of C IV EW with an increase in luminosity contrary to what has been observed.

Previous studies have linked coordinated variability across all velocities of a BAL to a change in ionization state, whereas a change in a portion of the BAL is often linked to a change in covering fraction (e.g., Capellupo et al. 2013). Although coordinated variability fits more naturally to a change in ionization parameter (a change in continuum would globally affect the outflow), this is not the case for saturated BALs, where covering fraction controls the apparent opacity of saturated lines (see Hall et al. 2011; Capellupo et al. 2014). Much like the case of Q1413 + 1143 observed by Capellupo et al. (2014), WPVS 007 shows a saturated C IV line with variability in C IV, Si IV, and P V, which indicates a change in ionization state of the gas cannot cause the variability, even if there are apparently coordinated changes across the BALs, because covering fraction is controlling the depth. Furthermore, though variations across the BALs appear to vary together, the changes across the BALs for WPVS 007 are not as coordinated as a function of velocity as they initially appear. Although the overall depths of the lines change, there are smaller changes along the BAL. As a function of velocity, not all components of the C IV line vary uniformly (e.g., the depth at the blue edge of the C IV line is changing separately from the red edge). For example, in 2015, most of the BALs seem to have gotten weaker in tandem. However, SimBAL analysis (which considered separately the high- and low-velocity BAL components) identified more absorption at low velocities in 2015 than previously modeled. Once the emission-line covering fraction was accounted for (with the low-velocity BAL gas not covering the emission lines), the low-velocity gas component did not disappear in 2015 as previously thought.

From fitting the continuum emission, we found that there is no significant change to either the power-law component or emission-line parameters between 2013 and 2017, but there is a significant change to the BAL absorption between these epochs. Furthermore, the continuum is different between 2015 and 2017, with the 2015 continuum being significantly redder, but the ionization parameter of the BAL gas between 2015 and 2017 is consistent. Although the Swift photometry shows variation epoch-to-epoch, these photons are longer in wavelength than the continuum photons that would be ionizing the gas in the outflow. We measured the continuum flux levels during the 2013 and 2017 HST observations by sampling the median flux in a line-free region of the continuum (1450 Å). We observe no change within errors between 2013 and 2017, while at the same time we observe a dramatic change in the BAL lines. The 2010 and 2015 epochs show the strongest and weakest continuum flux levels, respectively, consistent with the Swift photometry (see Figure 3). This indicates that there are changes in the continuum over time, but those changes are not driving the BAL variability between all epochs. For WPVS 007, we can therefore reject scenario (1), a cause in ionization parameter, and now consider the cloud-crossing scenario.

In previous variability studies, a change in covering fraction has been synonymous with the cloud-crossing scenario, where an absorber passes across the continuum-emitting region along the line of sight (Hamann et al. 2008; Gibson et al. 2010; Hall et al. 2011; Vivek et al. 2012a; Filiz et al. 2012; Capellupo et al. 2013; McGraw et al. 2015; Vivek et al. 2016; McGraw et al. 2017). In particular, Capellupo et al. (2014) found a case of a P V BAL showing variability in the P V, C IV, and Si IV lines and determined the cause of the variability to be cloud crossing. A change in absorption-line structure due to cloud crossing would produce a change in partial-covering parameter, $\mathrm{log}(a)$, and possibly also a change in column density in SimBAL. Although we have ruled out a change in column density as the main driver of the variability, the results of the regression analysis discussed in Section 4.2.1 imply that column density may also be a contributing factor determining the change in depth of the BALs. The cloud-crossing scenario for the case of WPVS 007 is shown in the middle panel of Figure 11.

Zoom In Zoom Out Reset image size

We can calculate the expected time for the BAL to completely cross the continuum, assuming Keplerian motion of the gas around the black hole. Using a radius of 0.07 pc, the speed of the gas is about 500 km s⁻¹, from ${v}_{\mathrm{cross}}=\sqrt{{GM}/R}$ where G is the gravitational constant, M_BH is the mass of the black hole, and R_BAL is the outflow radius. Given the size of the accretion disk at 1550 Å is estimated to be R₁₅₅₀ = (3.4 − 7.9) × 10⁻⁴ pc (Leighly et al. 2015), it would take between ∼10 months and ∼1.6 yr to cross the continuum at R. This crossing time is shorter than the time between the luminosity peak in 2010 and the occultation event in 2015. We can use the cloud-crossing method outlined in Capellupo et al. (2013) to calculate the crossing speed from the change in depth of the BAL, assuming the absorber and continuum are uniform and circular. Considering all bins and epochs, the crossing speed would range from 6.6–423 km s⁻¹, which leads to an outflow radius of 0.1–400 pc. At the low end, this is only marginally higher than the 0.07 pc radius determined using SimBAL, making the cloud-crossing scenario a possible cause of the variability observed in WPVS 007.

Leighly et al. (2015) speculated that the outflowing gas may be ablated material driving from the torus the same way "wind ablates spray from the crest of a wave." The outflow would therefore be the upper edge of the torus; this is consistent with the constraints on R_BAL from SimBAL and R_torus from Leighly et al. (2015). It is widely accepted that the torus must be clumpy (e.g., Netzer 2015). Given this picture, Leighly et al. (2015) proposed that the variability is caused by the changing scale height of the torus that on occasion more completely eclipses our line of sight. The changing torus scale-height model assumes that our line of sight skims the edge of the torus, and as the torus rotates, our view of the quasar is either relatively unobscured when the torus has a low scale height (in the case of 2010 where the BAL was strong, the photometry blue, and the emission lines broad), or highly obscured when the torus has a high scale height (in the case of 2015 where the BAL was weak, the photometry red, and the emission lines narrow). The outflow is clumpy and embedded in a larger region of gas and dust that controls the variability on longer timescales. A changing covering fraction would be consistent with the Leighly et al. (2015) picture. The changing torus scenario is shown in Figure 11 and in Figure 5 of Leighly et al.(2015).

Alternatively, if the outflow is located beyond the location of the torus, when the scale height of the torus is low (2010, 2013, 2017), it does not block any of the continuum, and therefore the angular size of the continuum is larger from the point of view of the outflow. (Note that we are assuming that the torus is completely opaque in this picture). In contrast, a high torus scale height (2015) blocks part of the continuum. In this case, the changing covering fraction results from the changing angular size of the continuum, and not a cloud crossing the observer's line of sight. When the angular size of the continuum as viewed from the outflow is larger, more of the outflow can cover the continuum, and so the covering fraction is larger; the changing scale-height model is the geometry that allows this to happen.

Alternatively, a change in covering fraction without a change in column density could be consistent with the formation and dissipation of clumps of material along the line of sight. Due to the way that SimBAL uses power-law partial covering, a change in power-law partial covering could be interpreted as a change in how diffuse clumps are along our line of sight (see Figures 12 and 20 in Leighly et al. 2019). We have used this idea from Leighly et al. (2019) to construct a diagram of what could be occurring in WPVS 007 (see the bottom panel of Figure 11). In this picture, the smaller, dense portions of clumps will have high opacity and lower covering fraction. However, the steep opacity profile would mean low opacity for other (larger) regions within the wind; we propose that this could be the case for the wind in 2015 when the covering fraction of the high- and medium-velocity gas is low. With a shallow opacity profile across clouds, we would see increased optical depth and therefore a higher covering fraction across all velocities (this would be the case for 2017, as we demonstrate in Figure 11).

We calculated the dissipation timescale of a clump of material, assuming that the high-velocity, medium-velocity, and low-velocity gas all represent distinct single "clouds" within the outflow. The diameters of the clouds are estimated using R ∼ N_H /n_H , and the dissipation timescale is assumed to be t ∼ R/v, where v is the maximum velocity of the gas (e.g., Hamann et al. 2013). We found that the dissipation timescale has a large variation between all epochs and velocity groups. The shortest dissipation timescales occur for the high-velocity gas in 2017, where dissipation occurs on the order of 3 hr. The medium-velocity gas across all epochs shows dissipation timescales on the order of hours for the 2013 observations, which is inconsistent with the lack of observed variability between June and December observations. The longest timescales we find are for the low-velocity gas in 2017 with dissipation timescales of almost 2 yr. Considering that we see little variation for the low-velocity gas from 2013–2017, this is also inconsistent. The assumption that entire portions of the outflow consist of single large clumps of gas is very limiting and likely not physical. SimBAL itself uses a partial-covering approach that assumes that the wind is made up of an ensemble of small clouds. However, the assumption that the clouds are small would further reduce the dissipation time of the clumps. MHD simulations of disk winds have shown evidence for the formation and dissipation of high-density knots that form from instabilities in the wind every ∼3 yr (Proga et al. 2000); however, these simulations were for line-driven accretion-disk winds in a much larger 10⁸ M_⊙ system. The outflow of WPVS 007 is highly unusual in its dramatic variability and may make an ideal candidate for evaluating future models of AGN outflows.