Large-scale Dynamics of Winds Originating from Black Hole Accretion Flows. I. Hydrodynamics

Can Cui; Feng Yuan; Bo Li

doi:10.3847/1538-4357/ab6e6e

1. Introduction

Observational evidence of winds from accretion disks has been accumulating for both cold and hot accretion flows over the past two decades. Cold accretion disks produce high-velocity winds, which have been widely observed in luminous active galactic nuclei (AGNs; e.g., Crenshaw et al. 2003; Tombesi et al. 2010, 2014; Liu et al. 2013; King & Pounds 2015) and black hole X-ray binaries (e.g., Neilsen & Homan 2012; Díaz Trigo & Boirin 2016; Homan et al. 2016). Detections of winds from hot accretion flows are more challenging because they are usually fully ionized. Nevertheless, wind observations from low-luminosity sources have gradually built up in recent years through various approaches. These detections include low-luminosity AGNs and radio galaxies (Tombesi et al. 2010, 2014; Crenshaw & Kraemer 2012; Cheung et al. 2016), the supermassive black hole in the Galactic center (Sgr A*; Wang et al. 2013; Ma et al. 2019), and black hole X-ray binaries in a hard state (Homan et al. 2016; Munoz-Darias et al. 2019).

Theoretically, the launching of winds has been extensively studied. For thin disks, three driving mechanisms have been proposed (Proga 2007), namely the thermal mechanism (e.g., Begelman et al. 1983) either in the context of circumstellar disks (e.g., Font et al. 2004) or in black hole accretion disks (e.g., Luketic et al. 2010; Waters & Proga 2012; Clarke & Alexander 2016), radiation (line force; e.g., Murray et al. 1995; Proga et al. 2000; Proga & Kallman 2004; Nomura & Ohsuga 2017), and the magnetic field (e.g., Blandford & Payne 1982; Contopoulos & Lovelace 1994; Romanova et al. 1997; Cao 2014; Bai et al. 2016). It is likely that all three mechanisms play a role, although most existing works focus only on one mechanism because of the technical difficulties. Two examples of exceptions are Proga (2003) and Waters & Proga (2018), in which both the thermal and magnetic mechanisms are taken into account.

Despite the relative rarity of observational evidence for winds launched from hot accretion flows, the theoretical understanding is more advanced than in the case of thin disks. This is partly because radiation is dynamically unimportant in hot accretion flows and because it is technically easier to simulate geometrically thick flows. It has long been suspected that strong winds exist in hot accretion flows (Narayan & Yi 1994; Blandford & Begelman 1999; Stone et al. 1999; Stone & Pringle 2001). This speculation was later confirmed by detailed numerical simulations (Narayan et al. 2012; Yuan et al. 2012a, 2012b; Li et al. 2013) and analytical studies (Bu et al. 2016a, 2016b; Bu & Mosallanezhad 2018). Yuan et al. (2015; hereafter Y15) study the properties of wind originating from hot accretion flows based on three-dimensional (3D) general relativistic (GR) magnetohydrodynamical (MHD) simulations, such as velocity and mass flux. They analyze data with a virtual particle trajectory approach that can effectively distinguish real wind from turbulent motions. This approach ensures that the wind properties obtained are valid.

Almost all the works mentioned above focus on the accretion flow scale. On the much larger scale, winds are now widely invoked to play a critical role for the interaction between the AGNs and the host galaxy, resulting in their coevolution (e.g., Ciotti et al. 2010, 2017; Ostriker et al. 2010; Choi et al. 2012; Eisenreich et al. 2017; Weinberger et al. 2017; Yoon et al. 2018, 2019; Yuan et al. 2018). For example, Yuan et al. (2018) recently comprehensively include feedback by wind and radiation from AGNs in both cold and hot feedback modes and find that wind plays a dominant role in both modes, although radiative feedback cannot be neglected.

Because the scales of wind launching and feedback differ by orders of magnitude, important questions then emerge: can wind produced on the accretion disk scale escape the gravitational potential of the central black hole, and how far can it propagate outward in the combined potential of the black hole and galaxy? What is the detailed dynamics of winds when they propagate outward, i.e., how do the velocity and density of winds change with radius?

For winds from hot accretion flows, answering these questions is somewhat pressing because recent cosmological simulations find that to overcome some serious problems in galaxy formation, e.g., reducing star formation efficiency in the most massive halos, winds launched from hot accretion flows must be invoked to interact with the interstellar medium on the galaxy scale (e.g., Weinberger et al. 2017). However, because a comprehensive study of the wind launching on the accretion flow scale was performed only recently (i.e., Y15), the large-scale dynamics of winds has not been investigated yet. This is the primary goal of our present work.

The crucial factor of studying the large-scale dynamics of winds is to employ correct boundary conditions because the hydrodynamical equations controlling the wind dynamics is a set of differential equations. Many works have studied winds from thin disks, either around young stars (e.g., Font et al. 2004; Bai et al. 2016) or black holes (e.g., Contopoulos & Lovelace 1994; Romanova et al. 1997; Proga et al. 2000; Proga 2003; Proga & Kallman 2004; Luketic et al. 2010; Waters & Proga 2012; Cao 2014; Clarke & Alexander 2016; Nomura & Ohsuga 2017; Waters & Proga 2018), and some works have already extended to very large radii. However, because these works fail to take into account the three driving mechanisms described above, they could not supply realistic boundary conditions. Therefore we need to revisit the large-scale dynamics with realistic boundary conditions.

In addition to these issues, another important factor is that we should take into account the gravitational potential of the host galaxy because we are interested in the dynamics of wind well beyond the accretion flow scale. The effect of the galaxy potential is hard to estimate without performing detailed calculations. Most previous works do not include this ingredient.

In the present paper, we aim at systematically studying the dynamics of winds launched from both hot accretion flows and cold thin disks. Realistic boundary conditions are analyzed and adopted, and the galaxy potential is included. The inner boundary is set at ∼10³r_g. Radiation force is neglected. This assumption is reasonable for wind from hot accretion flows because the radiation of an accretion flow is weak and the gas is fully ionized. For wind from a thin disk, however, radiation force is very likely to play an important role. In addition to simplifying calculation, we hope our thermal assumption can provide a lower limit to the dynamics of wind. The role of the magnetic field is unclear because the observational constraints on the magnitude and configuration are poor. In this paper we focus on the case without a magnetic field; the magnetic field will be taken into account in our next work.

The paper is structured as follows. In Section 2 we present analytical solutions of winds from hot accretion flows and thin disks. We start with an analogy with the solar wind to illustrate why we only look for transonic or supersonic solutions (Section 2.1). After a brief description of the model setup (Section 2.2), we present the gravitational potential of the galaxy (Section 2.3) and the equations we use (Section 2.4). In Section 2.5 we present a detailed discussion of how we adopt proper boundary conditions in the cases of cold and hot accretion flows. The results are described in detail in Section 2.6. In Section 3 we present 1D and 2D hydrodynamical numerical simulations and compare with the analytic solutions. We summarize and discuss our results in Section 4.

2. Analytical Solutions

2.1. Analogy with the Solar Wind

