Numerical Simulation of Hot Accretion Flows (III): Revisiting wind properties using trajectory approach
Abstract
Previous MHD simulations have shown that wind (i.e., uncollimated outflow) must exist in black hole hot accretion flows. In this paper, we continue our study by investigating the detailed properties of wind, such as mass flux and poloidal speed, and the mechanism of wind production. For this aim, we make use of a three dimensional GRMHD simulation of hot accretion flows around a Schwarzschild black hole. The simulation is designed so that the magnetic flux is not accumulated significantly around the black hole. To distinguish real wind from turbulent outflows, we track the trajectories of the virtual Largrangian particles from simulation data. We find two types of real outflows, i.e., a quasirelativistic jet close to the axis and a subrelativistic wind subtending a much larger solid angle. We confirm that the mass flux of wind is very significant and most of the wind originates from the surface layer of the accretion flow. The radial profile of the wind mass flux can be described by , with being the mass accretion rate at the black hole horizon and being the Schwarzschild radius. The poloidal wind speed almost remains constant once they are produced, but the fluxweighted wind speed roughly follows , with being the Keplerian speed at radius . The mass flux of jet is much lower but the speed is much higher, . Consequently, both the energy and momentum fluxes of the wind are much larger than those of the jet. We find that the wind is produced and accelerated primarily by the combination of centrifugal force and magnetic pressure gradient, while the jet is mainly accelerated by magnetic pressure gradient. Finally, we find that the wind production efficiency , in good agreement with the value required from largescale galaxy simulations with AGN feedback.
Subject headings:
galaxies: AGN  BH: feedback  BH physics: reradiation1. Introduction
Black hole accretion models can be divided into two classes based on the temperature of the accretion flow. One is cold accretion models such as the standard thin disk (Shakura & Sunayev 1976; Pringle 1981), the other is the hot accretion flow such as advectiondominated accretion flow (Narayan & Yi 1994, 1995; see Yuan & Narayan 2014 for the recent review of the theory of hot accretion flow and its applications). Hot accretion flows are believed to exist in lowluminosity active galactic nuclei, which likely reside in majority of galaxies, and hard/quiescent states of black hole Xray binaries. One important question in the study of hot accretion flows in the recent years is related to winds, i.e., uncollimated mass outflow. On the one hand, whether winds are present or not is a fundamental question in the dynamics of accretion, and the presence of wind help explain many observations, including the spectrum of black hole sources (e.g., Yuan, Quataert & Narayan 2003), emission lines from accretion flow (e.g., Wang et al. 2013), the Fermi bubbles in the Galactic center (Mou et al. 2014), and perhaps even some direct observations on outflow (e.g., Crenshaw et al. 2003; Crenshaw & Kraemer 2012; Yuan, Bu & Wu 2012). On the other hand, AGN feedback is widely believed to play a crucial role in galaxy formation and evolution (Fabian 2012; Kormendy & Ho 2013), while winds produced by the AGN accretion flow is one of the most important ingredients for such a feedback because they could interact with the ISM in the host galaxy by exchanging momentum and depositing energy (King 2003; Debuhr et al. 2010; Ciotti, Ostriker & Proga 2010; Ostriker et al. 2010; Novak et al. 2011; Gan et al. 2014). The aim of the present work is to understand the detailed wind properties from accretion flow, which will provide important input for studies of AGN feedback.
The study of winds from hot accretion flows can be traced back to Stone, Pringle & Begelman (1999; see also Igumenshchev & Abramowicz 1999, 2000). They performed the first global hydrodynamical numerical simulation of hot accretion flow and calculated the following timeaveraged radial profiles of inflow and outflow rates,
(1) 
(2) 
where the angle brackets represent time averages (and also average over the azimuthal angle in the case of 3D simulations). We emphasize this point because the order of doing timeaverage and the integral will make significant differences, as we will show later in this paper. Note also that the outflow rate calculated by eq. (2) does not necessarily represent the mass flux of “real outflow”, because the positive radial velocity may just come from the turbulent motion of the accretion flow. The most important result they obtained is that the inflow rate based on eq. (1) follows a powerlaw function of radius,
(3) 
Here is the mass inflow rate at the outer boundary . The radial dynamical range of this simulation is rather small, spanning about two orders of magnitude in radius. But the results were later confirmed by simulations with a much larger radial dynamical range of four orders of magnitiude (Yuan, Wu & Bu 2012). Moreover, MHD simulations yield very similar results that typically (e.g., Stone & Pringle 2001; Hawley & Balbus 2002; Igumenshchev et al. 2003; Pen et al. 2003; see review in Yuan, Wu & Bu 2012).
It is exciting to note that the predicted inward decrease of accretion rate has soon been confirmed by two observations, both are on Sgr A*. One is the detection of radio polarization at a level of (e.g., Aitken et al. 2000; Bower et al. 2003; Marrone et al. 2007). Such high polarization requires that the mass accretion rate close to the black hole horizon must be within a certain range, which is two orders of magnitude lower than the Bondi rate obtained from Chandra observations. The other evidence is from the Chandra observation of the iron emission lines originated from the hot accretion flow (Wang et al. 2013). The modeling to the lines indicates a flat radial density profile near the Bondi radius, which confirms that the mass accretion rate decreases with decreasing radius. This is because, if the mass accretion rate were a constant of radius, the density profile would be much steeper.
Two competing models have been proposed to explain the above numerical simulation result. In the adiabatic inflowoutflow solution (ADIOS), the inward decrease of mass accretion rate is due to the mass lost in the wind (Blandford & Begelman 1999; 2004; Begelman 2012). In the early works of Blandford & Begelman (1999; 2004), the value of parameter in eq. 3 is a free parameter. But in the more recent work of Begelman (2012), the value of is argued to be close to unity. The other model is the convectiondominated accretion flow (CDAF) model. In this model, the accretion flow is assumed to be convectively unstable. The inward decrease of accretion rate is explained as more and more gas is locked in convective eddies (Narayan et al. 2000; Quataert & Gruzinov 2000; Abramowicz et al. 2002; Igumenshchev 2002). For a long time, it is unclear which scenario is physical.
Three numerical simulation works have been conducted to investigate this problem (Narayan et al. 2012; Yuan, Bu & Wu 2012; Li, Ostriker & Sunyaev 2013). Both Narayan et al. (2012) and Yuan, Bu & Wu (2012) found that the hot accretion flow is convectively stable. This indicates that the CDAF model may not apply, leaving outflow/wind as the only possible solution. The fundamental question is, how strong the wind is. Narayan et al. (2012) calculated the outflow rate based on eq. (2), except that they move the average inside the integral. Although this calculation underestimates the mass flux of real outflow, as we will show later in this paper, it eliminates contributions from turbulent motion and produces substantially lower outflow rate than eq. 2. In fact, only upper limit was reported since the outflow rate was found to not converge with time. On the other hand, Yuan, Bu & Wu (2012) systematically compared the properties of inflow and outflow, such as angular momentum and temperature, and found that they are quite different. They therefore concluded that systematic outflow must exist and the outflow rate must be a significant fraction of that indicated by eq. (2). They argued that the rather weak outflow rate obtained in Narayan et al. (2012) is because outflow is intrinsically instantaneous. The outflow stream can wander around in 3D space thus will be cancelled if the timeaverage is performed first. Their work indicates that the mass lost via the wind is the reason for the inward decrease of the accretion rate (eq. 3). The hydrodynamical simulations by Li, Ostriker & Sunyaev (2013) obtained a similar conclusion as Yuan, Bu & Wu (2012). Begelman (2012) and Gu (2014) studied why winds should exist in hot accretion flows.
The aim of the present work is two fold. First, Yuan, Bu & Wu (2012) showed the existence of outflow only based on some indirect arguments, so it is necessary to show the existence of wind in a more direct way. Second, we want to quantitatively calculate the properties of wind, including the mass flux, angular distribution, and velocity. These properties are especially important to determining how effective the interaction between wind and ISM is. The amount of mass flux of wind is also useful to resolve the discrepancy on the mass flux of wind between Narayan et al. (2012) and Yuan, Bu & Wu (2012). Sadowski et al. (2013) have also studied the properties of wind and jet under various parameters such as black hole spin. Same with Narayan et al. (2012), their calculations are based on the timeaveraged quantities thus may only give a lower limit.
We use a “trajectory” approach in this work for our abovementioned aims. That is, we use the numerical simulation data of hot accretion flow to follow the trajectories of some “test particles” so as to see whether the particles can really escape or simply have turbulent motions, and further calculate the properties of the wind. The paper is organized as follows. In §2, we will describe the simulation data based on which we perform the analysis, and the “trajectory” method we use. The results are presented in §3. §4 is devoted to the analysis of the acceleration mechanism of wind. We then summarize and conclude in §5.
2. Simulation data and the trajectory method
2.1. Simulation data of hot accretion flow models
We have considered two different simulations to perform the “trajectory” analysis. One is the two dimensional MHD simulation described in Yuan, Bu & Wu (2012), and the other is the three dimensional 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 , where . 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, Bu & Wu (2012) 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 and , respectively. The simulation domains ranges from close to the black hole to . 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 and achieves inflow equilibrium (i.e., accretion has reached a steady state) out to a radius . 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 .
2.2. Trajectory method
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 timeconsuming than streamlines, but they loyally reflects 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 outmost radius where inflow equilibrium is achieved, , at a given snapshot at time 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 . 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 must be sufficiently small. This timescale depends on how fast the velocity of the particles change during their motion, which is the Keplerian timescale. So the Keplerian timescale must be properly timeresolved. 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 , which is much shorter than the Keplerian timescale of most of the region of the accretion flow. Taking 100 particles at various radius as example, we have tested two time steps, with one being the Keplerian timescale at and another the Keplerian timescale at . We found that the results such as the particle trajectories are largely indistinguishable.
3. Results
Before describing our results, we first define some terminologies. We use “outflow” to describe any flow with a positive radial velocity , i.e., flowing outward. This includes both “turbulent outflow” and “real outflow”. The difference between them is that in the former case the test particle will eventually return and join the accretion flow after flowing outward for some distance, while in the latter case the test particle continues to flow outward and eventually escapes the outer boundary of the simulation domain^{1}^{1}1In some cases, e.g., if the initial radius of test particles are small, some particles could not escape beyond the outer boundary by the end of the simulation. On the other hand, they keep moving radially outward without any signs of return (the poloidal speed does not decrease outward, refer to Fig. 8), we regard them as representing real outflow.. “Real outflow” consists of two components, i.e., disk jet and wind. Here, “disk jet” is different from the BlandfordZnajek jet in several aspects, as we have summarized in Yuan & Narayan (2014). Very briefly, “disk jet” originates from the innermost region of the accretion flow. It is quasirelativistic and matterdominated. The BlandfordZnajek jet is powered by the spin of the black hole. It is relativistic and Poynting fluxdominated. In the present paper, we do not have a BlandfordZnajek jet since the simulation data we use is for a Schwarzschild black hole. As we will describe in detail later, we find that the disk jet is confined in a region of away from the polar axis while the wind is located between the jet boundary and the surface of the accretion flow. The velocity of the jet is much higher than that of the wind.
Note that the definition of wind we adopt here is different from that adopted in some literature (e.g., Narayan et al. 2012; Sadowski et al. 2013), where they require that the Bernoulli parameter of wind must satisfy . The Bernoulli parameter of our wind can be of any value, negative or positive, at any radius. Our argument is that for nonsteady accretion flow, is not constant along trajectories, but usually increases outward (refer to Fig. 9). This means that even though at a certain radius, the wind particles can still escape to infinity. In fact, as we will describe later, we find that no matter what value is, the poloidal speed of wind does not decrease when they propagate outward until a radius within which turbulence is developed (refer to Fig. 8 and relevant discussions later). In addition, technically the value of of wind should depend on the initial condition of the simulation because of energy conservation. In many simulations including the current one, the initial condition is a bound torus thus is negative. But in reality, the accretion flow comes from much further away so is more likely positive. This implies that the value of obtained in simulation should be regarded as an lower limit.
3.1. Overall result: confirmation of the wind
Figure 1 shows the trajectories of sample test particles originated from locations and . The left and right panels show the trajectories at the three and two ( plane) dimensional 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., and , 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 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, Bu & Wu (2012). The wind region almost overlaps with the usual “coronal” region above the main body of accretion flow (refer to Fig. 4 in Yuan & Narayan 2014 for the structure of accretion flow), 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., , 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 originated 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 Fig. 16).
3.2. The mass flux of outflow and inflow
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 originated from radius . Distinguishing them is crucial for calculating the mass fluxes of the inflow and outflow correctly, as well as analyzing the wind properties.

