Time-Dependent AGN Disc Winds I - X-ray Irradiation

The basic equations for single fluid radiation hydrodynamics are

where $\rho$ is the fluid density, $\mathbf{v}$ the velocity, $\sf{P}$ a diagonal tensor with components $P$ the gas pressure. The total gas energy is $E=\frac{1}{2}\rho|\mathbf{v}|^{2}+\mathcal{E}$ where $\mathcal{E}=P/(\gamma-1)$ is the internal energy and $\gamma$ the gas constant. The radiation momentum and energy source terms are $\mathbf{G}$ and $G^{0}$ respectively and receive contributions from the X-ray (see §2.2) and UV bands (see §2.3). $\mathcal{L}$ is an optically thin radiative heating term (see §2.4) and $\mathbf{g}_{\rm{grav}}$ is the gravitational acceleration (see §2.5)

The temperature is $T=(\gamma-1)\mathcal{E}\mu m_{\rm{p}}/\rho k_{\rm{b}}$ where $\mu$ is the mean molecular weight and other symbols have their standard meaning.

The radiation source terms $\mathbf{G}$ and $cG^{0}$ are assumed to receive contributions from a UV and an X-ray band

One can alternatively think of these components as representing emission from a compact central source (the X-ray component) and emission from a disk that is extended in radius (the UV component). In principle, there could be X-rays from the disk and UV from the central source, but the approximation used here is that the ionizing radiation (X-rays) are centrally concentrated and the radiation that potentially drives outflow comes from an extended disk alone. We will therefore model the ionizing central source radiation (X-ray component) via the radiation hydrodynamics module in Athena++, assuming this component is only incident into the domain on the inner radial boundary. We will model the disk emission (UV component) as an optically thin radiation force computed using the modified CAK formalism.

This work is meant to be the first in a series of papers to explore the impact of ionizing radiation on the launching of line driven winds in AGN. In order to connect with earlier work, we have endeavored to match the assumptions and methodological approach of PK04 as nearly as possible. This includes the treatment of the modified CAK line driving prescription, the heating rate, boundary conditions, and gravity, which are described below. As PK04 note, their prescription was a conservative one in that they made a number of unfavorable assumptions (e.g. ignoring the X-ray contribution to the radiation force) that were not conducive to driving an outflow.

The key difference between this work and PK04 is that we leverage the radiation hydrodynamics module of Athena++ to evolve the X-ray contribution to the radiation field by solving the time dependent radiation transfer equation directly. In contrast, PK04 solved the radially attenuated X-ray radiation field with an equation of the form

where $\tau=\int\rho\kappa dr$ is the optical depth integrated along radial directed rays. We adopt a prescription that allows us to closely emulate this pure attenuation treatment of X-ray flux while also allowing us to include scattered and reradiated radiation, as described below.

The X-ray radiation field is treated by directly solving the gray (frequency averaged) time dependent radiation transport equation using the implicit implementation in Athena++ (Jiang, 2021). We refer the reader to Jiang (2021) for the specific formulation but the equation solved is equivalent to

where $I$ is the frequency integrated specific intensity and the source term

where $\kappa_{s}$ is the scattering opacity, $\kappa_{F}$ is the absorption contribution to the flux mean opacity, $\kappa_{P}$ the Planck mean, $\kappa_{E}$ the energy mean opacity and $\Gamma=\gamma(1-\mathbf{n}\cdot\mathbf{v}/c)$ is a frame transformation factor with $\gamma=1/\sqrt{1-v^{2}/c^{2}}$ the Lorentz factor. We use a subcript "c" to denote variables evaluated in the comoving frame, including

the angle averaged comoving frame mean intensity. The X-ray momentum and energy source terms are then

The form of equation (5) is designed so that in the limit $v\rightarrow 0$

For this general scheme outlined above, a common approximation for the opacities is to set $\kappa_{E}$ equal to $\kappa_{P}$ so that radiative equilibrium corresponds to $E_{r}=aT^{4}$. Similarly, $\kappa_{F}$ is often set so that the sum $\kappa_{F}+\kappa_{s}$ corresponds to Rosseland mean opacity, with $\kappa_{F}$ being the contribution from absorption opacity and $\kappa_{s}$ is set to the electron scattering opacity $\kappa_{\rm es}$, which is nearly constant in the Thomson limit. This scheme is usually a sensible choice when the radiation variables represent the total frequency integrated radiation field, but here we consider a case where the radiation hydrodynamics solver is only being used to model the X-ray (ionizing) component of the radiation field. Since X-rays alone do not determine the net heating/cooling or radiation force, we need to consider an alternative prescription that accounts for this complexity.

