NUMERICAL SIMULATION OF HOT ACCRETION FLOWS. III. REVISITING WIND PROPERTIES USING THE TRAJECTORY APPROACH

We have considered two different simulations to perform the "trajectory" analysis. One is the 2D MHD simulation described in Yuan et al. (2012a), and the other is the 3D GRMHD simulation described in Narayan et al. (2012). In the former, the initial condition is a rotating torus with constant specific angular momentum. The density maximum is located at $r=100\;{{r}_{{\rm s}}}$, where ${{r}_{{\rm s}}}\equiv 2\;{{r}_{g}}\equiv 2\;GM\;{{c}^{2}}$. The initial magnetic field is poloidal, a single set of loops confined to the interior of the torus, and the loops are parallel to the density loops. The simulation is performed using the ZEUS code. The readers are referred to Yuan et al. (2012a) for details. In the second simulation, the initial condition is again a rotating torus, but the details of the torus are different. It has inner and outer edges at ${{r}_{{\rm in}}}=10\;{{r}_{g}}$ and ${{r}_{{\rm out}}}=1000\;{{r}_{g}}$, respectively. The simulation domain ranges from close to the black hole to $\sim {{10}^{5}}\;{{r}_{g}}$. The initial magnetic field is again purely poloidal. Different from the first simulation, the magnetic field is broken into eight poloidal loops of alternating polarity. Each loop carries the same amount of magnetic flux, so that the black hole is unable to acquire a large net flux over the course of the simulation. The simulation is performed using the HARM code. It is run for a time of $2\times {{10}^{5}}\;GM/{{c}^{3}}$ and achieves inflow equilibrium (i.e., accretion has reached a steady state) out to a radius $\sim 90\;{{r}_{g}}$. The readers are referred to Narayan et al. (2012) for details. By using the trajectory approach to study the wind properties using data from both simulations, we find that the main results are very similar. Therefore, in this paper, we choose to only focus on the latter simulation (3D GRMHD). Throughout this paper, we use the spherical coordinate. The physical quantities in the present paper, if not specified, are in units of $G=M=c=1$.

Trajectory is related to the Lagrangian description of fluid, obtained by following the motion of fluid elements at consecutive times. Streamline is associated with the Euler description of fluid, obtained by connecting the velocity vectors of adjacent fluid elements at a given time. Trajectory is only equivalent to the streamline for strictly steady motion, which is not the case for accretion flow since it is always turbulent. Streamlines are easy to obtain and are widely used in the literature. On the other hand, obtaining Lagrangian trajectories is much more time-consuming than streamlines, but they loyally reflect the motion of fluid elements. For our purpose, we should consider trajectories rather than streamlines.

To get the trajectory, we first choose a set of "test particles" in the simulation domain within the outermost radius where inflow equilibrium is achieved, $\sim 90\;{{r}_{g}}$, at a given snapshot at time t of the simulation. They are not real particles, but a collection of spatial coordinates as the starting point for the trajectory calculation. With their velocities interpolated from simulation data, we can then obtain their location at time $t+\delta t.$ This process is then continued until the end of the simulation or when test particles leave the simulation domain. We use the "VISIT" software to perform particle trajectory calculation, which can do interpolation with a controlled precision. Obviously, to obtain robust particle trajectories, the time step of the simulation data output $\delta t$ must be sufficiently small. This timescale depends on how fast the velocity of the particles changes during their motion, which is the Keplerian timescale. So the Keplerian timescale must be properly time-resolved. We have tested in our calculations using different time resolutions and compared the corresponding trajectory to see whether the results converge. If not, we shorten the time step of the simulation data output. The time step we actually use in obtaining the trajectory is roughly the Keplerian timescale at $r\approx 6\;{{r}_{g}}$, which is much shorter than the Keplerian timescale of most of the region of the accretion flow. Taking 100 particles at various radii as an example, we have tested two time steps, with one being the Keplerian timescale at 6 ${{r}_{g}}$ and another the Keplerian timescale at 4 ${{r}_{g}}$. We found that the results such as the particle trajectories are largely indistinguishable.