To study the large-scale dynamics of wind from accretion disks, the solar wind can be a valuable reference. A canonical model describing the dynamics of solar wind is the Parker model (Parker 1958, 1960), in which distinct types of solutions are found. Requiring the wind velocity to be subsonic in the near-Sun region, only two solutions are allowed, namely, the transonic solution and the subsonic solution. A transonic solution represents a transition from a subsonic to a supersonic flow by passing through the critical (sonic) point, and the radial velocity tends to increase monotonically with radius. A nonmonotonical dependence of the radial speed on heliocentric distance occurs when multiple critical points exist as a result of the super-radial expansion of the flow tube, for example (e.g., Yeh & Pneuman 1977; Cuperman et al. 1990; Li et al. 2011). However, a transonic solution is always possible whereby the solution chooses the proper sonic point from the available critical points. The subsonic solution, also called the breeze solution in solar wind problems, has its velocity decreasing at large distances and never passes across the critical point.

The subsonic solution has been discarded conventionally for solar winds mainly because (1) at large radii the solution yields a finite pressure and density, which cannot be matched to the interstellar gas properties; (2) in situ measurements of the near-Earth solar wind show that the wind speed far exceeds the local sound speed (e.g., Abbo et al. 2016); and (3) the subsonic solution is proved to be unstable to low-frequency acoustic perturbations (Velli 2001). Consequently, we restrict ourselves to transonic solutions that are physical for the disk wind problem, although the subsonic solutions are presented for comparison.

2.2. Model Setup

We solve for steady-state 1D hydrodynamical equations in spherical polar coordinates (r, θ, ϕ) following Abramowicz & Zurek (1981). The inner boundary is set to be r₀ = 10³r_g in order to obtain reliable boundary conditions from the small-scale simulations of Y15. Here the gravitational radius is defined as ${r}_{g}\equiv {GM}/{c}^{2}$ and M is the black hole mass (r_g ∼ 10⁻² pc for M = 10⁸ M_⊙). We prescribe wind to propagate along a direction with constant angle θ = 45° from the equatorial plane, although the specific angular momentum is nonzero in the model. This simplification is well justified in 2D simulations in this paper (Section 3.2). We have experimented with different θ angles and found that this did not alter the results qualitatively.

2.3. Gravitational Potential

In this work we deal with winds extending to galaxy scales so the gravitational potential of both the central black hole and the host galaxy should be taken into account:

$\begin{eqnarray}&&{\rm{\Phi }}={{\rm{\Phi }}}_{\mathrm{BH}}+{{\rm{\Phi }}}_{\mathrm{galaxy}}.\end{eqnarray} \tag{ 1 }$

We adopt the usual pseudo-Newtonian gravitational potential to describe the black hole potential (Paczyński & Wiita 1980),

$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{BH}}=-\displaystyle \frac{{GM}}{r-{r}_{s}},\end{eqnarray} \tag{ 2 }$

where ${r}_{s}\equiv 2{{GM}}_{\mathrm{BH}}/{c}^{2}=2{r}_{g}$ is the Schwarzschild radius.

Observations show that the stellar and dark matter components of galaxies are distributed so that their total mass profile is well described by a density distribution ρ ∝ r⁻² over a large radial range (e.g., Treu & Koopmans 2002, 2004; Rusin & Kochanek 2005; Gavazzi et al. 2007; Czoske et al. 2008; Dye et al. 2008). This leads to the simple assumption that the galaxy potential is treated as a flat rotation curve with a constant velocity dispersion parameter σ. Then, the difference of the potential of the host galaxy between r and r₀ is written

$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Phi }}}_{\mathrm{galaxy}}={\int }_{{r}_{0}}^{r}\displaystyle \frac{{\sigma }^{2}}{r}={\sigma }^{2}\mathrm{ln}r+C,\end{eqnarray} \tag{ 3 }$

where $r=\sqrt{{R}^{2}+{z}^{2}}$ denotes the spherical radius and C is a constant. We adopt a velocity dispersion of 200 $\mathrm{km}\,{{\rm{s}}}^{-1}$ , which is common to elliptical galaxies (Kormendy & Ho 2013).

2.4. Equations

The conservation of mass and energy are given by (Abramowicz & Zurek 1981)

$\begin{eqnarray}&&\rho {v}_{r}{r}^{2}=\dot{M},\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&\displaystyle \frac{1}{2}({v}_{r}^{2}+{v}_{\phi }^{2})+\displaystyle \frac{{c}_{s}^{2}}{\gamma -1}-\displaystyle \frac{{GM}}{r-{r}_{s}}+{\sigma }^{2}\mathrm{ln}r=E\end{eqnarray} \tag{ 5 }$

,where $\dot{M}$ denotes the mass flux of the wind, and Equation (5) is the Bernoulli integral with the constant Bernoulli parameter E. The rotational velocity v_ϕ is related to the specific angular momentum by l = v_ϕR. The sound speed is defined as ${c}_{s}^{2}\equiv \partial P/\partial \rho =\gamma K{\rho }^{\gamma -1}$ , where γ and K are the polytropic index and polytropic constant, respectively. We adopt γ = 4/3 throughout the analytical section and defer our discussion on the adoption of the value and its physical interpretation to Section 3.1.2. Substituting ρ by c_s leads us to rewrite Equation (4) as

$\begin{eqnarray}&&{c}_{s}^{\displaystyle \frac{2}{\gamma -1}}{v}_{r}{r}^{2}=\dot{M}^{\prime} .\end{eqnarray} \tag{ 6 }$

Note that here $\dot{M}^{\prime}$ differs from $\dot{M}$ by a constant coefficient. We drop the superscript prime for the rest of the derivations. After differentiating Equations (5) and (6), we arrive at

$\begin{eqnarray}&&\begin{array}{l}{r}^{3}({v}_{r}-{c}_{s})\left(\displaystyle \frac{{{dv}}_{r}}{{dr}}-\displaystyle \frac{2}{\gamma -1}\displaystyle \frac{{{dc}}_{s}}{{dr}}\right)\\ \quad =\,\displaystyle \frac{{l}^{2}}{{\sin }^{2}\theta }-{l}^{{\prime} 2}-{\sigma }^{2}{r}^{2}+2{r}^{2}{v}_{r}{c}_{s},\end{array}\end{eqnarray} \tag{ 7 }$

where ${l}^{{\prime} }={{GMr}}^{3/2}/(r-{r}_{s})$ . At the sonic point, the wind velocity equals the sound speed so that Equation (7) can be reduced to

$\begin{eqnarray}&&{c}_{s}^{2}={v}_{r}^{2}=\displaystyle \frac{1}{2{r}_{{\rm{c}}}^{2}}\left(-\displaystyle \frac{{l}^{2}}{{\sin }^{2}\theta }+{l}_{{r}_{{\rm{c}}}}^{{\prime} 2}+{\sigma }^{2}{r}_{{\rm{c}}}^{2}\right),\end{eqnarray} \tag{ 8 }$

where r_c is the locus of the sonic point. Inserting Equation (8) into Equations (5), we can solve for r_c as a function of $(E,{l}^{2})$ , or more easily, we can solve for the angular momentum

