CLOUD FORMATION AND ACCELERATION IN A RADIATIVE ENVIRONMENT

The thermal stability of astrophysical gases has been extensively investigated, both analytically and numerically, in many different physical contexts, e.g., in star formation regions, in gas accreted by forming galaxies, and in interstellar and intergalactic media. The physical basis for gas in thermal equilibrium to be subject to thermal instability (TI) was introduced in a classic paper by Field (1965), who developed the linear theory accounting for thermal conduction. Later Balbus (1986) generalized Field's theory by relaxing the equilibrium assumption and permitting the background flow to be time-dependent.

One specific instance where gas can be prone to TI is when both cooling and heating are due to radiative processes. This situation occurs, for example, in active galactic nuclei (AGNs) where, for a temperature of about ${{10}^{5}}$–${{10}^{6}}$ K, the cooling is dominated by free–free emission while the energy gain is dominated by X-ray absorption and Compton scattering of energetic photons produced by an accreting black hole (e.g., Krolik et al. 1981; Lepp et al. 1985). In this situation, the radiation field is almost in the free-streaming limit and consequently the gas will receive a non-zero net push in a direction away from the source of radiation. Our objective here is to demonstrate that this push can be substantial due to a large increase in the bound–free and bound–bound opacities of the cold gas during a nonlinear phase of the TI.

In this paper, we present results from our detailed study on the nonlinear dynamics of thermally unstable gas that is accelerated by the same photons that establish its thermal environment. In Section 2, we describe the theoretical underpinnings on which this paper is premised. In Section 3, we outline the equations solved and the explicit form of both the energy and momentum source terms. In Section 4, we present our methods and results from one-dimensional (1D) and two-dimensional (2D) simulations. We finish in Section 5 where we discuss our findings.

To facilitate a simple comparison between our simulations and Field's theory (i.e., the linear growth rates), we include thermal conduction, assume the gas to be initially stationary, homogenous, and non-magnetized, and we do not include gravity. We model the formation and evolution of optically thin clouds whose thermal equilibrium state is controlled by impinging radiation by adopting a realistic prescription for heating and cooling appropriate for gas in an AGN environment. Specifically, we use a net cooling function, ${\mathcal{L}}$, that includes the following radiative processes: (1) Compton and inverse-Compton scattering, (2) photoionization and recombination, (3) bremsstrahlung, and (4) line-emission. In an optically thin case, ${\mathcal{L}}$ depends on the gas temperature, T, and photoionization parameter, $\xi =4\pi {{{\mathcal{F}}}_{{\rm X}}}/n$, where ${{{\mathcal{F}}}_{{\rm X}}}$ is the X-ray flux and n is the gas number density.¹

The following theoretical picture forms the basis of our expectations and motivates our simulation setup. For a given ξ, gas in equilibrium at temperature T_eq will have heating balancing cooling (i.e., ${\mathcal{L}}({{T}_{{\rm eq}}},\xi )=0)$. Therefore, it will occupy one point on the radiative equilibrium curve (the solid line which is hereafter denoted the S-curve) shown in the top panel of Figure 1. Suppose that initially this point is at a stable location on the S-curve such as at position 1 in Figure 1, but some physical event transpires that results in a reduction of ${{{\mathcal{F}}}_{{\rm X}}}$ so that the gas finds itself out of equilibrium at location (2) with, for example, ${{\xi }_{2}}=190$. It will quickly cool to reside on the S-curve at location (3) (with ${{\xi }_{3}}={{\xi }_{2}}$). However, at this position, gas is thermally unstable to perturbations with constant pressure p, as it violates Field's criterion for stability, ${{[\partial {\mathcal{L}}/\partial T]}_{p}}\gt 0$. More generally, only the region below the dashed line is thermally stable according to Balbus' criterion ${{[\partial ({\mathcal{L}}/T)/\partial T]}_{p}}\gt 0$, with the consequence that gas occupying the shaded gray region in Figure 1 will evolve to points on the S-curve outside of the gray region. Continuing with our example, the slightest density perturbation present in the gas at location (3) will grow exponentially at the linear theory rate and at nearly constant pressure until the perturbation becomes nonlinear. Rapid nonlinear growth (still at nearly constant pressure) will then commence (e.g., Burkert & Lin 2000) and result in the formation of a cloud with a core density approximately determined by the intersection of the dotted line and the S-curve (position 4), which is where the gas in the core is stable.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 1 illustrates how various gas opacities depend on ξ. Solid black and red lines show the bound–free opacity, ${{\sigma }_{{\rm X}}}$, and the total line opacity, ${{M}_{{\rm max} }}$, respectively. The solid blue line represents the opacity of the most opaque line. One can see that for ${{\xi }_{3}}$, the bound–free and bound–bound opacity is negligible. (The vertical line in the gray region marks the location ${{\xi }_{3}}=190$.) However, once the thermally unstable gas forms a dense cloud, ${{\sigma }_{{\rm X}}}$ and ${{M}_{{\rm max} }}$ are significantly increased, as indicated by the left vertical line at ${{\xi }_{4}}\approx 73$. Therefore the cloud can be accelerated by the same source that heats it.

To confirm and quantify our expectations, we numerically solve the equations of hydrodynamics:

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