In our current setup, we choose to turn off the direct coupling of the X-ray radiation field with the gas by setting $\mathbf{G}_{\rm{X}}=cG^{0}_{\rm{X}}=0$. The impact of the radiation on heating/cooling the gas is instead accounted for via the optically thin prescription described in § 2.4. In principle, we could still retain the momentum coupling given the X-ray contribution to the radiation force, but we drop this term to be consistent with previous work. Hence, only the radiation force from the UV disc component contributes directly to the radiation force. Instead, the radiation transfer of the X-rays is included solely to determine the associated ionization parameter

where $n=\rho/\mu m_{p}$ is the gas number density with $m_{p}$ the proton mass and $\mu=0.6$ the mean molecular mass. Since both the heating prescription and the UV line driving force are sensitive to the value of $\xi$, the X-ray component still indirectly impacts the radiation force and radiative heating/cooling.

These considerations also motivate an alternative prescription for the opacities that accounts for the fact that some fraction of the absorbed X-ray radiation field may be thermalized and reemitted outside the X-ray (ionizing) band. Determining the fraction reradiated in the X-rays would require a detailed non-LTE photoionization calculation that is beyond the scope of this work (and current numerical capabilities, at least for fully coupled radiation hydrodynamics). We intend to explore more sophisticated approximations in future work but for this study we consider a restricted class of models where we parameterize reradiation as a constant fraction $f_{\sigma}$ of the absorbed X-ray radiation that is remitted in the X-ray band. We assume this reradiated radiation is isotropically emitted in the comoving frame.

We implement the prescription described above by setting $\kappa_{s}=\kappa_{\rm es}$, $\kappa_{F}=\kappa_{a}$, and $\kappa_{E}=\kappa_{a}(1-f_{\sigma})$. Here, $\kappa_{a}$ is an attenuation opacity that we can vary to see how the radiation field responds under different assumptions. We also set $\kappa_{P}=0$ to decouple the reemitted radiation from gas temperature, which is instead set by our heating prescription. Our source term then reduces to

With these assumptions $\kappa_{a}$ is a form of pure absorption opacity that when $f_{\sigma}=0$, simply attenuates the radiation field. When $f_{\sigma}>0$, a fraction $f_{\sigma}$ of this attenuated radiation is reradiated isotropically in the comoving frame. When $f_{\sigma}=1$ our treatment of absorption and reemission is equivalent to our treatment of scattering.

Finally, we assume a constant radiation flux at the inner radial boundary directed along outgoing radial rays. To mimic the effects of a very compact source at the origin, only the most nearly radial outgoing rays are given a non-zero intensity on the inner boundary. This is accomplished by using a spherical polar discretization of the unit sphere (see e.g. White et al., 2023) with a polar axis corresponds to a unit vector in the local $r$ direction and selecting the ring of ray with smallest polar angle. The intensity of the selected rays is chosen such that the total flux in X-rays is a fraction $f_{X}=0.10\Gamma_{D}$ of the disc Eddington fraction.

The UV flux is assumed to be time-independent and the wind is assumed to be optically thin to the UV continuum. The momentum transfer can be broken up into contributions from electron scattering and radiation pressure on spectral lines,

where the surface integral is over the disc, assumed to be radiating like a Shakura-Sunyaev disc (Shakura & Sunyaev, 1973) with a constant accretion rate $\dot{M}$. Even though mass is injected into the domain, there is no feedback on the assummed mass flow rate, which remains constant. Hence, the radiation at the disk surface is modeled as

where $r_{*}=6GM/c^{2}$, $x$ is the re-radiation factor ($x=1$ in this work), and $\Gamma_{D}=\dot{M}/\dot{M}_{\rm Edd}$, where

We assume an efficiency $\eta=1/12$ in this work. In this model, the spectrum is blackbody, with a temperature

The fraction of the intensity in the UV is then found by integrating the Planck function

where $\lambda_{1}=200$ Å and $\lambda_{2}=3200$ Å and $\sigma$ is the Stefan-Boltzmann constant.

The strength of the spectral lines, relative to electron scattering is quantified via the force multiplier Owocki et al. (1988)

with $\alpha=0.6$ the ratio of optically thick to optically thin lines,

the number of lines and $\tau_{\rm{max}}=t\eta_{\rm{max}}$ with the parameter

determining the maximum number of lines available and the optical depth parameter

