Clumpy AGN Outflows due to Thermal Instability

Thermal instability (TI; Field 1965) was long ago recognized as a viable mechanism for producing multiple phases in active galactic nuclei (AGN) winds (e.g., Davidson & Netzer 1979; Krolik & Vrtilek 1984; Shlosman et al. 1985). Such winds are observed, for example, in some Seyfert galaxies, where the ultraviolet (UV) and X-ray absorbers have similar velocities, strongly suggesting that very different ions are nearly cospatial and therefore that different temperature regions coexist (e.g., Shields & Hamann 1997; Crenshaw et al. 1999; Gabel et al. 2003; Longinotti et al. 2013; Ebrero et al. 2016; Fu et al. 2017; Mehdipour et al. 2017, and references therein).

While TI is well understood in a local approximation (e.g., Balbus 1995; Waters & Proga 2019), it has proven challenging to quantitatively model clump formation in a dynamic flow. Only in recent years has it become clear how the in situ production of multiphase gas can be triggered using global time-dependent hydrodynamical simulations—and only in the context of accretion flows (Barai et al. 2012; Gaspari et al. 2013; Mościbrodzka & Proga 2013, hereafter MP13; Takeuchi et al. 2013) or stratified atmospheres (e.g., McCourt et al. 2012; Sharma et al. 2012). In an outflow regime, previous work has focused on highlighting the importance of considering the effects of clumpiness, but only on a qualitative basis (e.g., Nayakshin 2014; Elvis 2017).

Here, we present the first global simulations of an outflow that is multiphase due to TI. Specifically, we show that there exists a range of the so-called hydrodynamic escape parameter (HEP), small yet relevant, for which the outflow can develop regions with significant over- and under-densities and that the over-densities can survive their acceleration over a relatively large distance. We identify the physical effects that contribute to making the HEP range small, and thus explain why a clumpy outflow has not been identified in any previously published results from the simulations of thermally driven outflows, neither in 1D nor 2D simulations (e.g., Woods et al. 1996; Proga & Kallman 2002; Luketic et al. 2010; Higginbottom & Proga 2015; Dyda & Dannen 2017; Waters & Proga 2018). We detail our numerical methods in Section 2. We present the results from our calculations in Section 3 and a discussion in Section 4.

We employ the magnetohydrodynamics code Athena++ J. M. Stone et al. (2020, in preparation) to solve the equations of non-adiabatic gas dynamics

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{\boldsymbol{v}}+{\boldsymbol{P}})=-\rho {\rm{\nabla }}{\rm{\Phi }}+{{\boldsymbol{F}}}_{\mathrm{rad}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot [(E+P){\boldsymbol{v}}]=-\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}{\rm{\Phi }}-\rho { \mathcal L }+{\boldsymbol{v}}\cdot {{\boldsymbol{F}}}_{\mathrm{rad}}\end{eqnarray} \tag{ 3 }$$

where ρ is the fluid density, ${\boldsymbol{v}}$ is the fluid velocity, ${\boldsymbol{P}}\,=\,p\,{\boldsymbol{I}}$ with p the gas pressure and ${\boldsymbol{I}}$ the unit tensor, Φ = −GM_BH/r is the gravitational potential due to a black hole with mass M_BH, $E=\rho { \mathcal E }+1/2\rho | {\boldsymbol{v}}{| }^{2}$ is the total energy with ${ \mathcal E }={\left(\gamma -1\right)}^{-1}p/\rho$ the gas internal energy, ${{\boldsymbol{F}}}^{\mathrm{rad}}={F}^{\mathrm{rad}}\,\hat{{\boldsymbol{r}}}$ is the radiation force, and ${ \mathcal L }$ is the net cooling rate. All of our calculations assume γ = 5/3. Although we have computed models with the radiation force due to both electron scattering and line driving, in this Letter we present only models with the force due to electron scattering, i.e., F_rad = GM_BHρ Γ/r², where Γ = L/L_Edd is the Eddington fraction. Irradiation is thus assumed to be due to a point source, as is appropriate when considering parsec scales in relation to the X-ray coronae and the UV-emitting regions of AGN disks.