$$\begin{eqnarray}&&\frac{\partial (\rho {\boldsymbol{v}} )}{\partial t}+\nabla \cdot (\rho {\boldsymbol{vv}} +p\;{\sf I})={{{\boldsymbol{f}} }_{{\rm rad}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\frac{\partial E}{\partial t}+\nabla \cdot [(E+p){\boldsymbol{v}} ]=-\rho {\mathcal{L}}+{{\kappa }_{{\rm eq}}}{{\nabla }^{2}}\;T.\end{eqnarray} \tag{ 3 }$$

In the above, ρ is the gas density, ${\boldsymbol{v}}$ is the fluid velocity, p is the gas pressure, ${\sf I}$ is the unit tensor, ${{{\boldsymbol{f}} }_{{\rm rad}}}$ is the radiation force (see below), and $E=e+\rho {{v}^{2}}/2$ is the total energy density, with e the internal energy density. This system is closed with an ideal gas equation of state, $e=p/(\gamma -1)$, using an adiabatic index $\gamma =5/3$. We estimate the conduction coefficient ${{\kappa }_{{\rm eq}}}$ using the value for a fully ionized plasma (Spitzer 1962) evaluated at the equilibrium temperature and density of position 3 in Figure 1 (hereafter a suffix "eq" denotes evaluation of a quantity at $[{{\xi }_{3}},{{T}_{{\rm eq}}}]$).

We denote the functional dependence of the net cooling function as ${\mathcal{L}}={\Lambda }-{\Gamma }$, with

$$\begin{eqnarray}&&{\Lambda }=\frac{n}{\mu {{m}_{p}}}({{L}_{{\rm ff}}}+{{L}_{{\rm bb}}})\;[{\rm erg}\;{{{\rm g}}^{-1}}\;{{{\rm s}}^{-1}}]\ {\rm and}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\Gamma }=\frac{n}{\mu {{m}_{p}}}({{G}_{C}}+{{G}_{{\rm X}}})\;[{\rm erg}\;{{{\rm g}}^{-1}}\;{{{\rm s}}^{-1}}],\end{eqnarray} \tag{ 5 }$$

where μ is the mean molecular weight (set to 1.0 in this work) such that $n=\mu {{m}_{p}}\rho$ (with m_p the proton mass), ${{L}_{{\rm ff}}}$ and ${{L}_{{\rm bb}}}$ are the cooling rates due to free–free and bound–bound transitions, and G_C and G_X are the heating rates due to Compton and X-ray heating, respectively (all four rates are in units of ${\rm erg}\;{\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}$). For ${{L}_{{\rm ff}}}$ and G_C, we use the well-known analytic formulas based on atomic physics; see the Appendix. For ${{L}_{{\rm bb}}}$ and G_X, meanwhile, we use the analytical fits given by Blondin (1994), who found good agreement (to within 25%) between his approximate rates and those resulting from detailed photoionization calculations. Blondin's formulae, also provided in the Appendix, assume an optically thin gas of cosmic abundances illuminated by a ${{T}_{{\rm X}}}=10\;{\rm keV}/{{k}_{{\rm B}}}$ bremsstrahlung spectrum (with k_B Boltzmann's constant).

To evaluate ${{{\boldsymbol{f}} }_{{\rm rad}}}$, we consider the source of radiation to be located at $x=-\infty$, giving

$$\begin{eqnarray}&&{{{\boldsymbol{f}} }_{{\rm rad}}}=\frac{\rho {{\sigma }_{{\rm tot}}}{{{\mathcal{F}}}_{{\rm tot}}}}{c}\hat{x}\;[{\rm g}\;{\rm c}{{{\rm m}}^{-2}}\;{{{\rm s}}^{-2}}],\end{eqnarray} \tag{ 6 }$$

where ${{{\mathcal{F}}}_{{\rm tot}}}$ is the total continuum flux, c is the speed of light, and ${{\sigma }_{{\rm tot}}}$ is an effective total opacity, individual contributions to which can be self-consistently estimated as we describe below.

The X-ray flux is set by the photoionization parameter ξ. Therefore, we can specify the product ${{\sigma }_{{\rm tot}}}{{{\mathcal{F}}}_{{\rm tot}}}$ by assuming that there is nonzero continuum opacity only for X-rays and nonzero line opacity for UV photons. With this assumption we parametrize the UV flux ${{{\mathcal{F}}}_{{\rm UV}}}$ in terms of the X-ray flux using ${{f}_{{\rm UV}}}\equiv {{{\mathcal{F}}}_{{\rm UV}}}/{{{\mathcal{F}}}_{{\rm X}}}$ and use the following expression to estimate the opacity and flux product:

$$\begin{eqnarray}&&{{\sigma }_{{\rm tot}}}{{{\mathcal{F}}}_{{\rm tot}}}={{\sigma }_{e}}\left[ (1+{{f}_{{\rm UV}}})+{{\sigma }_{{\rm X}}}+{{f}_{{\rm UV}}}{{M}_{{\rm max} }} \right]{{{\mathcal{F}}}_{{\rm X}}}.\end{eqnarray} \tag{ 7 }$$

Here, ${{\sigma }_{e}}$ is the mass scattering coefficient of the free electrons. The term $(1+{{f}_{{\rm UV}}})$ in Equation (7) is due to electron scattering, while the term ${{\sigma }_{{\rm X}}}$ is an effective X-ray opacity in units of ${{\sigma }_{e}}$. We compute ${{\sigma }_{{\rm X}}}$ as $(4\pi {{G}_{{\rm X},h}}/\xi )/{{\sigma }_{e}}$, where ${{G}_{{\rm X},h}}$ is the heating part of G_X (see the Appendix). The last term in Equation (7), ${{M}_{{\rm max} }}$, is the maximum force multiplier characterizing the line contribution to the total opacity. We evaluate ${{M}_{{\rm max} }}$ following Stevens & Kallman (1990, SK 90 hereafter) with a modification, as described below. Note that in AGNs, ${{f}_{{\rm UV}}}\gt 1$.

The maximum force multiplier parametrizes the increase in the scattering coefficient due to spectral lines when all of the lines are optically thin. In this optically thin limit, the multiplier is just a sum of opacity contributions from all the lines and depends only on the gas composition, ionization, excitation and oscillator strengths (e.g., see Equations (10) and (11) in Castor et al. 1975; CAK hereafter). Following Owocki et al. (1988), who used a modified CAK method, one can evaluate this maximum multiplier as

$$\begin{eqnarray}&&{{M}_{{\rm max} }}={{k}_{{\rm CAK}}}(1-\alpha )\eta _{{\rm max} }^{\alpha },\end{eqnarray} \tag{ 8 }$$