with the thermal velocity of the gas $v_{\rm{th}}=4.2\times 10^{5}$ cm/s. The disc as a source of radiation is taken to extend outside the computational domain to $R_{d}=1500r_{g}$.

Finally, the work done by the UV radiation force is then

We assume a locally optically thin radiative heating and cooling due to Compton, X-ray photoionization, Bremsstrahlung and lines. For a 10keV Bremsstrahlung, Blondin (1994) found a heating function of the form

where the contributions in cgs units, $\left[\rm{erg\ cm^{-3}\ s^{-1}}\right]$, due to Compton, X-ray photoionization, Bremsstrahlung and lines are respectively

with the line temperature $T_{L}=1.3\times 10^{5}$ and $\delta=1$ the parameter that controls the efficiency of line cooling.

We treat the disc as a mass reservoir of fixed density and fixed radiation intensity. Thus, we are implicitly assuming that our disc boundary is at a height $h$ above the disc midplane. For radiation dominated discs (Shakura & Sunyaev, 1973) the photosphere should be roughly a few times the Eddington fraction in units of $r_{\rm{ISCO}}$. We take the shift to be

above the disc midplane. In this work we choose $f_{h}=0.85$, where $f_{h}=1$ would correspond to a displacement by one Shakura & Sunyaev scale height. In these shifted coordinates, the gravitational acceleration acquires additional non-radial components

where

and

Follwoing PK04 we consider a central black hole mass $M_{BH}=10^{8}M_{\odot}$, corresponding to an innermost stable circular orbit (ISCO) $r_{*}=6r_{g}=8.8\times 10^{13}\rm{cm}$. We report some results in terms of gravitational radii $r_{g}=GM/c^{2}=1.4\times 10^{13}\rm{cm}$ and the orbital period at the inner boundary $t_{0}=2\pi((60r_{g})^{3}/GM)^{1/2}=1.44\times 10^{6}s$

We impose inflow (outflow) boundary conditions at the inner (outer) radial boundaries and axis boundary conditions along the $\theta$ = 0 axis. We assume a reflection symmetry about the $\theta=\pi/2$ midplane. We use a vacuum boundary condition for the radiation along the disc midplane and outer radial boundaries and keep the X-ray flux fixed at the inner radial boundary. After every full time step we reset $\rho_{d}=10^{-8}\rm{gcm^{-3}}$, $v_{r}=0$ and $v_{\phi}=v_{K}=\sqrt{GM/r}$. We also impose that the vertical velocity component $v_{\theta}$ is unchanged due to resetting density.

We choose a domain size $n_{r}$ × $n_{\theta}$ = 96 × 140. The radial domain extends over the range $10\ r_{*}<r<500\ r_{*}$ with logarithmic spacing $dr_{i+1}/dr_{i}=1.05$. The polar angle range is $0<\theta<\pi/2$ and has logarithmic spacing $d\theta_{j+1}/d\theta{j}=0.938$ which ensures that we have sufficient resolution near the disc midplane to resolve the acceleration of the flow.

Initially, the cells along the disc are set to have $\rho=\rho_{d}$, $v_{r}=v_{\theta}=0$, $v_{\phi}=v_{K}$. In the rest of the domain $\rho=10^{-20}\ \rm{gcm^{-3}}$ and $v_{r}=v_{\theta}=v_{\phi}=0$. Everywhere the temperature is constant along vertical cylinders corresponding to the Shakura-Sunyaev disc temperature at the base given by Shakura & Sunyaev (1973).

We take the disc Eddington number $\Gamma_{D}=0.5$ and re-radiation factor $x=1$ (see eq 14). This choice of $\Gamma_{D}$ corresponds to an accretion rate $\dot{M}=1.3\;\rm M_{\odot}/yr$ for $M_{BH}=10^{8}M_{\odot}$. We assume the luminosity of the central source radiation in X-rays is 10% of the corresponding disk luminosity. We assume the central source does not emit in the UV. The UV radiation from the disc is assumed to extend all the way from the ISCO (effectively outside the simulation domain) with $r_{*}\leq R_{d}\leq 1500r_{*}$. A sphere of radius $r_{*}$, effectively the black hole and nearby corona, is assumed to be optically thick for purposes of shielding the wind from the backside of the disc. . We impose a density floor $\rho_{\rm{floor}}=10^{-22}\rm{gcm^{-3}}$ which adds matter to stay above this floor, while conserving momentum, if the density ever drops below it. In addition, we have a temperature floor at the same cylindrical radius.