$\begin{eqnarray}&&\begin{array}{l}{l}^{2}(E,{r}_{{\rm{c}}})=\displaystyle \frac{2{R}_{{\rm{c}}}^{2}}{\tfrac{1}{2}-n}\\ \quad \cdot \,\left(E-\displaystyle \frac{\left(\tfrac{1}{2}+n\right)({l}^{{\prime} 2}+{\sigma }^{2}{r}_{{\rm{c}}}^{2})}{2{r}_{{\rm{c}}}^{2}}+\displaystyle \frac{{GM}}{{r}_{{\rm{c}}}-{r}_{s}}-{\sigma }^{2}\mathrm{ln}{r}_{{\rm{c}}}\right).\end{array}\end{eqnarray} \tag{ 9 }$

Similarly, inserting Equation (8) into (6) yields the angular momentum as a function of $(\dot{M},{r}_{{\rm{c}}})$

$\begin{eqnarray}&&{l}^{2}(\dot{M},{r}_{{\rm{c}}})=({l}_{{r}_{{\rm{c}}}}^{{\prime} 2}+{\sigma }^{2}{r}_{{\rm{c}}}^{2}-{\dot{M}}^{\displaystyle \frac{1}{2n+1}}2{r}_{{\rm{c}}}^{\displaystyle \frac{4n}{2n+1}}){\sin }^{2}\theta .\end{eqnarray} \tag{ 10 }$

Inspecting Equations (9) and (10), we find that four unknowns $l,E,\dot{M},{r}_{{\rm{c}}}$ are present in two equations, requiring us to specify two of these variables and to solve for the other two. We estimate l and E based upon the boundary conditions (Section 2.5). Once constants l and E are specified, the loci of the sonic points r_c can be computed via Equation (9) and then mass fluxes $\dot{M}$ via Equation (10), hence the velocity and temperature (or sound speed) profiles of the transonic solution via solving Equations (5) and (6).

Our differential equations have three variables, so usually, we should supply three boundary conditions to solve the equations, namely radial velocity (v_r0), rotational velocity ( ${v}_{\phi 0}$ ), and temperature (T₀). They in turn determine the values of E, l, and $\dot{M}$ . However, for a set of arbitrarily given boundary conditions the solution in general does not pass through the sonic point, i.e., it is not a transonic solution. The sonic condition provides an additional constraint that requires a specific combination of the three quantities at the boundary. In other words, fine-tuning is required to obtain a transonic solution, i.e., only two of them are free to choose for obtaining a transonic solution.

2.5. Boundary Conditions

In this section we discuss the boundary conditions employed for winds from hot accretion flows and thin disks, i.e., the values of ${v}_{r0},{v}_{\phi 0}$ , and T₀ at ${r}_{0}={10}^{3}{r}_{g}$ .

2.5.1. Hot Accretion Flows

The dynamics of hot accretion flows around black holes are well studied (see Yuan & Narayan 2014 for a review). Wind properties from hot accretion flows are investigated in Y15 based on 3D GRMHD simulations. It is found that winds can be produced from r ∼ 30r_g up to the outer boundary of the accretion flow, implying that winds at r₀ = 10³r_g are a combination of those originated from r ≤ 10³r₀. Furthermore, winds originating from various radii have different velocities, and they remain almost constant when the wind propagates outward, implying that the velocities at r₀ are diverse. The mass flux is proportional to the radius as ${\dot{M}}_{\mathrm{wind}}\propto {r}^{s}$ , with s ≈ 1, indicating that the wind properties are dominated by the locally generated ones rather than by components originating from smaller radii. The mass flux-weighted poloidal velocity is described by ${v}_{p}\approx 0.2{v}_{k}(r)$ , with v_k(r) the Keplerian velocity at r.

Y15 demonstrate that winds occupy the region of 0° ≲ θ ≲ 45° and inflows are in the region of 45° ≲ θ ≲ 90°. Figure 1 shows the azimuthal distributions of some wind properties, such as velocity and temperature, at 10³r_g obtained in Y15. Given the diversity of these quantities at different θ, we consider various boundary conditions at r = 10³r_g (Table 1). For example, the boundary conditions of the "hot_bhg_sup" model, v_r0 = 3v_k, T₀ = 10¹⁰ K, and ${v}_{\phi 0}=0.7{v}_{k0}$ , correspond to winds at θ = 20°. Note that this set of physical quantities will result in a Mach number that exceeds unity at the boundary, thereby a corresponding supersonic solution is obtained, as we discuss later. Most of the models for hot winds are assigned with a rotational velocity of $0.5{v}_{k0}$ . We have confirmed that moderate deviation from this value ( $0.4\mbox{--}0.6{v}_{k0}$ ) did not alter the results qualitatively.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Physical quantities of winds as a function of θ at r = 10³r_g from the 3D GRMHD simulation of black hole accretion in Y15. Top panels: density, temperature, and Mach number. Bottom panels: radial, meridional, and rotational velocities normalized by Keplerian velocity ${v}_{k}({r}_{0})$ . The value of the density is shown in code units. The density and temperature are time-averaged values, while all three velocity components are taken from one snapshot.
Download figure:
Standard image High-resolution image

**Figure 1.** Physical quantities of winds as a function of θ at r = 10³r_g from the 3D GRMHD simulation of black hole accretion in Y15. Top panels: density, temperature, and Mach number. Bottom panels: radial, meridional, and rotational velocities normalized by Keplerian velocity ${v}_{k}({r}_{0})$ . The value of the density is shown in code units. The density and temperature are time-averaged values, while all three velocity components are taken from one snapshot.
Download figure:
Standard image High-resolution image

2.5.2. Thin Accretion Disks

Wind properties from thin disks are different from those from hot accretion flows. Because we still lack a comprehensive theoretical model for the production of winds from thin disks, we incorporate both observational and numerical simulation results to determine the wind properties at the boundary.

Observationally, a sample of over 50 AGNs is collected from Suzaku observations to study the properties of disk winds (Gofford et al. 2015). These winds are detected in a spatial range of 10²–10⁴r_s from the central black hole. A power-law relation between the AGN bolometric luminosity and the wind velocity is found,

$\begin{eqnarray}&&{v}_{w}=2.5\times {10}^{4}\left(\displaystyle \frac{{L}_{\mathrm{BH}}}{{10}^{45}\ \mathrm{erg}\ {{\rm{s}}}^{-1}}\right)\,\mathrm{km}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$

When we assume ${L}_{\mathrm{BH}}={10}^{44}\mbox{--}{10}^{45}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ (0.01–0.1L_Edd for a black hole mass of M_BH = 10⁸ M_⊙), this implies wind velocities in the range of ${v}_{w}=0.25{v}_{k}-2.5{v}_{k}$ or 0.008c–0.08c, where c is the speed of light. We adopt values that consecutively increase from 0.2v_k to 2v_k as our boundary conditions of different models. However, given the diversity of various observations of quasar wind, values beyond this range are also tested.