For the unobscured AGN spectral energy distribution (SED) of NGC 5548 from Mehdipour et al. (2015), shown in the top panel of Figure 1, D17 have tabulated the net cooling rate ${ \mathcal L }={ \mathcal L }(T,\xi )$, which is function of the gas temperature, T, and the ionization parameter, defined as ξ = L_X /nr² where L_X is the luminosity integrated between 1–1000 Ry, and $n={\left(\mu {m}_{p}\right)}^{-1}\rho$ is the gas number density (with m_p the proton mass; we set μ = 0.6 in this work). The ξ dependence on distance plays an important role in determining the thermal stability of the flow, as discussed in detail in Section 4. Also, the ratio between the hard and soft X-ray energy bands is an important characteristic affecting properties of TI (e.g., Kallman & McCray 1982; Krolik 1999; Mehdipour et al. 2015; Dyda & Dannen 2017).

Zoom In Zoom Out Reset image size

The solid line in the bottom panel of Figure 1 is the S-curve corresponding to this SED, i.e., the contour where ${ \mathcal L }(T,{\rm{\Xi }})=0$, where ${\rm{\Xi }}\equiv {p}_{\mathrm{rad}}/p=\xi /(4\pi {{ck}}_{{\rm{B}}}T)$ is the pressure ionization parameter, with p_rad the radiation pressure (equal to F_X/c in our models). For the adopted SED, gas is unstable by the isobaric criterion for local TI in the two zones where the slope of the S-curve is negative (these are the locations where Field's criterion for TI, ${[\partial { \mathcal L }/\partial T]}_{p}\lt 0$, is satisfied). Gas is actually thermally unstable everywhere left of the dashed line (hereafter referred to as the Balbus contour), as this region of parameter space satisfies Balbus' generalized criterion for TI, ${[\partial ({ \mathcal L }/T)/\partial T]}_{p}\lt 0$ (Balbus 1986). Points corresponding to the maximum value of Ξ on the cold stable branch of the S-curve (i.e., for $\mathrm{log}\,{{\rm{\Xi }}}_{{\rm{c}},\max }=1.10$) are dynamically the most significant, as they mark the entry into a cloud formation zone and dramatic changes in the flow profiles can occur there.

We present the results of 1D models run at both medium and high resolution, all using a uniform grid spacing $[{\rm{\Delta }}r\approx ({r}_{\mathrm{out}}-{r}_{0})/{N}_{r}$]. The standard resolution runs (Models A-1x, B-1x, C-1x, and D-1x) have N_r = 552, chosen such that gradient scale heights ${\lambda }_{q}\equiv q/| {dq}/{dr}|$ (where q is any hydrodynamic variable) are adequately resolved with λ_q/Δr ≈ 3 for r ≲ 1.1 r₀ and λ_q/Δr > 10 at all other points in the wind except in the vicinity of the location where ${\rm{\Xi }}={{\rm{\Xi }}}_{c,\max }$, where this ratio dips to about 1. The high-resolution runs (Models A-8x, B-8x, C-8x, and D-8x) have N_r = 8 × 552, nearly an order of magnitude higher resolution everywhere. The other relevant length scale to resolve is that of the clumps, which have a size on the order of the local cooling length scale ${\lambda }_{\mathrm{cool}}\equiv {c}_{s}\,{t}_{\mathrm{cool}}$, where ${t}_{\mathrm{cool}}={ \mathcal E }{\left(n{ \mathcal C }/{m}_{p}\right)}^{-1}$ is the cooling time while ${ \mathcal C }$ is the cooling rate in units of $\mathrm{erg}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$. For γ = 5/3,

$$\begin{eqnarray}&&{\lambda }_{\mathrm{cool}}\approx 3\times {10}^{16}\,{T}_{5}^{3/2}\,{n}_{3}^{-1}\,{{ \mathcal C }}_{23}^{-1}\ \mathrm{cm},\end{eqnarray} \tag{ 4 }$$

where ${T}_{5}=T/{10}^{5}\,{\rm{K}}$, ${n}_{3}=n/{10}^{3}\,{\mathrm{cm}}^{-3}$, ${{ \mathcal C }}_{23}={ \mathcal C }/{10}^{-23}\,\mathrm{erg}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$, and ${T}_{5}={n}_{3}={C}_{23}=1$ are characteristic values of the clumps. With Δr ≈ 2 × 10¹⁵ cm for our "-1x" runs, this length scale is adequately resolved.

We apply outflowing boundary conditions at the inner and outer radii with the density fixed at the innermost active grid point to the value of ${\rho }_{0}$ in Table 1. The initial conditions are a simple expanding atmosphere with $\rho ={\rho }_{0}\,{\left(r/{r}_{0}\right)}^{-2}$ and a β-law velocity profile, $v={v}_{\mathrm{esc}}\sqrt{1-{r}_{0}/r}$, where the "0" subscript is used to denote values at the inner radius of the computational domain, r₀. The pressure profile is set according to the temperature, chosen so that the gas lies along the S-curve, having a value Ξ₀ at r₀ (see Section 3).

Table 1.  Summary of Key Parameters and Results

Model	HEP	r0	ρ0	${t}_{\mathrm{sc},0}$	${t}_{\mathrm{cool}}/{t}_{\mathrm{sc}}$	$\langle v\rangle$	$\langle \dot{M}\rangle$	Comment
		(1018 cm)	$({10}^{-18}$ g cm−3)	(1012 s)		(107 cm s−1)	(1024 g s−1)	1x runs	(8x runs)
A	13.3	1.01	2.68	3.4	0.47 (0.48)	6.4 (6.0)	1.2 (1.2)	steady	(steady)
B	11.9	1.13	2.15	3.9	0.54 (0.49)	3.4 (3.5)	3.2 (3.3)	unsteady	(unsteady)
C	9.1	1.48	1.26	5.0	0.41 (0.42)	2.5 (2.3)	5.9 (6.0)	unsteady	(unsteady)
D	8.1	1.67	0.98	5.7	0.41 (0.41)	2.4 (2.2)	6.6 (6.7)	quasi-steady	(unsteady)

Note. HEP is the hydrodynamic escape parameter. The inner radius is derived from the choice of HEP as ${r}_{0}=\left(1-{\rm{\Gamma }}\right){{GM}}_{\mathrm{BH}}{\mathrm{HEP}}^{-1}{{\rm{c}}}_{s,0}^{-2}$ (see Equation (5)), while the density at the base follows from the definition of ξ as ${\rho }_{0}=\mu {m}_{p}{L}_{{\rm{X}}}{\xi }_{0}^{-1}{r}_{0}^{-2}$. The sound crossing time is defined as ${t}_{\mathrm{sc},0}=({r}_{\mathrm{out}}-{r}_{0})/{c}_{{\rm{s}},0}$. The ratio ${t}_{\mathrm{cool}}/{t}_{\mathrm{sc}}$ is given at the location where ${\rm{\Xi }}={{\rm{\Xi }}}_{{\rm{c}},\max }$. The average mass flux and velocity through r_out are shown as $\langle \dot{M}\rangle$ and $\langle v\rangle$. Comment columns denote the state of the flow at late times. Values/comments in parenthesis denote results for the 8x-resolution runs.

Download table as:  ASCIITypeset image

The 2D model uses a 256 × 256 grid in r and θ, with logarithmic spacing in r such that ${{dr}}_{i+1}=1.01\,{{dr}}_{i}$ and uniform spacing in θ from 0 to π. At the inner boundary we assumed a density profile $\rho ={\rho }_{0}[1+0.001\sin (2\theta )]$ to break spherical symmetry. We apply reflecting boundary conditions at 0 and π.

For a given M_BH, just three parameters govern our solutions: Γ, Ξ₀, and the HEP, which sets the strength of thermal driving. We define HEP as the ratio of effective gravitational potential and thermal energy at r₀,

$$\begin{eqnarray}&&\mathrm{HEP}=\displaystyle \frac{{{GM}}_{\mathrm{BH}}(1-{\rm{\Gamma }})}{{r}_{0}{c}_{s,0}^{2}}.\end{eqnarray} \tag{ 5 }$$

We only present models with Ξ₀ = 3, Γ = 0.3, and ${M}_{\mathrm{BH}}={10}^{6}\,{M}_{\odot }$, as HEP is the main governing parameter. We note, however, that the choice of Ξ₀ is not unimportant, e.g., selecting one too large can cause the flow to miss relevant regions. The dependence of stable wind solutions on Γ for various SEDs was explored in 1D by D17. In Table 1, we list model parameters and summarize gross outflow properties.

After examining over 100 1D simulations, with the standard resolution that span this parameter space, we arrived at four qualitatively different 1D wind solutions that capture the general behavior seen across all runs. These four cases illustrate how the stability of the outflow depends on the HEP. We chose the values of HEP for our models A and D such that they closely bracket the parameter space leading to transonic, clumpy winds (here represented by models B and C). For each model, Figure 2 shows the models tracks on the phase diagram, while Figure 3 shows radial profiles of ρ, v, T, and Ξ. The phase diagrams reveal that all solutions pass through the lower TI zone, but only models B and C with the standard resolution actually trigger TI to become unsteady. Additionally, two general trends are evident: (i) the range of Ξ decreases with HEP and (ii) for large Ξ, the wind is not in thermal equilibrium; the wind temperature is nowhere near that of the stable Compton branch on the S-curve (the maximum temperature shown in these plots is more than an order of magnitude lower than T_C). The first trend is due to the fact that for thermal winds, the velocity is a decreasing function of HEP, and by mass conservation for radial flows, ξ ∝ v in a steady state (i.e., v ∝ 1/n r², as is ξ). Trend (ii) is due to adiabatic cooling (see D17), but note that the highest temperatures reached in all four cases occupy thermally stable regions of the (T–Ξ)-plane (namely, regions to the right of the Balbus contour). We now examine the dynamics of these solutions in detail to understand why Models B, C, and D-x8 are clumpy, while A and D-x1 are not.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Unsteady, Clumpy Wind Solutions

3.2. Smooth Wind Solutions Passing through a TI Zone

3.3. Clumpy Winds in 2D

We conclude the presentation of our results with one example of a 2D simulation, the counterpart to model B. The bottom panels in Figure 4 shows the temperature and density in cgs units, whereas the top panels show the relative difference between these quantities and their time-averaged values ($\langle T\rangle$ and $\langle \rho \rangle$) at a given location for a snapshot near the end of the simulation. Notice that the flow is very clumpy at $r\simeq 2\,{r}_{0}$ but appears quite smooth at larger radii where the density is overall smaller and clumps have begun to expand. This is apparent from the density profiles in our 1D simulations, although not altogether obvious. The maximum density/temperature contrast is never more than 10.

Zoom In Zoom Out Reset image size

Compared to model B in 1D, we found a significant scatter of the 2D wind profiles at a given radius, which reflects the loss of spherical symmetry. The amplitude of the scatter (i.e., the variability in the θ direction) is of the same order as the variation in the radial direction. This in turn shows that the perturbations grow to a similar level in both directions.

In this study, we identified a finite parameter space for which a thermally driven wind is clumpy. The origin of the clumps in our X-ray irradiated outflows is very similar to that of clumpy accretion flows studied by Barai et al. (2012) and MP13: quite simply, small perturbations are allowed to grow due to TI. However, we found that because the flow enters the TI zone near the cold phase, TI mainly serves to raise the temperature of perturbations. Therefore, the clumpiness is a result of the separation of heated layers of gas from the cold, dense layers near the base of the flow rather than the formation of dense clumps within a more tenuous background plasma. Importantly, this formation of heated layers from TI requires that a perturbation can stay in the TI zone long enough for it to become nonlinear and that stretching due to rapid radial acceleration does not stabilize the growth of such layers.

These two requirements represent necessary conditions for TI to operate in dynamical flows and are due to velocity gradients in the wind that are not accounted for in the classical theory of TI. We similarly identified another necessary condition that distinguishes dynamical TI from local TI, this one arising from pressure gradients: Ξ must decrease once gas enters a TI zone (see Section 3.1). Because Ξ is the ratio of two pressures, both of which are decreasing functions of radius, Ξ could be a non-monotonic function of radius. Specifically, the models here have ${\rm{\Xi }}\propto 1/({r}^{2}p)$, and even though p decreases with radius, this decrease could be slower than 1/r² and Ξ can decrease downstream.

To build intuition for when to expect a clumpy versus smooth outflow, we could just consider the geometric effects involved in the problem. For example, let us assume that the radiation flux scales as r^−q, where q is a constant. We then find the following trend: the more q is less than 2, the more solutions resemble those from D17, where the wind is steady and monotonic and Ξ shows no "back tracking." For 2.0 ≲ q ≲ 2.5, the wind is unsteady and the range of Ξ and T is reduced compared to the cases presented here. For q ≳ 2.5, the flux drops so fast that the heating is very weak. Consequently, Ξ and T are never large enough for the solution to even approach the TI region and as a result the flow is smooth.

We have explored many other test runs, e.g., with line driving turned on or using a different SED. We found that clumpy outflows develop as in the examples shown here but for somewhat different HEP. We will report on the results from these simulations in a follow-up paper.

Support for Program number HST-AR-14579.001-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This work also was supported by NASA under ATP grant NNX14AK44.