For our fiducial model, we utilize a pure attenuation ($\kappa_{a}=\kappa_{\rm es}$ and $f_{\sigma}=0$) treatment to emulate equation (3). We first test our results by using an angular grid with a single radial ray. This single ray approximation is the closest we come to equation (3) with the primary difference being that $\tau$ is accumulated over the light crossing timescale to the current radius rather than constructed from the current instantaneous value of the column. The wind properties obtained from this single ray approximation were consistent with PK04 results.

Since we are ultimately interested in the effects of scattering and reemission, we perform our simulations with $n=256$ rays uniformly covering the unit sphere. We checked for convergence of our model using $n_{r}=128,512$ rays and concluded this was a sufficiently resolved angular resolution.

We find radial, time variable disc wind solutions, broadly in agreement with the results in PK04. In Fig. 1 we show (going counter-clockwise) density (and gas velocity vectors), temperature, poloidal velocity (and unit gas velocity), force multiplier (and radiation force vectors), ionization parameter (and unit UV flux vectors) and X-ray radiation energy density ( and X-ray flux) at $t=196t_{0}$. On long time-scales, the outflow has two components, a faster, less dense flow at $\theta\sim 50^{\circ}$ and a denser, slower component at $\theta\sim 75^{\circ}$.

The fast component has a density $\rho\sim 10^{-18}\rm{g\ cm^{-3}}$ and velocities as high as 11 000 km/s. The slow component is an order of magnitude denser with $\rho\sim 10^{-17}\rm{g\ cm^{-3}}$ and velocities of $v\sim 5\ 000$ km/s.

The wind is optically thick with $N_{H}\gtrsim 10^{24}$ for $\theta\lesssim 60^{\circ}$. We can understand the wind launching via the X-ray radiation energy density. In the wind, $E_{r}$ is low, ensuring the ionization is also sufficiently low and the force multiplier large enough for the line driving to overcome gravity. From the force multiplier plot, we see the radiation force is sufficiently strong ($M(t)>1$, indicated by the green contours) in the very base of the disc, close to the central region. Within this contour the ionization $\xi\lesssim 10^{4}$ is sufficiently low so as to not suppress the force multiplier. In the limit of an optically thin wind, $\tau_{\rm{max}}\rightarrow 0$, the maximum force multiplier $M_{\rm{max}}=2$ for $\log\xi=2.23$. For our Eddington parameter $\Gamma=0.5$, this is roughly the minimum force multiplier to launch a wind, so our best case scenario is wind launching in the orange colored contours of the ionization plot.

The flow is time variable and we see periods with stronger and weaker outflows. In the top panel of Fig 2 we show the total mass flux out the outer boundary in $0^{\circ}\leq\theta\leq 85^{\circ}$. We excised the gas close to the disc because the disc density is quite high relative to the wind and dominated the signal. We have identified two epochs representative phases of the wind, a “weak" state (shaded in blue) and a “high" state (shaded in red). For each epoch, we show the time-averaged wind fluxes exiting the outer boundary as a function of $\theta$ (middle panels). We plot the column density $N_{H}=\int ndr$, the momentum flux $\rho v_{r}$ and the mass flux density $dm_{\theta}=4\pi\rho v_{r}r^{2}d\theta$. The shading shows the one standard deviation variations of the respective quantities with respect to time. The lower panels show the time-averaged density and velocity field.

During the low state, mass leaves the domain only at $\theta\sim 60^{\circ}$. During the high state, the mass flux is dominated by the more radial component where the velocities reach $\gtrsim 11\ 000\ \rm{km/s}$ but there is also a slower, less dense component at higher latitudes. The density profiles show that even in the high state the flow is unsteady. We see evidence of a failed wind with matter accreting both at latitudes just above the slow part of the wind.

Close to the disc, where the gravitational force scales as $|\mathbf{g}_{z}|\approx GMz/r^{3}$ matter can be supported against gravity via the radiation force due to electron scattering, which is approximately constant near the photosphere. As the matter is launched, the relative contribution of gravity and electron scattering changes and further acceleration requires line driving to take over. If line driving is too weak, due to overionization for example, the wind will fail and the launched gas will fall back towards the central object. The wind launching is a highly dynamic process, with gas being launched, failing, and then being swept back up in the outflow.