Because it is difficulty to obtain the temperature and rotational velocity of winds from observations, we use numerical simulations to determine the values of these quantities. Proga et al. (2000) have performed hydrodynamical simulations on radiation line-driven winds from geometrically thin and optically thick accretion disks of luminous AGNs. Their results suggest that the wind temperatures at 10³r_g range from 10⁴ to 10⁶ K. Motivated by the above discussions, we set the boundary conditions of the fiducial model of thin-disk winds to v_r0 = v_k, ${v}_{\phi 0}=0.5{v}_{k}$ , and T₀ = 10⁵ K. Various values of T₀ and v_r0 are adopted for other models (Table 1). The adoption of a relatively high rotational velocity of ${v}_{\phi 0}=0.5{v}_{k0}$ for thin-disk winds comes from the detailed radiation line-driven wind simulations (W.X. Wang et. al. 2020, in preparation) at 10³r_g. Given that thin-disk winds may be predominantly driven from the inner region of the disk, the angular velocity at r₀ = 10³r_g could be low, hence we explore a parameter space in which ${v}_{\phi 0}$ is allowed to vary from 0.01 to 0.5v_k0.

2.6. Solutions

We look for transonic solutions following the approach presented in Section 2.4 with boundary conditions given in Section 2.5. We begin with a general analysis of the topology to the solutions in Section 2.6.1, then present the detailed solutions for hot accretion flows in Sections 2.6.2 and 2.6.3, and those for thin disks in Section 2.6.4. Special attention is paid to the influence of including galaxy potential on the solutions.

2.6.1. Solution Topology

Based on Equation (9), Figure 2 shows the square of angular momentum l² as a function of the locus of the critical point r_c at a specific Bernoulli parameter E. The solid and dotted curves correspond to model "hot_bh_t" and "hot_bhg_tb," respectively. On each curve, a specification of angular momentum l marks the loci of the critical points, shown as intersections of the gray dashed curve or the solid and dotted curves.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The square of angular momentum l² as a function of the locus of the critical point r_c at a specific Bernoulli parameter E. The solid curve corresponds to "hot_bh_t" (Table 1), which only includes the black hole potential, and the dotted curve corresponds to "hot_bhg_tb", which includes the potential of the black hole and the galaxy. The horizontal dashed curve shows the angular momentum l² adopted in "hot_bh_t" and "hot_bhg_tb" for easy inspection. The intersections with the dashed curve reveal the loci of critical points.
Download figure:
Standard image High-resolution image

Note that not all the critical points correspond to transonic solutions. Only saddle points are related to transonic solutions, while the vortex points only have mathematical meaning. These two types of critical points are usually denoted by "X" point and "O" point in the studies of hydrodynamical winds (Liang & Thompson 1980; Lu & Abramowicz 1988). The classic Bondi accretion has one critical point, which is a saddle point (Bondi 1952). When taking angular momentum into account, three critical points are possible (Abramowicz & Zurek 1981). The point closest to the black hole corresponds to the accretion of material onto the black hole, which requires the gas to pass the sound speed eventually. Thus it is a saddle point and related to the transonic solution. The second point is a vortex point, and the third point is a saddle point that is relevant to the study of the wind here. In Figure 2 the intersection of the dashed curve and the black curve and the first intersection of the dashed curve and the dotted curve correspond to the critical point of the wind. The first and second critical points of both the solid and dotted curves are located at radii r < r₀, thus are not seen in the figure.

The second intersection between the dashed curve and the dotted curve is an additional point to these three critical points. Its presence is due to the inclusion of galaxy potential. The locus of this critical point is determined by the boundary condition. For example, we find that the loci of the "hot_bhg_ta" model are much larger than that of the "hot_bhg_tb" model. This point is of vortex type, as we discuss in the following paragraph.

The topology of solutions is further illustrated in Figure 3, which shows the Mach number of the wind as a function of radius. The red curve in the left panel is for the "hot_bh_t" model, which possesses the topology of a saddle point. It corresponds to the black curve in Figure 2 with the parameters listed in Table 1. The red curve in the right panel of Figure 3 is for the "hot_bhg_tb" model. This model includes the galaxy potential and therefore possesses both a saddle point and a vortex point, which corresponds to the dotted curve in Figure 2. In both panels, the black curves are drawn by slightly adjusting the Bernoulli parameter E to deviate from the transonic solution, but with $\dot{M}$ and l fixed. The black curves above and below the red curve represent supersonic and subsonic solutions, respectively.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Mach number of the wind as a function of radius. Different curves in each panel correspond to the same mass flux $\dot{M}$ and angular momentum l, but different Bernoulli parameter E. The red curve in the left panel is for the "hot_bh_t" model, and the red curve in the right panel for the "hot_bhg_tb" model. Subsonic and supersonic solutions are present below and above the red curves.
Download figure:
Standard image High-resolution image

**Figure 3.** Mach number of the wind as a function of radius. Different curves in each panel correspond to the same mass flux $\dot{M}$ and angular momentum l, but different Bernoulli parameter E. The red curve in the left panel is for the "hot_bh_t" model, and the red curve in the right panel for the "hot_bhg_tb" model. Subsonic and supersonic solutions are present below and above the red curves.
Download figure:
Standard image High-resolution image

2.6.2. Winds from Hot Accretion Flow with Black Hole Potential Alone

Figure 4 shows the radial profiles of velocity, density, Mach number, and temperature of the representative transonic (solid line) and subsonic (dashed line) solutions. For other transonic and subsonic solutions the results are similar once their boundary conditions are in the realistic range. The supersonic solution is not shown here because its profiles share great similarity with the transonic solution (see the dotted curves in Figure 6).

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Radial profiles of radial velocity, Mach number, temperature, and density of the subsonic solution "hot_bh_s" (dashed) and transonic solution "hot_bh_t" (solid) of a wind launched at ${r}_{0}={10}^{3}{r}_{s}$ , with the detailed boundary parameters listed in Table 1.
Download figure:
Standard image High-resolution image

**Figure 4.** Radial profiles of radial velocity, Mach number, temperature, and density of the subsonic solution "hot_bh_s" (dashed) and transonic solution "hot_bh_t" (solid) of a wind launched at ${r}_{0}={10}^{3}{r}_{s}$ , with the detailed boundary parameters listed in Table 1.
Download figure:
Standard image High-resolution image

It is clear from the transonic and supersonic solutions that the wind can escape from the black hole and propagate to very large radius. For the transonic solutions, we find that the Mach number increases rapidly with radius, and the radial velocity only slightly increases. The radial profiles of density and temperature can be derived from the continuity equation and the polytropic relation. They can be roughly approximated by

$\begin{eqnarray}&&{v}_{r}\approx {v}_{r0},\,\,\rho \approx {\rho }_{0}{(r/{r}_{0})}^{-2},\,\,T\approx {T}_{0}{(r/{r}_{0})}^{-2(\gamma -1)}.\end{eqnarray} \tag{ 12 }$