Figure 1 shows the trajectories of sample test particles that originated from locations $(r,\phi )=\left( 80\;{{r}_{g}},0 \right)$ and $\theta =0-\pi$. The left and right panels show the trajectories at the 3D and 2D ($r-\theta$ plane) space, respectively. The results are similar for test particles originating from different radii.

From the figure we clearly see that winds are evident. They are largely located in the polar region, i.e., $\theta \sim 0-\pi /4$ and $\theta =3/4\pi -\pi$, and are symmetric with respect to the equatorial plane. The trajectories of many test particles in this region are almost straight lines in the poloidal plane and $\theta \;\sim$ const., indicating that turbulence is weak. There are also many other more turbulent trajectories with varying θ values during outflowing motion. The specific angular momentum of wind particles is found to be larger than that of the inflow, consistent with Yuan et al. (2012a). The wind region almost overlaps with the usual "coronal" region above the main body of accretion flow (refer to Figure 4 in Yuan & Narayan 2014), and the boundary of the wind region is about the surface of the accretion flow defined by the density scale height. In other words, the disk corona is outflowing. The main disk body, i.e., $\theta =\pi /4-3/4\pi$, is the "inflow" region. Most of the test particles originating from this region move inward, and their motion is turbulent. This turbulent motion is due to the magnetorotational instability (MRI; Balbus & Hawley 1998), as usual. Interestingly, we also find that some test particles originating from this region first move vertically toward the coronal region and then escape outward radially as disk wind. Such vertical motion is present in almost any radius and is perhaps an indicator of the magnetic buoyancy. This supplies new gas from the disk body to the corona/wind. Overall, the structure of the accretion flow is that the inflow region corresponds to the main disk body, while the wind region corresponds to the disk corona. This picture is consistent with Sadowski et al. (2013; see their Figure 16).

Zoom In Zoom Out Reset image size

To calculate the mass flux of the inflow and outflow, and to analyze the wind properties, we distinguish between various types of particle trajectories. In Figure 2, we show characteristic types of characteristic particle trajectories originating from radius r. Distinguishing them is crucial for calculating the mass fluxes of the inflow and outflow correctly, as well as analyzing the wind properties.

• 1.  
  The red and black lines represent outflows and inflows, respectively.
• 2.  
  The red solid line represents a real outflow, where the particle keeps moving outward and never crosses the radius r during its motion.
• 3.  
  The red dashed and red dotted lines represent turbulent outflows, where the particles first move outward but will later return and cross the radius r during their motion. Although the particles eventually move outward and inward, they both belong to the "turbulent outflow" category when we calculate the mass flux.
• 4.  
  It is similar for the three black lines except that they all represent inflow. The solid line represents real inflow, while both the dashed and dotted lines represent turbulent inflow.

To calculate the mass flux (or mass flow rate) of the wind ${{\dot{M}}_{{\rm wind}}}(r)$ at a given time, we first choose test particles initially distributed at fixed radius r with different θ and ϕ and obtain their trajectories. The "real outflow rate" is then calculated by summing up the corresponding mass flux carried by test particles whose trajectories belong to the real outflow (the red solid line in Figure 2)

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm wind}}}(r)=\mathop{\sum }\limits_{i}{{\rho }_{i}}(r){{v}_{r,i}}(r){{r}^{2}}{\rm sin} (\theta )\delta {{\theta }_{i}}\delta {{\phi }_{i}}.\end{eqnarray} \tag{ 4 }$$

Here ${{\rho }_{i}}(r)$ and ${{v}_{r,i}}$ are the mass density and radial velocity at the location where test particle "i" originates, and $\delta {{\theta }_{i}}$ and $\delta {{\phi }_{i}}$ are the ranges of θ and ϕ the particles occupy. We can obtain the "real inflow rate" with a similar approach. In this case, we sum up the corresponding mass flux of test particles whose trajectories are analogous to the black solid line in Figure 2. The mass flux corresponding to the dotted and dashed red lines is the "turbulent outflow rate," while that corresponding to the dashed and dotted black lines is the "turbulent inflow rate."