In Table 1 we show a summary of our other wind models. We indicate the gas/radiation coupling parameters that were varied (absorption opacity $\kappa_{a}$, scattering opacity $\kappa_{s}$ and re-emission fraction $f_{\sigma}$) and the resulting wind properties including mass outflow rate $\dot{m}$, outflow velocity $v_{\rm{out}}$, inclination angle $\omega$ and time between high state outflows $\Delta t_{\rm{high}}$. The goal is to explore the effects of scattering and re-radiation on the fiducial attenuation model. As an initial condition, we use the wind from the fiducial run at $t=2\times 10^{8}\rm{s}\approx 140.5t_{0}$, where a wind has been well established. We list the main radiation variables varied, namely the absorption and scattering opacitites and the re-radiation fraction. For each set of parameters, the solution is described as either a stationary wind, an episodic wind or having no wind. In the case of an outflow, we list the time averaged mass outflow rate, the peak outflow velocity and the azimuthal position of this velocity peak. In the case of episodic winds, the wind parameters are expressed during the “high" state during which an outflow is occurring. For episodic winds we also list the average period between outflow maxima.

We next consider cases in the pure scattering limit, where $\kappa_{a}=\kappa_{P}=0$, $f_{\sigma}=0$ and $\kappa_{s}=\kappa_{es}$. The radiation field changes on timescales of $\sim 10\ t_{0}$ as radiation is scattered from the low density funnel into the higher density wind. In the wind launching region, the X-ray energy density increases by $\sim 4$ orders of magnitude, leading to a commensurate increase in the ionization parameter. This effectively turns off the line driving and the wind gradually fails as new gas does not launch and the already launched gas exits the simulation domain.

To facilitate wind launching, we try the alternate scenario where matter is allowed to enter the initially empty domain from the disc midplane, as we used to launch our fiducial run. Matter enters the simulation domain, forming a thin disc, but the X-ray energy density is even higher than when using our restart prescription. Our results suggest that launching AGN disc winds is not possible in a pure scattering regime.

We therefore turn our attention to a regime with both scattering and absorption, $\kappa_{a}=\kappa_{s}=\kappa_{es}$. Unlike our fiducial model, we find a quasi episodic wind. In the top panel of Fig 3 we plot the mass flux exiting the simulation domain as a function of time. We identify several outbursts and periods of quiescence in the outflow, and highlight one such low (high) state shaded in blue (red). For these respective states we plot the time-averaged values of the dynamical variables at the outer boundary, and the density and velocity fields. As with the fiducial model, the low state is characterized by a slower, less dense and higher inclination flow than the high state.

The quasi periodicity is due to the ionization structure of the launching region. After outburst, the wind is too ionized and gas accretes towards the central region. This causes the density near the inner radius to peak, and the ionization to decrease. This allows for a larger force multiplier and a new outflow is launched. The time scale of this process is thus set by a combination of the inflow time and outflow time of the accreting/outflowing gas.

Increasing the absorption opacity has the effect of making the wind stationary and strengthening the outflow. Increasing $\kappa_{a}=2\kappa_{es}$ causes the wind to become stationary, with a fourfold increase in the mass flux compared to the quasi episodic winds high state and a doubling of the peak velocity. Further increasing $\kappa_{a}=10\kappa_{es}$, the values explored in PSK00, yield winds with 15 times the mass flux and peak velocities of 20 000 km/s.

We consider models where a fraction of absorbed radiation is re-emitted in the X-rays with $\kappa_{E}=\kappa_{a}\left(1-f_{\sigma}\right)$ and $0<f_{\sigma}\leq 1$.

Firstly consider models with only an absorption opacity. With only 5% re-radiation, $f_{\sigma}=0.05$, we find a very weak, episodic wind. Qualitatively it resembles the $A1S1$ model, with slightly weaker $\sim 50\%$ lower mass flux but 50% greater peak outflow velocities. Increasing the re-radiation to 10%, the qualitative behaviour is unchanged and we again find a weak, episodic outflow with peak velocities of $\sim 5\ 000\rm{km/s}$. Unsurprisingly, allowing gas to re-radiate has a similar physical effect to the wind as scattering. The X-rays tend to isotropize, and penetrate into the wind acceleration region. This tends to weaken the wind as gas becomes over-ionized. Further increasing the re-radiation fraction eventually leads to the wind being fully suppressed, which we see for $f_{\sigma}=0.5$.

We now consider models with both scattering, and absorption opacity. We set the re-radiation to 50%. When $\kappa_{a}=\kappa_{s}=\kappa_{es}$ we find that no wind is launched. This is expected as the model with these parameters and no scattering failed to launch a wind, so scattering could serve only to further inhibit wind acceleration.

We therefore increase the absorption opacity to $\kappa_{a}=2\kappa_{es}$ and find that we again enter a regime of episodic outflows. The wind is both quantitatively and qualitatively similar to the model without re-emission, A2S1.