We may imagine that the radial velocity decreases as a result of the gravitational force of the black hole. To understand why the radial velocity does not decrease with radius, we decompose the Bernoulli parameter E into individual components (Equation (5)) along the wind trajectory; the result is shown in the left plot of Figure 5. It demonstrates that the increase in gravitational energy is mainly compensated for by the reduction of the specific enthalpy as well as the rotational energy. At large radii, both the enthalpy and potential energy nearly vanish, and the total energy is dominated by the radial kinetic energy. In other words, the gas pressure gradient and centrifugal forces overcome the gravitational force of the black hole and do work to accelerate the wind.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Decomposition of the Bernoulli constant E into individual components for models "hot_bh_t" (left) and "'cold_fid" (right), including the specific radial kinetic energy KE_r $(1/2{v}_{r}^{2})$ , the specific rotational energy KE_ϕ $(1/2{v}_{\phi }^{2})$ , the specific enthalpy $h={c}_{s}^{2}/(\gamma -1)$ , and the specific potential energy ${{\rm{\Phi }}}_{\mathrm{BH}}=-{GM}/(r-{r}_{s})$ . For the right plot, a curve showing the ${\sigma }^{2}\mathrm{ln}r$ term in Equation (5) is also shown. It almost overlaps the enthalpy curve.
Download figure:
Standard image High-resolution image

**Figure 5.** Decomposition of the Bernoulli constant E into individual components for models "hot_bh_t" (left) and "'cold_fid" (right), including the specific radial kinetic energy KE_r $(1/2{v}_{r}^{2})$ , the specific rotational energy KE_ϕ $(1/2{v}_{\phi }^{2})$ , the specific enthalpy $h={c}_{s}^{2}/(\gamma -1)$ , and the specific potential energy ${{\rm{\Phi }}}_{\mathrm{BH}}=-{GM}/(r-{r}_{s})$ . For the right plot, a curve showing the ${\sigma }^{2}\mathrm{ln}r$ term in Equation (5) is also shown. It almost overlaps the enthalpy curve.
Download figure:
Standard image High-resolution image

A constant or slightly increasing radial velocity resembles the velocity profile of wind on small accretion flow scales. On this scale, Y15 have found that the wind velocity almost remains constant along their trajectories, and it is also the specific enthalpy that compensates for the increase in potential energy. However, the Bernoulli parameter of the wind is found to increase along the wind trajectory rather than remaining constant. This is because the accretion flow is strongly turbulent, while conservation of the Bernoulli parameter holds only for a strictly steady and inviscid flow. Our large-scale wind is well beyond the accretion disk scale, so that turbulence can reasonably be assumed to be absent and the Bernoulli parameter should be constant.

The radial profiles of the subsonic solution ("hot_bh_s") strongly deviate from the transonic solution, as is shown by the dashed curves in Figure 4. The radial velocity rapidly decreases with radius, while the temperature and density drop slowly and tend to be constant at large distances. This subsonic solution is obtained by adjusting the Bernoulli parameter E and keeping the mass flux $\dot{M}$ and angular momentum l the same as in "hot_bh_t." In this case, the sonic point condition Equation (8) is not satisfied. With a slight adjustment of E, the physical parameters at the boundary deviate from the parameters adopted in the transonic solution "hot_bh_t." As we have explained in Section 2.4, this will result in the failure of obtaining a transonic solution, implying that a transonic solution requires a specific combination of physical parameters at the boundary.

The fine-tuning of the transonic solution may lead to the conclusion that subsonic winds should be regarded as more physically realizable in nature than the supersonic winds. However, the evidence from the detection of solar winds contradicts this conjecture; the solar winds that have been detected are always supersonic at the Earth's orbit. The classic Parker wind model argues that the nonvanishing pressure in the subsonic solution prevents the wind from smoothly transitioning to the interstellar medium at infinity, thereby excluding the subsonic solutions from physical solutions. Moreover, Velli (1994) states that the subsonic solutions are unstable to low-frequency acoustic perturbations, leaving the transonic solution the only plausible solution to connect to infinity. Our 1D simulations prove the prevalence of transonic solutions by the facts that only transonic solutions are found, whereas no subsonic solutions are ever obtained (Section 3.1.1).

2.6.3. Winds from Hot Accretion Flow with Both Black Hole and Galaxy Potentials

As the wind propagates toward large radii, the galaxy potential should be included. Figure 6 shows the radial profiles of two transonic solutions "hot_bhg_ta" and "hot_bhg_tb," and a supersonic solution "hot_bh_sup." The subsonic solution "hot_bhg_s" is not shown in the figure because it resembles the profiles of the "hot_bh_s" model. As we have mentioned in Section 2.5.1, the wind at r₀ is a combination of components originating from various radii satisfying r ≲ 10³r_g, hence having various velocities. The boundary conditions of the three solutions shown here mainly differ by their radial velocities, with the supersonic solution "hot_bhg_sup" having the highest velocity and "hot_bhg_tb" having the lowest one (Table 1). The velocity in the "hot_bhg_tb" model is extremely low and not realistic because we deliberately chose this value to manifest the effect of the galaxy potential. The radial velocity in "hot_bhg_ta" is close to the mass flux-weighted radial velocity of the wind obtained in Y15.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Radial profiles of the radial velocity, Mach number, temperature, and density of the transonic solutions" hot_bhg_ta," "hot_bhg_tb" and the supersonic solution "hot_bh_sup," with detailed boundary parameters as listed in Table 1.
Download figure:
Standard image High-resolution image

The radial profiles of the transonic solution "hot_bhg_ta" (solid line) are similar to those of the transonic solution of the black hole-potential-only case, i.e., the solid line in Figure 4. The profiles of radial velocity, density, and temperature are also well approximated by Equation (12). This is also the case of the supersonic solution "hot_bh_sup." This result indicates that the effect of the galaxy potential is limited when "realistic" boundary conditions are employed. We have estimated the largest distance that the wind can propagate outward for the "hot_bhg_ta" model and found that it is well beyond the scale of the galaxy. In reality, as a result of the interaction between the wind and the interstellar medium, this estimate should be regarded as an upper limit.

The effect of the galaxy potential is significant if the radial velocity of the wind is low at the boundary, as illustrated by the dashed curve in Figure 6, which delineates the behavior of the solution "hot_bhg_tb." The wind has been accelerated first, with the radial velocity and the Mach number increasing. Beyond $r={10}^{6}{r}_{g}$ , the acceleration stops and the wind begins to decelerate beyond 10⁷r_g. This is because the adopted radial velocity at r₀ is so low that the galaxy potential plays a relatively much more important role. With the increase in radius, the galaxy potential energy becomes higher, thus the kinetic energy has to be decreased according to the conservation of total energy. From this argument, we expect that the realistic transonic solution "hot_bhg_ta" shares the same behavior as "hot_bhg_tb" when our calculation domain is significantly extended because the increasing galaxy potential will eventually slow the wind down.

The effect of the galaxy potential on wind dynamics depends on the magnitudes of the gravitational force corresponding to this potential and also on the velocity of the wind itself. A low wind velocity and/or strong force will lead to prominent variations in the wind velocity. In our case, although the force from the galaxy can be estimated, the amplitude of the wind velocity requires detailed calculations with realistic boundary conditions. From the results presented above, we find that for most cases the wind velocity is not very low, and the effect of the galaxy potential is limited.