where ${{k}_{{\rm CAK}}}$ is proportional to the total number of lines, α is the ratio of strong to weak lines, and finally ${{\eta }_{{\rm max} }}={{\kappa }_{L,{\rm max} }}/{{\sigma }_{e}}$ is a dimensionless measure of opacity of the most opaque line (with ${{\kappa }_{L,{\rm max} }}$ being the line opacity coefficient of the thickest line).

SK 90 carried out detailed photoionization calculations for a radiative environment appropriate for X-ray sources and parametrized their results in terms of the above expression for ${{M}_{{\rm max} }}$ by allowing ${{k}_{{\rm CAK}}}$ to be T-dependent and ${{\eta }_{{\rm max} }}$ to be ξ-dependent. Instead of the fit for ${{k}_{{\rm CAK}}}(T)$ from SK 90, we use Equation (17) of Proga (2007) due to T. Kallman (2006, private communication), as this may better represent the increase in the number of lines with decreasing temperature in AGNs. This expression and that for ${{\eta }_{{\rm max} }}(\xi )$, which is Equation (19) of SK 90, are provided in the Appendix. Both of these fits were generated assuming $\alpha =0.6$. In the bottom panel of Figure 1, we plot ${{\sigma }_{{\rm X}}}$ along with ${{M}_{{\rm max} }}$. Notice that M_max can be roughly a few thousand for gas ionized by a weak radiation field, whereas it decreases asymptotically to zero for highly ionized gas as the radiation field becomes stronger.