Zoom In Zoom Out Reset image size

Figure 3 shows the mass flow rate at $r=80\;{{r}_{g}}$ per unit θ and ϕ as a function of θ. The left panel shows the mass flux at two different times but the same ϕ, while the right panel shows the mass flux at two different ϕ but the same time. A positive value is for real outflow, while a negative value is for total inflow, i.e., including both the "real inflow" and "turbulent inflow." By comparing the black and red lines in the left and right panels, we see that the specific values of θ at which real outflows reside change with ϕ and time. In other words, for a fixed θ, the flow can be inflow or outflow for different time t and ϕ.⁷ This result indicates that if we move the $t\phi$ average inside the integrals of Equation (2), i.e., to integrate max(${{\langle \rho {{v}_{r}}\rangle }_{t\phi }}$, 0) as in Narayan et al. (2012), significant mass flux of the real outflow will be canceled and we will substantially underestimate the mass flux of the real outflow. This is the main reason why Narayan et al. (2012) reported much weaker real outflow than Yuan et al. (2012a). This also explains why Narayan et al. (2012) reported that the mass flux of the real outflow does not converge with time since more mass flux will be canceled if the period of integration time is longer. Figure 4 shows the mass flow rate integrated over all ϕ and averaged from time t = 100,000 to 120,000 M at radius $r=80\;{{r}_{g}}$ as a function of θ. As with Figure 3, positive and negative values are for real outflow and total inflow, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 3 and 4 we see that, consistent with the qualitative result shown in Figure 1, the inflow primarily takes place in the main disk body, while most of the mass flux of the real outflow occurs in the range of $\theta$ = 30°–60°. But as we have pointed out in Section 3.1, for some values of ϕ at any radius, some real outflow also exists in the main disk body. This can be seen from Figure 4, which shows some significant mass flux of the real outflow even close to the equatorial plane, i.e., $\theta \sim 100{}^\circ$. We hardly see this wind in Figure 1 because these wind particles first quickly move vertically to the surface of the accretion flow before moving outward as outflowing wind.

The next important questions are, what is the mass flux of the real outflow, and how does it compare with the total outflow rate calculated by Equation (2)? By integrating the wind mass flux shown in Figure 4 over θ but without doing the time average, we can obtain the total real outflow rate at a given time and radius. Then we can obtain the radial profile of mass flux of the real outflow. The red dashed line in Figure 5 shows the result at $t=100,000\;M$.⁸ We find that the radial profile of the mass flux of wind can be well described by

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm wind}}}(r)\approx {{\dot{M}}_{{\rm BH}}}{{\left( \frac{r}{40\;{{r}_{g}}} \right)}^{s}}\equiv {{\dot{M}}_{{\rm BH}}}{{\left( \frac{r}{20\;{{r}_{{\rm s}}}} \right)}^{s}},\;s\approx 1,\end{eqnarray} \tag{ 5 }$$

where ${{\dot{M}}_{{\rm BH}}}$ is the mass accretion rate at the black hole horizon. So the mass flux of the wind at 40 ${{r}_{g}}$ is equal to the mass accretion rate to the black hole ${{\dot{M}}_{{\rm BH}}}$. Such a power-law distribution is likely valid up to the outer boundary of the accretion flow and then quickly decreases beyond the outer boundary (F. Yuan et al. 2015, in preparation). For comparison, we have also shown in the figure by the red and black solid lines the total outflow and total inflow rates calculated by Equations (1) and (2), and the mass flux of the real outflow calculated following the method in Narayan et al. (2012) by the red dotted line. We can see that the mass flux calculated by Equation (2) is equal to ${{\dot{M}}_{{\rm BH}}}$ at $30\;{{r}_{g}}$, while that calculated by the Narayan et al. (2012) method is much weaker, equal to ${{\dot{M}}_{{\rm BH}}}$ at $100\;{{r}_{g}}$. Sadowski et al. (2013) adopted the same method as Narayan et al. (2012) to study the mass flux of the real outflow but could not obtain a radial profile, because they found that the mass flux of the real outflow is too weak and thus difficult to fit using any formula. We also note that the power-law index in Equation (5) is in good agreement with that obtained from the analytical study by Begelman (2012).