From these results, we can obtain the following propagation scenario of wind that is launched from a hot accretion flow. The winds close to the axis have the highest velocity at the boundary, therefore they will propagate to the largest distance in the galaxy, while the winds close to the equatorial plane have the lowest velocity and therefore will stop at smaller distances.

2.6.4. Winds from Thin Disks

From Table 1, we can see that the boundary conditions of winds from a thin disk are quite different from that of a hot accretion flow. The temperature and consequently the Bernoulli parameter are much lower in this case. In addition, for all the models of thin-disk winds, the wind is already supersonic at 10³r_g, which is mainly attributed to the low temperature of the wind. This implies that all the wind solutions can be realized in nature.

Table 1. Parameters of Analytic Solutions

	Potential	E	r_c	r₀	T₀	v_r0	v_ϕ0	Mach₀	Branch
			(r_g)	(r_g)	(K)	(v_k0)	(v_k0)
Hot Accretion Flow

hot_bh_s	bh	1.76 × 10⁻⁴		10³	1.32 × 10⁹	0.21	0.5	0.44	sub
hot_bh_t	bh	1.66 × 10⁻⁴	3.83 × 10³	10³	1.30 × 10⁹	0.22	0.5	0.46	trans
hot_bhg_s	bh+g	1.79 × 10⁻⁴		10³	1.32 × 10⁹	0.21	0.5	0.44	sub
hot_bhg_ta	bh+g	1.69 × 10⁻⁴	3.86 × 10³	10³	1.30 × 10⁹	0.22	0.5	0.46	trans
hot_bhg_tb	bh+g	8.00 × 10⁻⁶	5.65 × 10⁵	10³	1.12 × 10⁹	5.8 × 10⁻⁴	0.5	1.33 × 10⁻³	trans
hot_bhg_sup	bh+g	1.31 × 10⁻²		500	10¹⁰	3.0	0.7	2.29	sup
Thin Disk

cold_fid	bh+g	−1.09 × 10⁻⁴		10³	10⁵	1.0	0.5	241	sup
cold_vr0.2	bh+g	−7.92 × 10⁻⁴		10³	10⁵	0.2	0.5	48.2	sup
cold_vr2	bh+g	7.54 × 10⁻³		10³	10⁵	2.0	0.5	482	sup
cold_vphi0.3	bh+g	−2.2 × 10⁻⁴		10³	10⁵	1.0	0.3	241	sup
cold_vphi0.1	bh+g	−2.8 × 10⁻⁴		10³	10⁵	1.0	0.1	241	sup
cold_vphi0.01	bh+g	−2.9 × 10⁻⁴		10³	10⁵	1.0	0.05	241	sup
cold_T4	bh+g	−1.10 × 10⁻⁴		10³	10⁴	1.0	0.5	762	sup
cold_T6	bh+g	−1.09 × 10⁻⁴		10³	10⁶	1.0	0.5	76.2	sup
cold_T4vr2	bh+g	2.02 × 10⁻³		10³	10⁴	2.0	0.5	76.2	sup

Download table as: ASCII Typeset image

Following the same approach as for the hot winds in studying wind dynamics, we find that in general the wind from a thin disk cannot propagate as far as in the case of hot accretion flows. The stop radii of different models depend on the value of the Bernoulli parameters, with a larger E corresponding to a larger "stop radius." Figure 7 shows the dynamics of two examples of models with different radial velocities taken from Table 1, "cold_fid" and "cold_vr2." The difference of the velocity causes different Bernoulli parameters, with the Bernoulli parameter of the "cold_fid" model being negative and that of the "cold_vr2" model positive. We can see from the figure that the wind in the "cold_fid" model cannot even escape from the black hole gravity and stops at $\lesssim {10}^{4}{r}_{g}$ . As we have emphasized before, because radiation and magnetic forces are neglected in our calculation, in reality, the stop radius should be larger. The value of the stop radius is sensitive to the value of the Bernoulli parameter at the boundary. For the wind in the "cold_vr2" model whose velocity at the boundary is only two times higher, it stops at a much larger radius that lies beyond the radial range shown in the figure. Given these results, the fact that winds from quasars have been widely observed as far as ∼15 kpc from the black hole (e.g., Liu et al. 2013) implies that either radiation and magnetic forces must continue to accelerate winds at large radii, or the velocity of the wind at the boundary must be relatively high. This can potentially be checked by examining observational data.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Radial profiles of the radial velocity, Mach number, temperature, and density of thin-disk winds for models "cold_fid" and "cold_vr2."
Download figure:
Standard image High-resolution image

The value of the stop radius can be roughly estimated as follows. From Equation (5) we can deduce that the radius where the wind stops is determined by equating the sum of two potential energy terms to the Bernoulli parameter because the kinetic energy approaches zero there while the enthalpy is always negligible for the cold wind. As an illustrative example, we draw the right plot of Figure 5, which shows the energy decomposition of the "cold_fid" model. At large radii, all the terms in Equation (5) are close to zero except for the Bernoulli parameter and the black hole potential terms. These two nonzero terms are equal to each other at ≲10⁴r_g, which is exactly the radius where the wind stops. We also vary the radial and rotational velocities as well as the temperature at the boundary for thin-disk winds (Table 1). Slower radial and rotational velocity and lower temperature result in a smaller stop radius because the Bernoulli parameter is consequently smaller. Figure 8 shows that when ${v}_{\phi 0}$ is lower than $0.1{v}_{k0}$ , its impact on the wind dynamics is weak because other parameters become dominant in determining the Bernoulli constant.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Radial profiles of the radial velocity, Mach number, temperature, and density profiles of different v_ϕ0 for thin-disk winds (Table 1). Curves denote angular velocities of $0.5{v}_{k0}$ (cold_fid), $0.3{v}_{k0}$ (cold_vphi0.3), $0.1{v}_{k0}$ (cold_vphi0.1), and $0.01{v}_{k0}$ (cold_vphi0.01) at the inner boundary.
Download figure:
Standard image High-resolution image

**Figure 8.** Radial profiles of the radial velocity, Mach number, temperature, and density profiles of different v_ϕ0 for thin-disk winds (Table 1). Curves denote angular velocities of $0.5{v}_{k0}$ (cold_fid), $0.3{v}_{k0}$ (cold_vphi0.3), $0.1{v}_{k0}$ (cold_vphi0.1), and $0.01{v}_{k0}$ (cold_vphi0.01) at the inner boundary.
Download figure:
Standard image High-resolution image

These calculation results indicate that if the winds produced from the accretion disk have different components with different values of the Bernoulli parameter, as is likely the case, they will stop at different radii. This is similar to the case of wind from hot accretion flows, except that the stop radius should be systematically smaller because of the lower enthalpy of the cold wind. Those with a small stop radius, i.e., the so-called "failed winds," are invoked to explain the origin of the broad line region of AGNs by Giustini & Proga (2019).

Another question is how the wind velocity changes with radius when the winds propagate outward. This information is useful to study AGN feedback because numerical simulations of feedback that have different spatial resolutions need to adopt the wind velocity at different radii. We can see from Figure 7 that before the wind stops, its radial velocity only slightly decreases at the beginning, and then almost remains constant throughout the radius. This result indicates that during the outward propagation of the wind, until it is close to the "stop radius," the rotational energy of the wind almost exactly compensates for the gravitational potential energy. This result is very similar to the results of the transonic solution "hot_bhg_ta" shown in Figure 6. All the physical quantities can therefore be approximated by Equation (12).