We solve Equations (1)–(3) in 1D and 2D using the CTU integrator, Roe flux, and explicit conduction module of the MHD code Athena (Stone et al. 2008). We modify the original version of the code by adding the momentum and heating and cooling source terms in the same way that Athena's built-in static gravitational potential source term is implemented to achieve 2nd order accuracy in both space and time. We use a less accurate method for integrating the conduction term in time in our 2D simulations than in our 1D simulations. Specifically, we use Athena's super time-stepping scheme (STS; see O'Sullivan & Downes 2006), although we note that a second order accurate in time STS algorithm does exist (Meyer et al. 2012) and we are testing it for future use.

4.1. Initial and Boundary Conditions

Given the atomic physics behind our S-curve and the opacities in Figure 1, the following free parameters govern our problem: the wavenumber and density amplitude of the TI perturbation, k and $\delta \rho$, as well as the ratio of the initial thermal and sound crossing times ${{t}_{{\rm th}}}/{{t}_{{\rm sc}}}$, which together determine the number of clouds and their formation time; the initial photoionization parameter ${{\xi }_{3}}$, which controls the intensity of the radiation field and the equilibrium temperature T_eq; the equilibrium pressure, namely the product ${{n}_{{\rm eq}}}{{T}_{{\rm eq}}}$, which sets the physical units of the cloud and its environment; and finally f_UV, which parametrizes the shape of the spectral energy distribution in a simple way.

Our initial conditions are full wavelength profiles of the TI condensation mode, found from Equations (11) to (14) in Field (1965), applied to density, velocity, and pressure. We adopt ${{n}_{{\rm eq}}}{{T}_{{\rm eq}}}={{10}^{13}}\;{\rm K}\;{\rm c}{{{\rm m}}^{-3}}$ in accordance with AGN observations and their modeling (e.g., Davidson & Netzer 1979; Krolik et al. 1981). The length scale of the perturbation is fixed by the adiabatic sound speed at position 3 in Figure 1, c_eq, and a choice for the initial sound crossing time, t_sc. We chose t_sc equal to the initial thermal time, ${{t}_{{\rm th}}}=({{e}_{{\rm eq}}}/{{\rho }_{{\rm eq}}})/{{{\Lambda }}_{{\rm eq}}}$, which results in near maximum linear growth rates for the condensation mode. These rates are obtained by solving the dispersion relation in Field (1965; Equation (15)). For both 1D and 2D simulations, we use periodic boundary conditions and set the domain size in the x-direction, ${\Delta }x$, equal to the perturbation wavelength ${{\lambda }_{x}}={{c}_{{\rm eq}}}{{t}_{{\rm sc}}}$. The amplitude of the density perturbation is $\delta \rho =5.0\times {{10}^{-5}}{{\rho }_{{\rm eq}}}$.

With this setup, the unstable region of Figure 1 is parametrized entirely by ${{\xi }_{3}}$, and varying this parameter leads to the formation of clouds with substantially different properties. A realistic value for AGNs is likely ${{\xi }_{3}}=500$, as this results in a pressure photoionization parameter ${\Xi }\equiv ({{{\mathcal{F}}}_{{\rm X}}}/c)/({{n}_{{\rm eq}}}{{k}_{{\rm B}}}{{T}_{{\rm eq}}})\approx 9.0$, and the AGN environment is expected to be hospitable to clouds for ${\Xi }\lesssim 10$ (e.g., Krolik et al. 1981). For ${{\xi }_{3}}=500$, the linear growth rate of the TI is comparatively small, taking more than 400 days for the density of the cold gas to double, and estimates based on Figure 1 indicate that a 1D (or 2D planar) cloud will form at ${{\xi }_{c}}\approx 20$ with a width ${{l}_{c}}\sim 0.5\;{\rm AU}$, a temperature ${{T}_{c}}\sim 4\times {{10}^{4}}\;{\rm K}$, a number density ${{n}_{c}}\sim 3\times {{10}^{8}}\;{\rm c}{{{\rm m}}^{-3}}$, and a density contrast of $\chi \sim 30$. The radiation force due to lines can be very powerful, with ${{M}_{{\rm max} }}$ at least 10³.

This more realistic cloud is, however, very optically thick for many strong UV lines. (Moreover, as we discuss in Section 5, it is challenging to accurately resolve in multi-dimensions.) Our present simulations are designed to explore the optically thin regime, as this regime allows us to focus on the purely hydrodynamical effects of cloud formation and acceleration without the further complications involved when solving the equations of radiation hydrodynamics (RHD; cf. Proga et al. 2014). This restricts ${{\xi }_{3}}$ to a very narrow range of values corresponding to larger, less realistic values of Ξ, as there is obviously an upper limit on the density of clouds whose acceleration we can accurately model without RHD.

To estimate this upper limit on the density, we assume the cloud forms with both a constant density ${{\rho }_{c}}$ and width l_c (which will be seen to be very nearly the case here), so that we can write ${{\eta }_{{\rm max} }}={{\tau }_{L,{\rm max} }}/{{\tau }_{s}}$, where ${{\tau }_{L,{\rm max} }}={{\kappa }_{L,{\rm max} }}{{\rho }_{c}}{{l}_{c}}$ is the optical depth of the thickest line and ${{\tau }_{s}}={{\sigma }_{e}}{{\rho }_{c}}{{l}_{c}}$ is the electron scattering optical depth of the cloud. Demanding that the cloud to be optically thin to all bound–bound transitions (i.e., ${{\tau }_{L,{\rm max} }}\lt 1$) requires ${{\eta }_{{\rm max} }}{{\tau }_{s}}\lt 1$, or

$$\begin{eqnarray}&&{{\rho }_{c}}\lt {{\left( {{\sigma }_{e}}{{l}_{c}}{{\eta }_{{\rm max} }} \right)}^{-1}}.\end{eqnarray} \tag{ 9 }$$

We chose the value ${{\xi }_{3}}=190$ since it produces a cloud with the highest density contrast that satisfies this inequality at all times in 1D. (In Section 4.4, we verify our optically thin assumption in both 1D and 2D using a more accurate estimate.) For ${{\xi }_{3}}=190$, ${{M}_{{\rm max} }}$ does not exceed 40. To explore the effects of the stronger line force, we set ${{f}_{{\rm UV}}}=10$, which is guided by observational results from Zheng et al. (1997) and Laor et al. (1997).

4.2. Simulations

We have performed over 30 simulations exploring a variety of parameters and numerical setups. Here we report in some detail on a set of four simulations that differ only by the applied radiation force: RF (electron scattering only), RFX (electron scattering plus X-ray absorption), RFLD (electron scattering plus line-driving), and RFLDX (electron scattering, line-driving, and X-ray absorption). We will especially focus on run RFLDX, as this case is most representative of the physical conditions in AGNs. Compared to a more realistic AGN cloud described above, for our adopted value of ${{\xi }_{3}}=190$ the cloud growth is much faster, taking only 2.4 days to double in density, while ${{l}_{c}}\approx 6.5\times {{10}^{10}}\;{\rm cm}$ ($0.004\;{\rm AU}$), ${{T}_{c}}\approx 7.0\;\times {{10}^{4}}\;{\rm K}$, ${{n}_{c}}\approx 1.4\times {{10}^{8}}\;{\rm c}{{{\rm m}}^{-3}}$, $\chi \approx 8.0$, and ${\Xi }\approx 19$. The physical units corresponding to our numerical results are listed in Table 1.

Table 1.  Physical Units

Corresponding to $\xi =190$ and ${{t}_{{\rm th}}}/{{t}_{{\rm sc}}}=1$
Quantity	Value (cgs units)
${{\rho }_{{\rm eq}}}$	8.65 × 10−17 g cm−3
neq	5.17 × 107 cm−3
Teq	1.93 × 105 K
ceq	5.16 × 106 cm s-1
tsc	6.03 × 104 s
λx	3.11 × 1011 cm
${{{\mathcal{F}}}_{{\rm X}}}$	7.82 × 108 erg s−1 cm−2
κeq	1.58 × 107 erg s−1 K−1 cm−1
${{{\Lambda }}_{{\rm eq}}}$	3.97 × 108 erg g−1 s−1

Download table as:  ASCIITypeset image

4.3. Results of 1D Simulations

Despite the different radiation forces, the clouds in all four runs are formed at the same time and with the same density and temperature contrasts, and the gas therefore traces the same "tracks" on the $T-\xi$ plot in Figure 1. The over-plotted red and blue dots in Figure 1 show the tracks for RFLDX at $t=120\;{{t}_{{\rm sc}}}$, which represents the time where the flow reached a thermal steady state (see below for more details). Note that the red tracks do not reach the radiative equilibrium curve (i.e., ${\mathcal{L}}=0),$ but rather an equilibrium curve given by $\rho {\mathcal{L}}={{\kappa }_{{\rm eq}}}{{\nabla }^{2}}\;T$.

Also note that there are tracks occupying an unstable (according to Balbus' criterion) region in Figure 1, namely, the tracks that are above the dashed line. Given that the gas in the cloud core occupies location (4) and is in pressure equilibrium with the medium, it must be the case that some portion of the gas occupies this unstable region in order for the density and temperature to be continuous everywhere. These "unstable" tracks correspond to the gas in the conductive interfaces of the cloud. In Figure 2, we plot profiles of the solution overplotted in Figure 1, and the width of the interfaces can be judged from the density panel. The local Field length in the interfaces is close to the initial equilibrium value, ${{\lambda }_{F}}=2\pi \sqrt{{{\kappa }_{{\rm eq}}}{{T}_{{\rm eq}}}/({{\rho }_{{\rm eq}}}{{{\Lambda }}_{{\rm eq}}})}\;{{\lambda }_{x}}\approx 0.19\;{{\lambda }_{x}}$, while the interface width is two or three times smaller than this. Gas is permitted to be thermally unstable at regions smaller than the Field length when the stabilizing influence of conduction (measured by the heat flux that is shown in the fifth panel of Figure 2) is large enough.

Zoom In Zoom Out Reset image size

We find that in all of these runs the gas arrives at a state of thermal equilibrium by $t=90\;{{t}_{{\rm sc}}}$. The second from bottom panel in Figure 2 shows how this equilibrium state is possible. Both the net cooling function (red curve) and the conduction term (black curve) are positive at the interfaces of the cloud and negative in the hot medium. These terms are of opposite signs in Equation (3) and therefore balance each other. The compression term (cyan) is negligible at this time, but it is a dominant term when the cloud is forming.

The radiation force prevents the gas from reaching a mechanical equilibrium state. Rather, in each case the cloud core undergoes dynamical changes (i.e., the pressure and velocity profiles adjust) to permit nearly uniform acceleration. To show this, we plot the net flow acceleration in the bottom panel of Figure 2 (cyan line), which is the sum of the other curves displayed. Line driving operates almost uniformly throughout the cloud core. As can be seen by either the acceleration or the pressure panel, the response of the gas pressure is to exert a nearly constant drag force on the cloud to compensate for the driving force, while the medium is pushed along since it has nowhere else it can go in 1D. The adjustment of the forces is obtained shortly after the cloud is formed, with the acceleration profiles resembling those shown in the bottom panel at around $t=55\;{{t}_{{\rm sc}}}$ for runs RFLD and RFLDX. Some profiles (especially velocity) continue to undergo changes until $t=120\;{{t}_{{\rm sc}}}$; the shapes shown in Figure 2 are maintained as the cloud continues to accelerate to high Mach numbers.

The results from our 1D simulations confirm our basic picture for cloud formation and acceleration. As a next step, we perform 2D simulations where destructive processes may change our results. Our primary concern is the Rayleigh–Taylor (RT) instability, and subsequently the Kelvin–Helmholtz (KH) instability. Figure 2 shows that the left interface is RT stable, as the density increases in the direction of acceleration. However, the right interface is likely unstable due to the adverse density arrangement (e.g., Krolik 1977, 1979; Mathews & Blumenthal 1977; Jacquet & Krumholz 2011; Jiang et al. 2013).

4.4. Results of 2D Simulations

The simplest extension of our 1D simulations to 2D is to form a cloud with a planar slab configuration. The initial conditions and overall setup is as before. However, to fully explore 2D effects, we now break the uniformity in the y-direction by introducing a perturbation with a wavelength the size of the domain in the y-direction (i.e., ${{\lambda }_{y}}={{\lambda }_{x}}/2$) and $\delta \rho =5\times {{10}^{-7}}{{\rho }_{{\rm eq}}}$.

In Figure 3, we present a comparison of our four 1D runs and their 2D counterparts, and we also verify our optically thin assumption. The maximum density of the cloud versus time is plotted in the top panel. It is clear that there is no significant difference in any of the runs during the cloud formation process, which ends at $t\approx 50\;{{t}_{{\rm sc}}}$. As mentioned in Section 4.1, for runs RFLD and RFLDX to be optically thin to UV photons, we require ${{\tau }_{L,{\rm max} }}\lt 1$. We estimate ${{\tau }_{L,{\rm max} }}$ as²

$$\begin{eqnarray}&&{{\tau }_{L,{\rm max} }}={{\sigma }_{e}}\int _{0}^{{{\lambda }_{x}}}{{\eta }_{{\rm max} ,90}}(x)\rho (x)dx,\end{eqnarray} \tag{ 10 }$$

where ${{\eta }_{{\rm max} ,90}}$ denotes only those values of ${{\eta }_{{\rm max} }}$ in the range $[0.9\;{{\eta }_{{\rm max} }},{{\eta }_{{\rm max} }}]$. We use this range to be able to identify gas at a constant opacity in our numerical representation of ${{\eta }_{{\rm max} }}$.

Zoom In Zoom Out Reset image size

The black curves in the top panel of Figure 3 show this estimate for ${{\tau }_{L,{\rm max} }}$ in both 1D (dashed line) and 2D (solid line). In the latter calculation, we consider a ray through the center of the cloud at $y=0.25{{\lambda }_{x}}$. Overall, the cloud can indeed be considered optically thin at all times, except possibly at its center during a very short period of the acceleration phase in 2D, which as will be made clear below, coincides with when the cloud is significantly lengthened by the onset of disruptive processes.

The bottom panel of Figure 3 shows the average velocity of the cloud versus time, where we define the cloud as being the gas to the left of the gray region in Figure 1 (i.e., $\xi \lt 121.2$, which corresponds to $\rho \gt 1.57{{\rho }_{{\rm eq}}}$). It indicates that significant acceleration only takes place after the cloud formation process has ended, as line opacity is only activated once the cloud forms. While 1D runs of RFLD and RFLDX uniformly accelerate to supersonic speeds, the acceleration is suddenly halted in 2D around $t\approx 66\;{{t}_{{\rm sc}}}$.

In Figure 4, we plot snapshots of run RLFDX to illustrate the very different fate of a rapidly accelerated cloud in 2D. The first frame shows that the slab formed at $t=45\;{{t}_{{\rm sc}}}$. The initial perturbation in the y-direction grows into a slight over-density region in the center of the slab, which then undergoes greater acceleration than its surroundings and causes a small bulge to appear around $t\approx 55\;{{t}_{{\rm sc}}}$ (not shown). Viewing this bulge as a perturbation along the surface of the slab, the basic criteria for the RT instability is satisfied: heavy fluid is pushing against light fluid. The same conclusion applies to runs RFX and RFLD, although for run RFX it will take much longer (several hundred t_sc) for the bulge to grow due to the weak acceleration. Run RF, however, which has a constant acceleration due to Thomson opacity, evolves identically in 1D and 2D for all time; any density perturbation in the y-direction receives the same push as any other point in the flow.

Zoom In Zoom Out Reset image size

The remaining frames in Figure 4 reveal how the breakup of the cloud ensues as the RT instability develops and soon becomes accompanied by the KH instability. First the bulge evolves to become mushroom-shaped, forming a structure resembling the classic RT "spike," which features prominent KH "rolls" at $t=66\;{{t}_{{\rm sc}}}$ made possible by the increased relative velocity between the cloud and medium. The halting of the acceleration happens around this time. The connecting plume then disperses (i.e., its gas is heated) as the spike separates further from the slab. The slab would likely disperse also, due to the mass lost to the spike, but instead it is somewhat thickened by the approach of the spike from the backside. This collision perspective is shown in the final panel at $t=75\;{{t}_{{\rm sc}}}$.

The cloud is eventually either completely dispersed back into the medium, or coexists with it in a disordered manner in what could be called a clumpy flow once the vertical symmetry is lost. The simulations presented here do not let us make any definitive statements because we noticed that our solutions become resolution dependent at about time $73\;{{t}_{{\rm sc}}}$ for run RFLDX. This loss of convergence is shown in Figure 5, where for four different resolutions, we plot the mass fraction of the three components of the gas: (i) the cloud (blue), again defined as a gas with $\xi \lt 121.2$ or $\rho \gt 1.57\;{{\rho }_{{\rm eq}}}$; (ii) the interface (black), defined as the portion of the gas in the gray region that is unstable, i.e., the tracks above the dashed line in Figure 1 with $121\leqslant \xi \leqslant 278$ or $0.68\;{{\rho }_{{\rm eq}}}\leqslant \rho \leqslant 1.57\;{{\rho }_{{\rm eq}}}$; and (iii) the medium (red), defined as the (stable) gas with $\xi \gt 278$ or $\rho \lt 0.68\;{{\rho }_{{\rm eq}}}$. The mass fraction is defined as the mass of each component divided by the total mass and is given by $m={{({{N}_{x}}{{N}_{y}})}^{-1}}{{\sum }_{ij}}{{\rho }_{ij}}/{{\rho }_{{\rm eq}}}$, where N_x and N_y are the number of grid zones in the x- and y-directions and the sum ranges over all zones (i, j) that satisfy one of the criteria (i)–(iii). Once the mass fractions for resolutions N_x = 1024 and N_x = 2048 differ, we cannot claim to accurately follow the cloud's evolution.

Zoom In Zoom Out Reset image size

The mass fractions provide a complementary description of the cloud evolution depicted in Figure 4. The cloud appears fully formed by $t\approx 50\;{{t}_{{\rm sc}}}$, which coincides with the peak mass fraction of the medium, but mass keeps piling on until $t\approx 63\;{{t}_{{\rm sc}}}$. Indeed, we observe that the velocity arrows at $t=50\;{{t}_{{\rm sc}}}$ in Figure 4 point toward the cloud, indicating continued growth at the expense of the medium. The overall fraction of gas occupying interface regions is a minimum at $t\approx 60\;{{t}_{{\rm sc}}}$, and this is despite the overall increase in the size of interface region (due to the bulge) because the interfaces are narrower. The cloud mass fraction reaches a maximum at $t\approx 63\;{{t}_{{\rm sc}}}$ when the RT spike has become fully nonlinear, and thereafter monotonically decreases, with the mass being taken up entirely in the interfaces, until the RT spike becomes detached from the slab. During this time, the medium continues losing mass to the interfaces. The loss of convergence is likely due to the appearance of small scale structures as the cloud is disrupted.

4.5. 2D Convergence Study

The resolution dependent late-time dynamics encountered above warrants demonstrating that our results are converged at earlier times, especially since an acceleration scheme such as STS was found to be necessary in 2D due to the strict time constraint imposed by the CFL condition (${\Delta }t\propto {\Delta }{{x}^{2}}$ due to conduction) combined with the long duration of the runs. To show that our 2D results using STS are converged, in the top panel of Figure 6 we plot L₁ errors for run RFLDXa at $t=45\;{{t}_{{\rm sc}}}$. This run is a modified version of run RFLDX where we suppressed the contribution of the force due to electron scattering by a factor of 11 (by setting the term $(1+{{f}_{{\rm UV}}})$ in Equation (7) equal to 1) to minimize the effects of advection errors not associated with the dominant term ${{f}_{{\rm UV}}}{{M}_{{\rm max} }}$.

Zoom In Zoom Out Reset image size

For a uniform mesh the L₁ error norm is the quantity ${{({{N}_{x}}{{N}_{y}})}^{-1}}{{\sum }_{ij}}{{R}_{ij}}$, where ${{R}_{ij}}=|{{q}_{ij}}-{{q}_{{\rm ref}}}|$ is the local absolute error in some quantity q at the location of zone (i, j) with respect to some reference solution q_ref evaluated at the same location, and the sum ranges over all zones in the domain. As a reference solution we used the highest resolution simulation we could reasonably afford, namely $[{{N}_{x}},{{N}_{y}}]=[2048,1024]$. The dashed line in Figure 2 marks the slope corresponding to a convergence rate of 2, showing that the solutions for density and velocity converge to the reference solution in line with Athena's second order accurate spatial discretization. However, pressure shows a reduced convergence rate, which may be due to the fact that the STS scheme is only first order accurate in time. Interestingly, we observed that reducing the Courant number to 0.125 will result in pressure also showing order convergence, but our results presented here used a Courant number 0.5.

In Figure 7 we plot relative errors for pressure, defined as ${{R}_{ij}}/|{{p}_{{\rm ref}}}|$, for a horizontal slice through the domain at j = 0. This plot also demonstrates self-convergence and, not surprisingly, it shows that the largest errors are found at the conductive interfaces of the cloud.

Zoom In Zoom Out Reset image size

As the next planned phase of our ongoing investigation of cloud acceleration initiated in Proga et al. (2014), this paper considered the cloud formation and acceleration processes simultaneously for the first time. In particular, we have extended the basic theory of the nonlinear outcome of TI by self-consistently solving for the dynamics of optically thin gas in the presence of a strong radiation field. Our resulting simulations have made it possible to study in detail the evolution of the isobaric condensation mode in thermally unstable gas from an initial perturbation to a dense, high velocity cloud. The cloud forms in a radiation pressure dominated environment, but the radiation force has practically no effect on the cloud formation process in either 1D or 2D. The reason is simply because there is no momentum transfer from the radiation to the gas unless there is also sufficient opacity, and the primary sources of opacity are not activated until the cloud is formed, a point that has been previously noted by, for example, Krolik (1988).

The initial motivation for this work was simply to demonstrate that the nonlinear phase of TI leads to a natural mechanism to produce fast clouds via acceleration due to lines. In so doing we were led to the inescapable conclusion that accelerated clouds undergo rapid deformation and are ultimately destroyed, thereby confirming long-standing assertions about the inevitability of cloud destruction (e.g., Mathews 1986; Krolik 1999). However, we find that optically thin clouds can survive long enough to be accelerated to relatively high velocities and travel a significant distance of many cloud sizes, in contrast to investigations exploring pre-existing optically thick clouds (Proga et al 2014 and references therein). Consequently, the best hope for cloud-based models of AGNs is to identify robust mechanisms for continually producing clouds (e.g., in "outflows from inflows" as illustrated in simulations presented by Kurosawa & Proga 2009; Mościbrodzka & Proga 2013). It is therefore important to thoroughly study how clouds form via TI and how they evolve before disruptive processes take hold.

We found that a rich set of dynamics unfolds during the cloud formation process. For example, for run RFLDX (radiation force due to electron scattering, X-rays, and lines) in 1D, there are three nonlinear phases of the TI: (i) initial growth and saturation $(t\approx 40-46\;{{t}_{{\rm sc}}})$; (ii) evolution toward a uniformly accelerating solution $(t\approx 46-55\;{{t}_{{\rm sc}}})$; and (iii) evolution toward a thermal equilibrium state $(t\approx 46-90\;{{t}_{{\rm sc}}})$. Figure 4 shows that phases (i) and (ii) have both completed before the RT instability can develop, so these phases also take place in 2D, while phase (iii) has a different outcome in 2D and the uniform acceleration cannot be maintained. These phases will be described in more detail in T. Waters & D. Proga (2015, in preparation).

Our follow up paper will also further investigate the growth of the RT and KH instabilities. This is an important matter since the appearance of these instabilities governs the end phase of TI and therefore dictates the ultimate fate of the cloud. These instabilities develop after the TI saturates. Therefore, it should be possible to confirm the theoretical linear growth rates by conducting a careful numerical perturbation analysis (e.g., by introducing sinusoidal perturbations along the interfaces of a 2D planar slab initialized using a 1D solution). That said, we did not find it possible to make a meaningful comparison of the growth rate of the bulge and the classical RT rate in the present setup. Recalling Figure 3, the bulge is already borderline nonlinear by $t=60\;{{t}_{{\rm sc}}}$, implying that the slab is still evolving thermally (i.e., in phase (iii) of nonlinear TI evolution) during the linear RT growth regime. This dynamical complication combined with the simplifying assumptions inherent in the linear theory for RT (such as constant acceleration everywhere) warrant using a more controlled approach for isolating the development of the individual instabilities.

The numerical setup presented here can be used for several different exploratory studies of cloud formation and acceleration. We chose the initial perturbation amplitude $\delta \rho$ in the y-direction to be 10⁻² that in the x-direction in order to arrive at a slab configuration. With an equal ratio (and equal growth rates), a round cloud will be formed in 2D instead. Multiple clouds can be formed using higher wavenumber perturbations. Cloud fragmentation can be studied by decreasing the ratio ${{t}_{{\rm th}}}/{{t}_{{\rm sc}}}$. Classical evaporation (Cowie & McKee 1977; Balbus 1985) can be explored by arranging for the Field length to exceed the size of the cloud. This large range of initial configurations would be difficult or impossible to construct otherwise, which illustrates an obvious advantage of studying cloud acceleration via the formation process, namely that the internal gas dynamics are self-consistently treated. Indeed, we find that pressure equilibrium with the surrounding medium is naturally maintained, cloud interfaces form with a width determined by the conductivity, the radiation and drag forces reach a balance so that hydrostatic equilibrium is established in the reference frame of the cloud, and the cloud reaches a thermal equilibrium state in which heating by conduction is balanced by line cooling, in agreement with Begelman & McKee (1990). Moreover, the equilibrium location on the S-curve largely determines the cloud density, temperature, and opacity before it is accelerated, while plotting the evolutionary tracks on the $T-\xi$ plane has shown itself to be a useful tool both for understanding the time evolution of the cloud and for characterizing the components of the gas.³

Taking into consideration the numerical requirements involved in this study, multi-dimensional simulations will likely be constrained to only explore values of ${{\xi }_{3}}\lesssim 300$ for the S-curve used here, thereby limiting the overall density contrast of clouds to $\chi \approx 18$. This limitation arises due to the need to resolve interfaces, as simulations of clouds formed via TI will not be converged unless the conductive interfaces are resolved (Koyama & Inutsuka 2004). Previous numerical studies using pre-existing clouds without thermal conduction and with unresolved interfaces were able to explore much higher density contrasts such as $\chi =50$ (McCourt et al. 2014) and even $\chi \sim {{10}^{4}}$ (Krause et al. 2012). We found that a realistic $\kappa \propto {{T}^{5/2}}$ would require upwards of ${{N}_{x}}=4,096$ zones (i.e., at least three levels of refinement if adaptive mesh refinement is employed) even for our modest density contrast of $\chi \approx 8$, as would values of $\chi \gtrsim 30$ with continued use of a constant conductivity. This rapid steepening of the interfaces with either increased χ or realistic $\kappa (T)$ results in ever smaller transition regions between the cloud and the medium but does not imply that the role of thermal conductivity becomes less important. As a consequence, the conduction time step would be so small at these resolutions that even STS schemes would become impractically slow, necessitating the use of implicit techniques. Such extensions to this work may be needed to assess, for example, if the the timescale for cloud destruction is sensitive to the slope of the interfaces, i.e., if it decreases with steeper density gradients, as was found to be the case when a cloud is disrupted by the passage of a shock (Nakamura et al. 2006).

Future work is also needed to address important aspects of cloud formation and acceleration neglected here. Magnetic fields can play a dominant role in both, as they can prevent the gas from condensating in the first place (e.g., Mathews & Doane 1990), they may help accelerate clouds through confinement (e.g., Rees 1987; Arav & Li 1994), and they can significantly effect cloud dynamics due to the affects of anisotropic conduction (Choi & Stone 2012). Other environmental factors such as the presence of Coriolis forces, gravity, winds, shocks, etc. could be incorporated to build a more global model of cloud dynamics for comparison with observations. However, it is crucial to model 3D effects, as both compression and expansion of clouds can be enhanced by a large factor, and this can affect the timescale for cloud destruction. A direct extension of this study designed to form initially spherical clouds in 3D would violate our optically thin assumption and therefore requires a full RHD treatment using methods that can accurately capture shadows (e.g., Davis et al. 2012; Jiang et al. 2012). Such methods, which have already been employed by Proga et al. (2014), would also allow investigation of the so-called line-deshadowing instability (Owocki & Rybicki 1984) that could lead to changes in the structure and time-dependence of the clouds.

This work was supported by NASA under ATP grant NNX11AI96G. Our simulations were performed on the Eureka cluster at UNLV's National Supercomputing Center for Energy and the Environment. We thank J. Stone for developing the Athena code and making it public, as well as for discussions about this work. We are grateful to the referee for constructive comments that helped us improve this work. T.W. thanks Chad Meyer for discussions regarding the use and accuracy of STS schemes.

The net cooling function ${\mathcal{L}}$ that defines our S-curve is comprised of four heating and cooling rates. The analytic expressions for ${{L}_{{\rm ff}}}$ and G_C are

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{L}_{{\rm ff}}} & = & \frac{{{2}^{5}}\pi {{e}^{6}}}{3\;h{{m}_{e}}{{c}^{3}}}\sqrt{\frac{2\pi {{k}_{{\rm B}}}T}{3{{m}_{e}}}}{{Z}^{2}}{{{\bar{g}}}_{B}} \\ {} & = & 3.3\times {{10}^{-27}}\sqrt{T}\;[{\rm erg}\ {\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}],\;{\rm and} \\ \end{array}\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{G}_{C}} & = & \frac{{{k}_{{\rm B}}}{{\sigma }_{e}}}{4\pi {{m}_{e}}{{c}^{2}}}\xi \;{{T}_{{\rm X}}}\left( 1-4\frac{T}{{{T}_{{\rm X}}}} \right) \\ {} & = & 8.9\times {{10}^{-36}}\xi {{T}_{{\rm X}}}\left( 1-4\frac{T}{{{T}_{{\rm X}}}} \right)[{\rm erg}\;\ {\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}], \\ \end{array}\end{eqnarray} \tag{ A2 }$$

where m_e and e are the mass and charge of an electron, h is Planck's constant, Z is the ion atomic number, ${{\bar{g}}_{B}}$ is an averaged Gaunt factor, and ${{T}_{{\rm X}}}=10\;{\rm keV}/{{k}_{{\rm B}}}$. The analytical fits from Blondin (1994) in our notation read

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{L}_{{\rm bb}}} & = & \delta \left[ 1.7\times {{10}^{-18}}{\rm exp} \left( -\frac{1.3\times {{10}^{5}}}{T} \right){{\xi }^{-1}}{{T}^{-1/2}}+{{10}^{-24}} \right] \\ {} & {} & \times [{\rm erg}\;{\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}]\ {\mkern 1mu} {\rm and} \\ \end{array}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{G}_{{\rm X}}} & = & 1.5\times {{10}^{-21}}{{\xi }^{1/4}}{{T}^{-1/2}}\left( 1-\frac{T}{{{T}_{{\rm X}}}} \right) \\ {} & {} & \times [{\rm erg}\;{\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}]. \\ \end{array}\end{eqnarray} \tag{ A4 }$$

Here, δ is a parameter introduced by Blondin (1994); setting $\delta \lt 1$ mimics reducing the strength of line cooling when relaxing his assumption of optically thin gas. We keep $\delta =1$ since we assume optically thin clouds. In Section 3, we refer to the heating part of G_X, which is ${{G}_{{\rm X},h}}=1.5\times {{10}^{-21}}{{\xi }^{1/4}}{{T}^{-1/2}}\;[{\rm erg}\;\ {\rm c}{{{\rm m}}^{3}}\;{{{\rm s}}^{-1}}]$. Finally, it is important to note that Blondin's photoionization calculations were revisited and independently verified using XSTAR by Dorodnitsyn et al. (2008) for an incident AGN power-law spectrum. Their analytical fits differ from the above only by a minor modification to Equation (A3), which Dorodnitsyn (2008) report had no significant dynamical effects on their simulation results.

For line-driving, we use Equation (17) from Proga (2007)

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {\rm log} {{k}_{{\rm CAK}}} \\ \quad =\,\left\{ \begin{array}{ccccccccccccccc} -0.383 & {\rm for}\;{\rm log} T\leqslant 4 \\ -0.630{\rm log} T+2.138 & {\rm for}\;4\lt {\rm log} T\leqslant 4.75 \\ -3.870{\rm log} T+17.528 & {\rm for}\;{\rm log} T\gt 4.75 \\ \end{array} \right. \\ \end{array}\end{eqnarray} \tag{ A5 }$$

and Equation (19) from Stevens & Kallman (1990):

$$\begin{eqnarray}{\rm log} {{\eta }_{{\rm max} }}=\left\{ \begin{array}{ccccccccccccccc} 6.9\;{\rm exp} (0.16\;{{\xi }^{0.4}}) & {\rm for}\;\;{\rm log} \xi \leqslant 0.5 \\ {} & {} \\ 9.1{\rm exp} (-7.96\times {{10}^{-3}}\xi ) & {\rm for}\;{\rm log} \xi \gt 0.5. \\ \end{array} \right.\end{eqnarray} \tag{ A6 }$$