Zoom In Zoom Out Reset image size

Also shown in the figure by the red dot-dashed line is the mass flux of the disk jet. We find that when $r\lesssim 400\;{{r}_{g}}$, the mass flux of the disk jet can be described by

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm jet}}}(r)=\frac{1}{20}{{\dot{M}}_{{\rm wind}}}(r).\end{eqnarray} \tag{ 6 }$$

The mass flux increases with radius, indicating that the jet is gradually supplied by matter from the wind. When $r\gtrsim 400\;{{r}_{g}}$, however, the mass flux of the jet is found to be almost saturated.

In this subsection, we consider the evolution of poloidal velocity in the wind. This is the dominant component of the wind velocity at large disk radii once magnetic field becomes subdominant. Figure 6 shows the poloidal speed of the real outflow as a function of θ for various ϕ at three different radii, $r=40,80$, and $160\;{{r}_{g}}$, and t = 100,000 M. To obtain these three plots, we first choose some test particles at various θ and ϕ at the three radii and obtain their trajectories. We then select those particles corresponding to real outflow and obtain their poloidal speed. We can see that the poloidal speed as a function of θ has a sharp jump at $\theta \sim 15{}^\circ$ away from the rotation axis. The poloidal speed of outflow close to the axis is $\gtrsim 0.3\;c$, much larger than that away from the axis, which is $\lesssim 0.05\;c$. We thus can naturally identify the real outflow within $\theta \sim 15{}^\circ$ to the axis as the "disk jet," while the real outflow out of this range is wind. Note that the simulation data we use are for a non-rotating black hole, so the presence of the disk jet is irrelevant to black hole spin, although the spin of the black hole may strengthen the disk jet (e.g., Sadowski et al. 2013). The disk jet originates from the inner disk region and is powered by the rotation energy of the accretion flow. This is different from the Blandford–Znajek jet originating from the black hole horizon (Blandford & Znajek 1977), which is powered by the black hole spin energy. Other differences between the two types of jet include that the disk jet is sub-relativistic and matter dominated, while the Blandford–Znajek jet is relativistic and Poynting flux dominated (see Yuan & Narayan 2014 for a summary).

Zoom In Zoom Out Reset image size

By comparing the three plots in Figure 6, we see that the poloidal wind velocity seems to decrease with increasing radius. To quantify, we calculate the mass-flux-weighted poloidal speed of the real outflow as a function of radius. We distinguish the two types of the real outflow, i.e., wind and disk jet, in our calculations. For this purpose, when calculating the poloidal wind velocity, the integration is only over the range of $15{}^\circ \lesssim \theta \lesssim 165{}^\circ$, while for the entire outflow, we integrate over all θ. The results are shown by blue (total outflow) and red (wind only) dots in Figure 7. The poloidal velocity of the total outflow is slightly higher than that of wind, as expected. The lines are fitting functions, approximately given by

$$\begin{eqnarray}&&{{v}_{{\rm p},{\rm wind}}}(r)\approx 0.21{{v}_{k}}(r).\end{eqnarray} \tag{ 7 }$$

Here ${{v}_{k}}(r)\left( \equiv {{({\rm GM}/r)}^{1/2}} \right)$ is the Keplerian speed at radius r. We note that this equation describes the poloidal velocity of wind at radius r. On the one hand, the wind can be launched from any radius $\lesssim r$. On the other hand, we see in Figure 5 that the wind mass flux increases rapidly with radius. Therefore, the mass-flux-weighted wind velocity (Equation (7)) primarily reflects wind launched close to radius r.

Zoom In Zoom Out Reset image size

We have also calculated the evolution of the poloidal wind velocity along individual test particle trajectories. The results are shown in Figure 8. In the left panel, six representative test particles are shown, initially located at $r=80\;{{r}_{g}}$ but different θ. Specifically, the red dot–dashed line corresponds to $\theta =10{}^\circ$, i.e, within the disk jet region, while others have $\theta \gtrsim 20{}^\circ$. The three red lines have a positive Bernoulli parameter at 80 ${{r}_{g}}$, while the three blue lines have $Be\lt 0$ at 80 ${{r}_{g}}$ (refer to Figure 9). The acceleration of the red dot-dashed line is the most significant, indicating strong acceleration in the jet region. This is confirmed by our detailed analysis of the acceleration mechanism in Section 4. We see that the blue dashed and red solid lines in the left panel of the figure also show strong acceleration. This is likely because these two test particles later enter the "jet region" although they are initially located out of this region. In fact, we find that the θ values of many test particles change significantly as they travel outward. For other lines (particles always in the wind region), while there are fluctuations, the poloidal wind velocity roughly remains constant along the particle trajectories, extending from $r=80\;{{r}_{g}}$ to $r\sim 800\;{{r}_{g}}$, regardless of the sign of their initial Be. Beyond $r\sim 800\;{{r}_{g}}$, the poloidal velocity seems to decrease with radius. This is related to the value of Be. We will argue in the next subsection that such a decrease is likely not reliable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the right panel of Figure 8, we show the averaged poloidal velocity of eight test particles along their trajectories. The initial locations of these particles are uniformly distributed in θ. The red and black lines correspond to particles originating from $r=40$ r_g and 160 ${{r}_{g}}$, respectively. The initial poloidal velocities of the red lines are larger than those of the black lines, consistent with Figure 7. We further distinguish the dashed lines, which correspond to particles that are initially located in the wind region ($15{}^\circ \lesssim \theta \lesssim 165{}^\circ$), and the solid lines, corresponding to particles that always stay in the wind region along their trajectories. We see that the particle poloidal velocity either remains constant or slightly increases outward along the trajectory. Solid and dashed lines differ significantly because we are not using the mass-flux-weighted average. The increasing or constant behavior of the poloidal velocity with radius strongly suggests that additional acceleration forces must operate to compensate for the change in gravitational energy, which will be discussed in Section 4.

From Figure 8, we deduce that the asymptotic terminal poloidal wind velocity originating from radius r can be approximated by

$$\begin{eqnarray}&&{{v}_{{\rm p},{\rm term}}}(r)\approx (0.2\sim 0.4){{v}_{k}}(r).\end{eqnarray} \tag{ 8 }$$

Note the different meaning between this equation and Equation (7). This result also compliments the discussion following Equation (7): the measured wind at any given radius, say, r, is a mixture of wind launched from smaller disk radii and wind that is produced more locally. The former typically has higher velocity but carries smaller mass flux, while the latter carries higher mass flux with smaller velocity. Overall, Equations (7) and (8) are approximately consistent with each other, and they are also consistent with values estimated in Yuan et al. (2012a; Section 3.5).

We now discuss the Bernoulli parameter of the wind. Following Penna et al. (2013a), we adopt in the present paper the following definition of the Bernoulli parameter with the magnetic term included:

$$\begin{eqnarray}&&Be=-\frac{\rho {{u}_{t}}+{\Gamma }u{{u}_{t}}+{{b}^{2}}{{u}_{t}}}{\rho }-1,\end{eqnarray} \tag{ 9 }$$

where ${\Gamma }=5/3$ is the adiabatic index, u_t is the time component of the four-velocity, u is the internal energy, and b is magnetic field strength in the fluid frame. The rest mass energy is subtracted. Far from the black hole, the above equation reduces to the Newtonian quantity, i.e., the sum of kinetic energy, gas enthalpy, and magnetic enthalpy.

The evolution of the Bernoulli parameter along the trajectories of representative wind particles is shown in Figure 9. The three red lines have a positive Be at the starting point, while the three blue lines have a negative initial value of Be. From this figure, we find the following results.

First, the value of Be is not a constant along particle trajectories. This is not unexpected because conservation of Be holds only when the accretion flow is strictly steady and inviscid, while real accretion flow is always turbulent. As a result, it is inappropriate to use the sign of Be to judge whether the flow can escape to infinity, especially when the outflow is still within the radius at which turbulence is well developed.

Second, for outflow originating from smaller θ, Be increases to very large values at large radii. On the other hand, the value of Be varies much less significantly and roughly remains constant for the wind originating from closer to the disk surface. Comparing with Figure 8, we see that changes in Be correlate positively with changes in poloidal velocity. For the three trajectories whose Be increases outward, there is also significant acceleration in their poloidal velocities. This is because their θ values gradually decrease and enter into the jet region, thus experiencing strong acceleration.

It is interesting to note that for wind originating within $40{}^\circ \lesssim \theta \lesssim 50{}^\circ$, their value of Be becomes almost constant when $r\gtrsim 800\;{{r}_{g}}$. This corresponds to the slight decrease of the poloidal velocity beyond $\sim 800\;{{r}_{g}}$ shown in the left panel of Figure 8. The reason why Be does not change beyond $800\;{{r}_{g}}$ is because in this region turbulence has not well developed within our simulation time. Note that this radius is different from the inflow equilibrium radius, which is $\sim 90\;{{r}_{g}}$. Within $\sim 90\;{{r}_{g}}$, everything, especially the radial density profile, is fully reliable. Beyond this radius, the density profile is not reliable, but other properties, such as the level of turbulence and subsequently outflow properties, are still reliable up to a much larger "turbulence radius," the limiting radius of turbulence steady state. This radius can be estimated as follows. Turbulence in accretion flow is because of MRI. The fastest growth rate of MRI at radius r is $\sim {\Omega }(r)$. More precisely, it takes 3–4 orbits for MRI to develop and ∼10 orbits to saturate (Hawley et al. 1995). For our simulation time of ${{t}_{{\rm simulation}}}\sim 2\times {{10}^{5}}$, taking a timescale of 3 orbits, we can obtain that the "turbulence radius" is $\sim 500\;{{r}_{g}}$. This is close to the value of 800 r_g mentioned above. Another way to understand the "turbulence radius" is as follows. For a geometrically thick disk, the largest turbulence eddies have size of order r. The corresponding eddy turnover time is $r/\sigma (r)$, where $\sigma (r)$ is the rms turbulent velocity. Our simulation data show $\sigma (r)\sim 0.15{{v}_{k}}(r)$. If at a certain radius the eddy turnover time is substantially smaller than the duration of the simulation, then the local turbulence is likely to have reached quasi-steady state. Therefore, the "turbulence radius" should be some fraction of $\sigma (r){{t}_{{\rm simulation}}}$, which gives a similar result to the above estimation. Therefore, we think that the results beyond $\sim 800\;{{r}_{g}}$ are not reliable. It is very likely that Be will keep changing and the poloidal velocity still remains constant beyond $800\;{{r}_{g}}$. This implies that wind can at least escape beyond the outer boundary of accretion flows. Simulations with longer run times can check this point.

Although Be is not a constant along particle trajectories, it is still useful to decompose Be into individual physical terms and compare their contributions, as shown in Figure 10. The left and right panels correspond to particles originating from jet and wind regions, respectively. In the case of the disk jet, it is mainly the enthalpy that compensates for the increase of gravitational energy and kinetic energy. Magnetic energy also plays an active role at smaller radius ($r\;\lesssim \;100\;{{r}_{g}}$). In the case of wind, the role of magnetic energy appears unimportant. The increase of gravitational energy is mainly compensated by the reduction of specific enthalpy and kinetic energy. But we note that although the total kinetic energy decreases along the trajectory, the poloidal component does not. It usually keeps constant, as shown by Figure 8. This corresponds to the work done by the centrifugal force, as we will discuss in Section 4.

Zoom In Zoom Out Reset image size

Based on the trajectory analysis, we now calculate the energy and momentum fluxes from both the wind and jet as follows:

$$\begin{eqnarray}&&{{\dot{E}}_{{\rm jet}({\rm wind})}}(r)=\int \frac{1}{2}\rho (r,\theta ,\phi )v_{p}^{3}(r,\theta ,\phi ){{r}^{2}}{\rm sin} (\theta )d\theta d\phi ,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{\dot{P}}_{{\rm jet}({\rm wind})}}(r)=\int \rho (r,\theta ,\phi )v_{p}^{2}(r,\theta ,\phi ){{r}^{2}}{\rm sin} (\theta )d\theta d\phi .\end{eqnarray} \tag{ 11 }$$

Here v_p(r) is the poloidal velocity at radius r, for which we assume ${{v}_{p}}(r)\approx {{v}_{r}}(r)$ in our calculation because we find ${{v}_{\theta }}\ll {{v}_{r}}$. The integration over θ for wind and the disk jet is bounded by $\theta \approx 15{}^\circ$, according to Figure 6.

Figures 11 and 12 show the radial profiles of energy and momentum fluxes, respectively, calculated at T = 100,000 M. In both figures, the red solid and blue dashed lines are for the wind and disk jet, respectively. We see that the energy flux of the wind is $\gtrsim 3$ times stronger than that of the disk jet, while the contrast in the momentum flux between wind and jet is much larger. This is again mainly because of the low density in the disk jet. From Figure 11, we see that the energy flux of the wind rapidly increases at small radii and then becomes almost saturated at $r\gtrsim 40\;{{r}_{g}}$. The rapid increase corresponds to the rapid increase of the wind mass flux with radius when $r\lesssim 40\;{{r}_{g}}$ (see Figure 5). For $r\gtrsim 40\;{{r}_{g}}$, we have

$$\begin{eqnarray}&&{{\dot{E}}_{{\rm wind}}}(r)\approx \frac{1}{2}{{\dot{M}}_{{\rm wind}}}(r)v_{{\rm p},{\rm wind}}^{2}(r)\approx \frac{1}{1000}{{\dot{M}}_{{\rm BH}}}{{c}^{2}}.\end{eqnarray} \tag{ 12 }$$

In the above calculation, Equations (5) and (7) are used. This result indicates that the energy flux at large radius is roughly saturated, consistent with Figure 11. The main reason energy flux saturates is s = 1 in Equation (5). For the momentum flux of the wind, we have

$$\begin{eqnarray}&&{{\dot{P}}_{{\rm wind}}}(r)\approx {{\dot{M}}_{{\rm wind}}}(r){{v}_{{\rm p},{\rm wind}}}(r)\propto {{r}^{1/2}}.\end{eqnarray} \tag{ 13 }$$

This is consistent with the result shown in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The energy flux obtained in Equation (12) is in good agreement with that required in large-scale AGN feedback simulations (e.g., Ciotti et al. 2010; Gaspari et al. 2012). In these works, AGN feedback is involved to heat the intercluster medium to compensate for the rapid cooling rate in the systems (i.e., the cooling flow problem). It was found that to be consistent with observations of both isolated galaxies and galaxy clusters, the required "mechanical feedback efficiency," defined as $\epsilon \equiv {{\dot{E}}_{{\rm wind}}}/{{\dot{M}}_{{\rm BH}}}{{c}^{2}}$, must be in the range of $\sim {{10}^{-3}}-{{10}^{-4}}$. Our results provide a natural explanation for the required value of , at least when the AGN is in the hot accretion mode.

Our results highlight the importance of wind over jet on AGN feedback. On the other hand, one important caveat is that our calculation is based on the simulation of a Schwarzschild black hole. If the black hole is rapidly spinning, the power of the disk jet is expected to become stronger (Sadowski et al. 2013). In addition to the disk jet, a Poynting-flux-dominated jet (BZ jet) will be produced through the black hole horizon, powered by the black hole spin (Blandford & Znajek 1977; Komissarov 2001; McKinney 2005; Hawley & Krolik 2006; Tchekhovskoy et al. 2010, 2011; Penna et al. 2013b). The dependence of the power of the BZ jet on spin is clear (see references above), while the dependence of the disk jet and wind on spin remains to be investigated, perhaps less sensitive compared to the BZ jet.

In addition to the black hole spin, another parameter is the magnetic flux threading the inner region of the accretion flow. In our simulation, this flux is small and the accretion flow is called to be in the "SANE" state (Narayan et al. 2012). If the flux is large, the system enters the "magnetically arrested disk" (MAD) state. There have been some studies on the dependence of the jet and wind power on the magnetic flux. For example, it has been found that in the MAD state, the power of the BZ jet will dominate the disk jet (Narayan et al. 2003; Tchekhovskoy et al. 2010; Penna et al. 2013a; see review by Yuan & Narayan 2014). Sadowski et al. (2013) compared the power of the wind and the jet and found that, if the black hole spin and magnetic flux at the horizon are large, jet power will usually dominate the wind power. But note that as we have described in the previous section, their estimation of the mass flux of wind should be regarded as a lower limit. For a rapidly spinning black hole accreting in the MAD limit, the power of the jet is even larger than the accretion power (Tchekhovskoy et al. 2011, 2012). A systematic study on the dependence of the jet and wind power on magnetic flux and black hole spin is required. We plan to revisit this problem in the next work.

To study the production and acceleration mechanism of the wind and jet, we have calculated the forces at the jet and wind region at a single point and a given snapshot. The locations where we evaluate the forces correspond to real outflow based on our particle trajectory study. We note that forces are stochastic, but we found that our analysis can be regarded as representative, except in the cases we will mention below. Since we evaluate the forces in the frame co-rotating with the flow at the "evaluation location," we should also include the centrifugal force, in addition to the gravitational force, gas pressure gradient, and Lorentz force.⁹ We show the results for three representative points in the disk jet, wind, and main disk regions in Figure 13. For the wind, the main driving forces are the centrifugal force and magnetic pressure gradient. From the figure we notice that the gradient of the magnetic pressure is "downward," pointing toward the positive θ direction. This is somewhat surprising because we usually expect that magnetic field becomes weaker away from the main disk body toward the coronal region. This reflects the strong fluctuation of the accretion flow. In fact, if we choose another time or another location to do the force analysis, we very likely find that the gradient of the magnetic pressure becomes "upward." For the same reason, the direction of the gas pressure gradient also strongly fluctuates with time and location. But statistically, the gradients of both the gas and magnetic pressure are pointing along the positive r direction and thus are helpful to the acceleration of wind. Their magnitudes are also comparable to the centrifugal force, as shown by Figure 13.

Zoom In Zoom Out Reset image size

From the figure we can see that the magnitude of the centrifugal force is larger than (or at least comparable to in general) the gravitational force. This means that the specific angular momentum of the wind is larger than or close to the Keplerian value. For comparison, we also show in the figure the force analysis in the inflow region. We see that the magnitude of the centrifugal force is now smaller than the gravitational force.¹⁰ This is consistent with the result obtained in Yuan et al. (2012a) that the specific angular momentum of the outflow is systematically larger than that of the inflow. This implies that some angular momentum is transferred from some fluid element to another, likely by the magnetic field lines from the main disk body to the coronal region. Once the combination of the centrifugal force and the pressure gradient exceeds the gravitational force, wind will be accelerated.

According to the above analysis, the mechanism of the acceleration of wind is similar to the Blandford & Payne (1982) mechanism in the sense that the centrifugal force plays an important role. The differences are that the gradient of the pressure plays a comparable role compared with the centrifugal force. In addition, there is no large-scale magnetic field formed in our simulation, and the wind region is not force-free.

In the case of the disk jet, the acceleration mechanism is completely different from that of Blandford & Payne (1982). Here the dominant force is the gradient of the toroidal magnetic pressure, consistent with the result of Hawley & Krolik (2006). This mechanism is the so-called magnetic tower mechanism (Lynden-Bell 2003; see also Shibata & Uchida 1985; Kato et al. 2004). The reason for the much higher jet velocity compared with the wind velocity is related to the strong magnetic field and low density. In the case of both the wind and jet, the energy of the outflow mainly comes from the rotation energy of the accretion flow. The rotation energy is converted into the magnetic energy, which then is converted into the kinetic energy of the wind and jet.