3. Numerical Simulations

Using the Athena++ code, we perform 1D and 2D numerical simulations to examine the time-dependent solutions of hydrodynamical winds and compare with the analytical solutions. The conservation laws of mass, momentum, and energy in their conservative form read

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

$\begin{eqnarray}&&\displaystyle \frac{(\partial \rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{v}}{\boldsymbol{v}}+P{\boldsymbol{I}}\right)=-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot \left[(E+P){\boldsymbol{v}}\right]=0,\end{eqnarray} \tag{ 15 }$

where the total energy density is given by $E=P/(\gamma -1)+\rho {\rm{\Phi }}+\rho {v}^{2}/2$ and ${\boldsymbol{I}}$ is the identity tensor.

3.1. 1D Simulations

We first discuss the initial and boundary conditions and the polytropic index, two factors that deserve exceptional care, before we present the simulation results.

3.1.1. Initial and Boundary Conditions

For clarity, we start with 1D simulations without angular momentum (v_ϕ = 0). Thus, three variables must be specified at the inner radial boundary, namely, pressure, density, and radial velocity. For the first attempt, $P,\rho ,{v}_{r}$ are fixed in the ghost zone during the course of the simulation, which results in unphysical solutions. Discontinuities emerge in the first grid cell of the active zone, although it is connected with a smooth transonic solution at larger radii.

These solutions are interpreted as unphysical because the discontinuity clearly violates the conservation of mass. Mathematically speaking, this unphysical behavior arises because the number of dependent variables prescribed at the boundary exceeds the number of ingoing characteristics. Here we mean by ingoing that the characteristics point toward the computational domain (for details, see, e.g., Wu & Hu 1984; Grappin et al. 2000). For the 1D hydrodynamical problem, three characteristics exist that are related to the entropy wave (λ = v_r), the sound waves propagating forward and backward ( $\lambda ={v}_{r}+{c}_{s}$ , $\lambda ={v}_{r}-{c}_{s}$ ). Here, λ represents the eigenvalue and ${c}_{s}^{2}=\gamma P/\rho$ is the square of the adiabatic sound speed. For the injection of subsonic winds at the inner radial boundary, the entropy wave and the forward sound wave have positive eigenvalues (speeds), while the backward sound wave always moves in the opposite direction to the wind. Thus, only two variables can be arbitrarily prescribed at the inner boundary, with the third one being adjusted itself.

We adopt time-dependent boundary conditions for which the pressure and density are fixed while the radial velocity is allowed to be adjusted at each time step. Specifically, only at the first time step, we set the values of three quantities at both the ghost zone and the first active grid the values given in Table 1. At later times, only the pressure and the density are specified while the radial velocity is determined via ensuring mass conservation between the ghost zone and the first active grid cell, i.e., ${v}_{r}^{{\prime} }={\rho }_{0}{v}_{r0}{r}_{0}^{2}/{\rho }^{{\prime} }{r}^{{\prime} 2}$ , where the subscript "0" denotes quantities in the first active grid cell and primed quantities are in the ghost cell.

3.1.2. Polytropic Index

Mathematically, it can be proved that the set of equations comprised of the conservation laws of mass and momentum and the polytropic relation has transonic solutions when γ < 5/3. This criterion is obtained by requesting the critical point to be a saddle point, and is equivalent to having two eigenvalues of opposite signs for the differential equation. Moreover, outflowing gas always transits from subsonic to supersonic, thereby the slope of the radial velocity or the Mach number must be positive. This further confines the polytropic index to γ < 3/2 (Parker 1963). Our 1D simulations confirm these conclusions.

Physically, the polytropic index γ = 5/3 describes the adiabatic flows. A value of γ < 5/3 may result from the intricate interplay among thermal conduction, heating, and cooling. For instance, the rather flat radial profiles of the electron temperature in the near-Sun solar wind are usually attributed to the fact that the electron thermal conduction is rather efficient at temperatures exceeding ∼10⁶ K.⁴ As a result, empirical values of γ ≈ 1.1 were frequently adopted (see, e.g., Usmanov et al. 2000, and references therein). On the other hand, γ was empirically determined to be ∼1.46 in interplanetary space (e.g., Totten et al. 1995), indicating the existence of some weak but non-negligible heating (e.g., Marsch 2006).

Note that the criterion of γ < 3/2 holds only when the angular momentum of the wind is omitted. When the angular momentum is taken into account, the γ < 3/2 constraint is not necessarily valid. The inclusion of angular momentum complicates the set of equations, and there is no simple criterion on the polytropic index under this condition. We have tested a set of polytropic indices in cases without angular momentum and inspect if these γ values can result in transonic solutions when angular momentum is included. It is found that the angular momentum further confines the range of γ in which it is possible to obtain the transonic solution.

3.1.3. Simulation Results

In the simulation setup, we adopt the same boundary conditions as in the analytic solutions. The exact way has been described in Section 3.1.1. The wind is injected into the ghost zone cells, and the standard outflow boundary condition is employed at the outer boundary. We find that the transonic solution is the only steady solution regardless of which initial values of the radial velocity are specified. The subsonic solution that we obtain in the analytical solutions is not present in the numerical simulations. The explanation of its absence is based on the stability analysis (Velli 1994), which demonstrates that the subsonic solution is unstable to low-frequency acoustic waves. The power-law fittings of transonic solutions obtained in numerical simulations with rotational velocity are fully consistent with the analytical solutions.

3.2. 2D Simulations

Winds at different angles θ may interact with each other, which affects their dynamics. This motivates us to perform 2D axisymmetric simulations to mimic more realistic conditions. The polytropic index γ = 4/3 is adopted throughout, and the boundary conditions are employed from the small-scale 3D GRMHD simulation of black hole hot accretion flows (Y15) at a radius r = 10³r_g (Figure 1). The physical quantities are averaged over the azimuthal angle ϕ. We only take the time average for density and temperature over the time interval when the system reaches steady state. For velocities, we directly adopt values in one snapshot because the flow is turbulent in these simulations and the radial velocity frequently flips sign, so averaging over time will lead to a significant underestimate because positive and negative values cancel each other out. We emphasize that the data from θ ∈ [45°, 90°] are unphysical. The gas in this region is not wind but outflowing gas due to the outward angular momentum transport of the accretion flow. We include this region in 2D simulations but do not consider it when the simulation data are analyzed.