The red and black lines represent outflows and inflows, respectively.

The red solid line represents a real outflow, where the particle keeps moving outward and never cross the radius during its motion.

The red dashed and red dotted lines represent turbulent outflows, where the particles first move outward, but will later return and cross the radius 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.

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 at a given time, we first choose test particles initially distributed at fixed radius 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 Fig. 2)
(4) 
Here and are the mass density and radial velocity at the location where test particle “i” originates, and are the ranges of and the particle 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 of the black solid line in Fig. 2. The mass flux corresponding to the dotted and dashed red lines are “turbulent outflow rate”, while that corresponding to the dashed and dotted black lines are are “turbulent inflow rate”.
Fig. 3 shows the mass flow rate at 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. Positive value is for real outflow while 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 and ^{2}^{2}2The values of time in the left panel and the values of in the right panel are chosen so that the change of inflow and outflow for a fixed is significant.. This result indicates that if we move the average inside the integrals of eq. 2, i.e., to integrate max(,0) as in Narayan et al. (2012), significant mass flux of real outflow will be cancelled 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, Bu & Wu (2012). This also explains why Narayan et al. (2012) reported that a the mass flux of the real outflow does not converge with time since more mass flux will be cancelled if the period of integration time is longer. Fig. 4 shows the mass flow rate integrated over all and averaged from time at radius as a function of . Same with Fig. 3, positive and negative values are for real outflow and total inflow, respectively.
From Figs. 3 & 4 we see that consistent with the qualitative result shown in Fig. 1, the inflow primarily takes place in the main disk body while most of mass flux of the real outflow occurs in the range of . But as we have pointed out in §3.1, for some values of at any radius, some real outflow also exists in the main disk body. This can be seen from Fig. 4 which shows that some significant mass flux of real outflow even close to the equatorial plane, i.e., . We hardly see these wind in Fig. 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 real outflow? How is it compared with the total outflow rate calculated by eq. (2)? By integrating the wind mass flux shown in Fig. 4 over but without doing 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 real outflow. The red dashed line in Fig. 5 shows the result at .^{3}^{3}3We want to point out a caveat here. As we have stated in §2.1, the inflow equilibrium is reached only up to in our simulation, but the xaxis of the figure extends to . We find that the radial profile of the mass flux of wind can be well described by
(5) 
where is the mass accretion rate at the black hole horizon. So the mass flux of the wind at is equal to the mass accretion rate to the black hole . Such a powerlaw distribution is likely valid up to the outer boundary of the accretion flow, and then quickly decreases beyond the outer boundary (Yuan et al. 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 eqs. (2) and (1), and the mass flux of 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 eq. (2) is equal to at ; while that calculated by the Narayan et al. (2012) method is much weaker, equal to at . Sadowski et al. (2013) adopted the same method as Narayan et al. (2012) to study the mass flux of real outflow but could not obtain a radial profile, because they found that the mass flux of the real outflow is too weak thus difficult to fit using any formula. We also note that the powerlaw index in eq. (5) is in good agreement with that obtained from the analytical study by Begelman (2012).
Also shown in the figure by the red dotdashed line is the mass flux of disk jet. The mass flux of disk jet can be described by
(6) 
The mass flux increases with radius, indicating that the jet is gradually supplied by matter from the wind.
3.3. The poloidal speed
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. Fig. 6 shows the poloidal speed of the real outflow as a function of for various at three different radii, and , and . To obtain these three plots, we first choose some test particles at vavious 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 away from the rotation axis. The poloidal speed of outflow close to the axis is , much larger than those away from the axis, which is . We thus can naturally identify the real outflow within to the axis as “disk jet”, while the real outflow out of this range to be wind. Note that the simulation data we use is for a nonrotating 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 diskjet originates from the inner disk region and are powered by the rotation energy of the accretion flow. This is different from the BlandfordZnajek jet originated from the black hole horizon (Blandford & Znajek 1976), which is powered by the black hole spin energy. Other differences between the two types of jet include that the disk jet is subrelativistic and matterdominated, while the BlandfordZnajek jet is relativistic and Poyntingflux dominated (see Yuan & Narayan 2014 for a summary).
By comparing the three plots in Fig. 6, we see that the poloidal wind velocity seems to decrease with increasing radius. To quantify, we calculate the mass fluxweighted poloidal speed of 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 ; while for the entire outflow, we integrate over all . The results are shown by blue (total outflow) and red (wind only) dots in Fig. 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
(7) 
Here is the Keplerian speed at radius . We note that this equation describes the poloidal velocity of wind at radius . On the one hand, the wind can be launched from any radius . On the other hand, we see in Figure 5 that the wind mass flux increases rapidly with radius. Therefore, the massflux weighted wind velocity (7) primarily reflects wind launched close to radius .
We have also calculated the evolution of the poloidal wind velocity along individual test particle trajectories. The results are shown in Fig. 8. In the left panel, six representative test particles are shown, initially located at but different . Specifically, the red dotdashed line corresponds to , i.e, within the disk jet region; while others have . The three red lines have a positive Bernoulli parameter at while the three blue lines have at (refer to Fig. 9). The acceleration of the red dotdashed line is the most significant, indicating strong acceleration in the jet region. This is confirmed by our detailed analysis of the acceleration mechanism in §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 remain constant along the particle trajectories, extending from to , regardless of the sign of their initial . Beyond , the poloidal velocity seems to decrease with radius. This is related to the value of . We will argue in the next subsection that such a decrease is likely not reliable.
In the right panel of Fig. 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 originated from and , respectively. The initial poloidal velocities of the red lines are larger than that of the black lines, consistent with Fig. 7. We further distinguish the dashed lines, which correspond to particles that are initially located in the wind region (), 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 keeps constant or slightly increases outward along the trajectory. Solid and dashed lines differ significantly because we are not using mass fluxweighted average. The increasing or constant behavior of the poloidal velocity with radius strongly suggests additional acceleration forces must operate to compensate for the change in gravitational energy, which will be discussed in §4.
From Fig. 8, we deduce that the asymptotic terminal poloidal wind velocity originated from radius can be approximated by:
(8) 
Note the different meaning between this equation and eq. (7). This result also compliments the discussion following eq. (7): the measured wind at any given radius, say , is a mixture of wind launched from smaller disk radii, and the 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, eqs. (7) and (8) are approximately consistent with each other, and they are also consistent with values estimated in Yuan, Bu & Wu (2012; section 3.5).
3.4. The Bernoulli parameter of real outflow
We now discuss the Bernoulli parameter of the wind. Following Penna et al. (2013b), we adopt in the present paper the following definition of the Bernoulli parameter with the magnetic term included,
(9) 
where is the adiabatic index, is the time component of the fourvelocity, is the internal energy, 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 Fig. 9. The three red lines have a positive at the starting point while the three blue lines have a negative initial value of . From this figure, we find the following results.
First, the value of is not a constant along particle trajectories. This is not unexpected because conservation of 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 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 originated from smaller , increases to very large values at large radii. On the other hand, the value of varies much less significantly and roughly remains constant for the wind originated from closer to the disk surface. Comparing with Fig. 8, we see that changes in correlates positively with changes in poloidal velocity. For the three trajectories whose increases outward, there is also significant acceleration in their poloidal velocities. This is because that their values gradually decreases and enters into the jet region, thus experience strong acceleration.
It is interesting to note that for wind originated within , their value of becomes almost constant when . This corresponds to the slight decrease of the poloidal velocity beyond shown in the left panel of Fig. 8. The reason why does not change beyond 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 . Within , 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 is . More precisely, it takes orbits for MRI to develop and orbits to saturate (Hawley et al. 1995). For our simulation time of , taking a timescale of 3 orbits, we can obtain that the “turbulence radius” is . This is close to the value of 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 . The corresponding eddy turnover time is , where is the rms turbulent velocity. Our simulation data shows . 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 quasisteady state. Therefore the “turbulence radius” should be some fraction of , which gives a similar result to the above estimation. Therefore, we think that the results beyond are not reliable. It is very likely that will keep changing and the poloidal velocity still remains constant beyond . This implies that wind can at least escape beyond the outer boundary of accretion flows. Simulations with longer run time can check this point.
Although is not a constant along particle trajectories, it is still useful to decompose into individual physical terms and compare their contributions, as shown in Fig. 10. The left and right panels correspond to particles originated from jet and wind regions, respectively. In the case of 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 (). 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 Fig. 8. This corresponds to the work done by the centrifugal force, as we will discuss in §4.
3.5. The fluxes of energy and momentum of wind and jet
Based on the trajectory analysis, we now calculate the energy and momentum fluxes from both the wind and jet as follows.
(10) 
(11) 
Here is the poloidal velocity at radius , which we assume in our calculation because we find . The integration over for wind and diskjet is bounded by , according to Fig. 6.
Figs. 11 & 12 show the radial profiles of energy and momentum fluxes, respectively, calculated at . 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 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 diskjet. From Fig. 11, we see that the energy flux of the wind rapidly increases at small radii, and then becomes almost saturated at . The rapid increase corresponds to the rapid increase of the wind mass flux with radius when (see Fig. 5). For , we have
(12) 
In the above calculation, eqs. (5) and (7) are used. This result indicates that the energy flux at large radius is roughly saturated, consistent with Fig. 11. The main reason energy flux saturates is in Eq. 5. For the momentum flux of the wind, we have
(13) 
This is consistent with the result shown in Fig. 12.
The energy flux obtained in eq. (12) is in good agreement with that required in large scale AGN feedback simulations (e.g., Ciotti, Ostriker & Proga 2010; Gaspari, Brighenti & Temi 2012). In these works, AGN feedback is involved to heat the intercluster medium to compensate for 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 , must be in the range of . 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 fluxdominated jet (BZ jet) will be produced through the black hole horizon, powered by the black hole spin (Blandford & Znajek 1976; Komissarov 2001; McKinney 2005; Hawley & Krolik 2006; Tchekhovskoy, Narayan & McKinney 2010, 2011; Penna, Narayan & Sadowski 2013a). The dependence of the power of 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, Igumenshchev & Abramowicz 2003; Tchekhovskoy, Narayan & McKinney 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 previous section, their estimation of the mass flux of wind is should be regarded as a low 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, Narayan & McKinney 2011; Tchekhovskoy, McKinney & Narayan 2012). 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.
4. The mechanism of producing outflow
4.1. Mechanism of acceleration of wind
To study the production and acceleration mechanism of 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 corotating frame corotating with the flow at the “evaluation location”, we should also include the centrifugal force, in addition to the gravitational force, gas pressure gradient, and the Lorentz force^{4}^{4}4We only include the gradient of magnetic pressure. Since , the magnetic tension force is much weaker than the gradient of magnetic pressure so we neglect the tension 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”. Because of 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 direction thus are helpful to the acceleration of wind. Their magnitudes are also comparable to the centrifugal force, as shown by Fig. 13.
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^{5}^{5}5This confirms that the accretion flow is subKeplerian.. This is consistent with the result obtained in Yuan, Bu & Wu (2012) that the specific angular momentum of outflow is systematically larger than that of the inflow. This implies that some angular momentum is transferred from some fluid element to other, 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 Blandofrd & 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 largescale magnetic field formed in our simulation and the wind region is not force free.
4.2. Mechanism of the acceleration of jet
In the case of disk jet, the acceleration mechanism is completely different from the 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 socalled magnetic tower mechanism (LyndenBell 2003; see also Shibata & Uchida 1985; Kato, Mineshige & Shibata 2004). The reason for much higher jet velocity compared with the wind velocity is related with the strong magnetic field and low density. In both the case of 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.
5. Summary
Previous numerical simulations to the black hole accretion flow have shown that the mass accretion rate decreases inward (eq. (3)). It has been proposed that convection or outflow may result in this result. Stability analysis, however, has indicated that the accretion flow is convectively stable when the magnetic field is present. This excludes the possibility of convective instability and leave only the mass loss via outflow as the likely reason. However, previous theoretical works have obtained quite different result on the strength of outflow. In terms of the mass accretion rate onto the black hole horizon , Narayan et al. (2012) found that even at , the mass flux of real outflow is still very weak. On the other hand, Yuan, Bu & Wu (2012) argue that mass flux of outflow should be strong, i.e., a significant fraction of that described by eq. (2).
One of the main aims of the present work therefore is to investigate how strong the real outflow is. It is well known that accretion flow is turbulent, so there must be some gas moving outward in any snapshot as part of turbulent eddies. These “tubulent outflow” are not “real outflow” since they will turn back and join the accretion flow. The main difficulty of obtaining the mass flux of real outflow is therefore how to exclude the “contamination” of turbulent outflow in eq. (2). For this aim, instead of eq. (2), Narayan et al. (2012) move the average inside the integral in eq. (2).
We have adopted a different approach to investigate this problem. We use a “trajectory approach” to analyze the data of GRMHD numerical simulation of accretion flow. Different from the streamline analysis often adopted in accretion flow study, this approach can provide the trajectory of each “virtual test particles” in the accretion flow thus directly show whether the the flow is turbulent outflow or real outflow. The most important result of our analysis is that the mass flux of real outflow is found to be as high as (i.e., eq. (5); refer to Fig. 5). In other words, the mass flux of real outflow is equal to at . As comparison, the mass flux calculated by eq. (2) is equal to at . The reason why Narayan et al. (2012) found a much weaker outflow is that the real outflow is instantaneous. They wander around in 3D space, as shown by Fig. 3, thus will be cancelled if the timeaverage is done first.
Several other important results are as follows.

Most of the real outflow occur in the coronal region of the accretion flow. Within the disk body, it is mainly inflow (Fig. 1).

There are two distinct types of real outflow. One is within the region of away from the axis, and another is outside of this region. In the former region, the poloidal speed of outflow is as high as ; while in the latter region, the speed is much lower, (Fig. 6). We call the outflow in the first region “disk jet” while the outflow in the second region “wind”. “Disk jet” is different from the BlandfordZnajek jet in several ways, as summarized in Yuan & Narayan (2014). The most interesting point is that disk jet exists even though the black hole is nonrotating.

For a given test particle, the poloidal speed of disk jet is found to increase along their trajectory; while for wind, the poloidal speed almost remains constant (Fig. 8). This implies that wind can at least escape beyond the outer boundary of the accretion flows. On the other hand, the poloidal speed of disk jet and wind decreases with increasing radius where they are produced (refer to eq. 8). This implies that the wind has a mixture of poloidal speed depending on its original launching radius. But the mass fluxweighted poloidal speed of wind as a function of radius can be described by eq. (7) (Fig. 7).

The value of Bernoulli parameter of real outflow is not a constant along their trajectories (Fig. 9). The physical reason is that the accretion flow is not steady but turbulent. Because of this reason, the value of for a real outflow is not necessarily positive.

The poloidal speed of outflow does not decrease along the trajectory. This indicates that there must be some acceleration forces. We have analyzed the data and found that in the case of disk jet, the dominant acceleration force is the gradient of magnetic pressure. While for the wind, the centrifugal force and the gradient of gas and magnetic pressure play comparable roles (Fig. 13).
Acknowledgments
FY thanks Jim Stone for the valuable discussions. FY, ZG, and DB were supported in part by the National Basic Research Program of China (973 Program, grant 2014CB845800), the Strategic Priority Research Program “The Emergence of Cosmological Structures” of CAS (grant XDB09000000), and NSF of China (grants 11103059, 11103061, 11121062, and 11133005). RN was supported in part by NSF grant AST1312651 and NASA grant NNX14AB47G. AS acknowledges support for this work by NASA through Einstein Postdoctotral Fellowship number PF4150126 awarded by the Chandra Xray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS803060. XB is supported by NASA through Hubble Fellowship grant HSTHF251301.001A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 526555. The authors acknowledge computational support from NSF via XSEDE resources (grant TG AST080026N), from NASA via the HighEnd Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, and from SHAO Super Computing Platform.
References
 (1) Abramowicz, M. A., Igumenshchev, I. V., Quataert, E., & Narayan, R. 2002, ApJ, 565, 1101
 (2) Aitken, D. K., Greaves, J., Chrysostomou, A., et al. 2000, ApJ, 534, L173
 (3) Balbus, S. A., Hawley, J. F. 1998, Rev. Mod. Phys. 70, 1
 (4) Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1
 (5) Blandford, R. D., & Begelman, M. C. 2004, MNRAS, 349, 66
 (6) Blandford, R.D, & Payne, D.G. 1982, MNRAS, 199, 883
 (7) Blandford, R.D., & Znajek, R.L. 1977, MNRAS, 179, 433
 (8) Begelman, M. C. 2012, MNRAS, 420, 2912
 (9) Bower, G. C., Wright, M. C. H., Falcke, H., Backer, D. C. 2003, ApJ, 588, 331
 (10) Ciotti, L., Ostriker, J. P. & Proga, D. 2010, ApJ, 717, 708
 (11) Crenshaw, D. M., & Kraemer, S. B. 2012, ApJ, 753, 75
 (12) Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003, ARA&A, 41, 117
 (13) Debuhr, J., Quataert, E., Ma, C.P., Hopkins, P., 2010, MNRAS, 406, L55
 (14) Fabian, A. C. 2012, ARA&A, 50, 455
 (15) Gan, Z., Yuan, F., Ostriker, J. P., Ciotti, L., & Novak, G. S. 2014, ApJ, 789, 150
 (16) Gaspari, M., Brighenti, F., & Temi, P. 2012, MNRAS, 424, 190
 (17) Gu, W. M. 2014, ApJ, in press (arXiv:1411.3965)
 (18) Hawley, J. F., & Balbus S. 2002, ApJ, 573, 738
 (19) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
 (20) Hawley, J. F., & Krolik, J. 2006, ApJ, 641, 103
 (21) Igumenshchev, I. V. 2002, ApJ, 577, L31
 (22) Igumenshchev, I. V., Abramowicz, M. A. 1999, MNRAS, 303, 309
 (23) Igumenshchev, I. V., Abramowicz, M. A. 2000, ApJS, 130, 463
 (24) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042
 (25) Kato, Y., Mineshige, S., & Shibata, K. 2004, ApJ, 605, 307
 (26) King, A. 2003, ApJ, 596, L27
 (27) Komissarov, S. S. 2001, MNRAS, 326, L41
 (28) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
 (29) Li, J., Ostriker, J., & Sunyaev, R. 2013, ApJ, 767, 105
 (30) LyndenBell, D. 2003, MNRAS, 341, 1360
 (31) Marrone, D. P., Moran, J. M., Zhao, J. H., Rao, R. 2007, ApJ, 654, L57
 (32) McKinney J. C., 2005, ApJ, 630, L5
 (33) Mou, G., Yuan, F., Bu, D., Sun, M., & Su, M. 2014, ApJ, 790, 109
 (34) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2000, ApJ, 539, 798
 (35) Narayan, R., Sadowski, A, Penna, RF, & Kulkarni, AK. 2012, MNRAS, 426, 3241
 (36) Narayan, R., & Yi I. 1994, ApJ, 428, L13
 (37) Narayan, R., & Yi I. 1995, ApJ, 452, 710
 (38) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69
 (39) Novak, G. S., Ostriker J. P., & Ciotti L., 2011, ApJ, 737, 26
 (40) Ostriker, J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642
 (41) Pen, U.L., Matzner C. D., Wong S., 2003, ApJ, 596, L207
 (42) Penna, R. F., Narayan, R., Sadowski, A., 2013a, MNRAS, 436, 3741
 (43) Penna, R. F., Kulkarni, A., & Narayan, R. 2013b, A&A, 559, A116
 (44) Pringle, J.E., 1981, ARA&A, 19, 137
 (45) Quataert, E., & Gruzinov, A. 2000, ApJ, 539, 809
 (46) Sadowski, A, Narayan, R, Penna, R, & Zhu Y. 2013, MNRAS, 436, 3856
 (47) Shakura, N.I., & Sunyaev, R.A. 1973, Astron. Astrophys, 24, 337
 (48) Shibata, K., & Uchida, Y. 1985, PASJ, 37, 31
 (49) Stone, J. M., Pringle, J. E. 2001, MNRAS, 322, 461
 (50) Stone, J. M., Pringle, J. E., Begelman, M. C. 1999, MNRAS, 310, 1002
 (51) Tchekhovskoy, A., McKinney, J.C., & Narayan, R. 2012, JPhCS, 372, 012040
 (52) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2010, ApJ, 711, 50
 (53) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79
 (54) Wang, Q. D., Nowak, M. A, Markoff, S. B, et al. 2013, Science, 341, 981
 (55) Yuan, F., Bu, D. F., Wu M. C. 2012, ApJ, 761, 130
 (56) Yuan, F., Wu, M. C., Bu D. F. 2012, ApJ, 761, 129
 (57) Yuan, F., Narayan, R. 2014, ARA&A, 52, 529
 (58) Yuan, F., Quataert, E., & Narayan, R. 2003, ApJ, 598, 301