Figure 9 shows the radial profiles of the radial velocity, Mach number, temperature, and density of the wind at different θ. The wind reaches steady state up to r ∼ 3×10⁴r_g. The plot shows two types of behavior that are analogous to the analytical solutions. The first are flows at θ = 30° with a subsonic radial velocity at the inner boundary. They are accelerated to become supersonic. This resembles the transonic analytical solution such as "hot_bh_t" or "hot_bhg_ta." For easy comparison, we show the analytic solution of "hot_bhg_ta" by dashed gray curves. The profiles of the density, temperature, and velocity are also well approximated by Equation (12). Moving closer to the polar axis, supersonic solutions with radial velocity being supersonic at the boundary are shown, which correspond to "hot_bhg_sup." Its density, radial velocity, and temperature profiles also match the analytic solutions well. The similarity between 1D and 2D simulations indicates that the gas motions at different θ interact little with each other, hence a 1D simulation is a good approximation for wind dynamics. The results are also consistent with the conclusion of Y15 that the wind is largely laminar. We also find that the wind trajectories are well approximated by straight lines, which justifies the assumption adopted in our analytical solutions.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Radial profiles of the radial velocity, Mach number, temperature, and density profiles at several different θ angles in the 2D simulation. Colors denote angles of 5° (black), 10° (blue), 20° (dark green), and 30° (light green). The analytic solution of model "hot_bhg_ta" is shown by dashed curves.
Download figure:
Standard image High-resolution image

4. Summary and Discussions

Much work has been conducted to study the launch and acceleration of disk winds on the scale of black hole accretion flows. It remains unclear, however, how far the wind can propagate, and what the large-scale dynamics of these winds is. This is especially the case for winds from hot accretion flows because the study of the hot wind is relatively new and winds are invoked to play an important role in AGN feedback. As the first paper of this series to answer these questions, in this paper we study the large-scale dynamics of wind under the most basic thermal assumption. The role of the magnetic field will be investigated in our second paper. In this work, both analytical study and numerical simulation are performed. The inner boundary is set to be r₀ = 10³r_g. The boundary conditions of the wind, which play a crucial role in determining the wind dynamics, are discussed and presented in Table 1 (for analytical models and 1D simulations) and Figure 1 (for 2D simulations). For hot accretion flows, the wind production on the accretion scale has been well studied by 3D GRMHD simulations, so these boundary conditions are taken from these simulations. For thin disks, theoretical studies that take into account all physical driving mechanisms are still lacking because of technical difficulties, so we combine theoretical and observational studies to obtain realistic boundary conditions. Specifically, the velocity of the wind is taken from observations and the temperature and rotational velocity are adopted from numerical simulations. The potential energy of the black hole as well as the host galaxy are taken into account. We focus on transonic and supersonic solutions because they are physical solutions that can be realized in nature. Mathematically, subsonic solutions also exist, but the solution is unphysical, partly because the subsonic solution is unstable.

For winds from hot accretion flows, we find that winds can propagate to very large radii. Transonic, supersonic, and subsonic analytical solutions have all been obtained, depending on the specified boundary conditions. The physically viable transonic and supersonic solutions show similar radial profiles of dynamical quantities (Figures 4 and 6), and they are approximated by Equation (12). We find that when the galaxy potential is included, the solution does not deviate much from models that exclude it. The gravitational force exerted by the host galaxy is not able to modify the wind poloidal velocity significantly. This is mainly because the wind velocity from our detailed calculation is very high compared to the acceleration by the galaxy. The radial velocity of the wind remains constant with increasing radius because the enthalpy and rotational energy of the winds almost exactly compensate for the potential energy when the winds propagate outward (Figure 5).

For wind from thin disks, we find that the wind usually stops at a smaller radius than the wind from hot accretion flows. This is mainly because the temperature of the wind is much lower and accordingly, the Bernoulli parameter is smaller. The stop radius is determined by equating the Bernoulli parameter to the potential energy because all other terms in the Bernoulli equation can be neglected at the stop radius. Before the stop radius, however, as with the hot wind, the radial velocity of the cold wind also roughly remains constant and the radial profiles of the physical quantities can also be roughly described by Equation (12). During the propagation of the cold wind, different from the hot wind, only the rotation energy compensates for the increase in gravitational potential energy because the enthalpy is negligibly low.

We have performed 1D and 2D numerical simulations to study the dynamics of the hot wind and compared this with analytical solutions. When we started with a subsonic solution, we found that one of three characteristics had a negative velocity, which led to wave propagation out of the simulation domain. Hence, only two variables of $P,\rho$ and v_r can be given at the inner radial boundary, and the other variable should be self-consistently determined to satisfy the compatibility requirements. This is realized through time-dependent boundary condition where the radial velocity is allowed to be adjusted at each time step by ensuring mass conservation. We only find the transonic and supersonic solutions; the subsonic solution found in the analytical solution is not present because the solution is unstable. The profiles of the dynamical quantities for the transonic and supersonic solutions are well consistent with those from analytical solutions.

In this study, we have not taken into account the interaction between winds and the interstellar medium (ISM). An interesting question is then how far the wind can propagate when the ISM is included. To estimate the largest distance, we equate the ram pressure of the wind to the thermal pressure of the ISM. The thermal pressure of the wind can be ignored because at large radii, it is supersonic. We adopt the observed values of mass flux and velocity of thin-disk winds (Gofford et al. 2015). The mass flux is $2.5\,{M}_{\odot }\ {\mathrm{yr}}^{-1}$ and the wind velocity is ${v}_{r}=1.5{v}_{{\rm{K}}}$ at $r={10}^{3}{r}_{g}$ for a black hole mass ${M}_{\mathrm{BH}}={10}^{9}\,{M}_{\odot }$ and ${L}_{\mathrm{BH}}/{L}_{\mathrm{Edd}}=0.1$ . The solid angle of the wind is assumed to occupy ∼30% of the whole sphere (Tombesi et al. 2011). The terminal distance of the wind is estimated to be

$\begin{eqnarray}&&{D}_{\mathrm{term}}\approx 65\times {\left(\displaystyle \frac{{10}^{-3}{\mathrm{cm}}^{-3}}{{n}_{\mathrm{ISM}}}\right)}^{1/2}{\left(\displaystyle \frac{{10}^{7}{\rm{K}}}{{T}_{\mathrm{ISM}}}\right)}^{1/2}\,\mathrm{kpc},\end{eqnarray} \tag{ 16 }$

where n_ISM and T_ISM are the number density and temperature of the ISM, respectively. Taking ${n}_{\mathrm{ISM}}={10}^{-3}\ {\mathrm{cm}}^{-3}$ and ${T}_{\mathrm{ISM}}={10}^{7}\ {\rm{K}}$ results in a distance of about 65 kpc, consistent with observations of quasar-driven winds within a factor of a few (e.g., Liu et al. 2013).

We emphasize that Equation (16) should be regarded as the upper limit of the terminal distance. When winds propagate outward, they will be contaminated by the gas in the ISM, leading to a decrease in radial velocity as a result of conservation of momentum, and the terminal distance will become smaller. The magnitude of the decrease depends on the contrast between the momentum flux of the wind and the mass of the ISM gas that is picked up in the winds.

This work is supported in part by the National Key Research and Development Program of China (grant No. 2016YFA0400704), the Natural Science Foundation of China (grants 11573051, 11633006, 11650110427, 11661161012), the Key Research Program of Frontier Sciences of CAS (No. QYZDJSSW-SYS008), and the Astronomical Big Data Joint Research Center cofounded by the National Astronomical Observatories, Chinese Academy of Sciences and the Alibaba Cloud. This work made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory.