Magneto-hydrodynamical Numerical simulation of wind production from black hole hot accretion flows at very large radii
aa r X i v : . [ a s t r o - ph . H E ] M a r Draft version July 16, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MAGNETO-HYDRODYNAMICAL NUMERICAL SIMULATION OF WIND PRODUCTION FROM BLACKHOLE HOT ACCRETION FLOWS AT VERY LARGE RADII
De-Fu Bu and Feng Yuan , Zhao-Ming Gan , Xiao-Hong Yang Draft version July 16, 2018
ABSTRACTNumerical simulations of black hole hot accretion flows have shown the existence of strong wind.Those works focus only on the region close to black hole thus it is unknown whether or where thewind production stops at large radii. To address this question, Bu et al. (2016) have performedhydrodynamic (HD) simulations by taking into account the gravitational potential of both the blackhole and the nuclear star clusters. The latter is assumed to be ∝ σ ln( r ), with σ being the velocitydispersion of stars and r be the distance from the center of the galaxy. It was found that when thegravity is dominated by nuclear stars, i.e., outside of radius R A ≡ GM BH /σ , winds can no longerbe produced. That work, however, neglects the magnetic field, which is believed to play a crucialdynamical role in the accretion and thus must be taken into account. In this paper, we revisit thisproblem by performing magneto-hydrodynamical (MHD) simulations. We confirm the result of Buet al. (2016), namely wind can’t be produced at the region of R > R A . Our result, combinedwith the results of Yuan et al. (2015), indicates that the formula describing the mass flux of wind˙ M wind = ˙ M BH ( r/ r s ) can only be applied to the region where the black hole potential is dominant.Here ˙ M BH is the mass accretion rate at the black hole horizon and the value of R A is similar to theBondi radius. Subject headings: accretion, accretion disks − black hole physics − hydrodynamics INTRODUCTIONHot accretion flows are very common in the universe,ranging from low-luminosity active galactic nuclei, whichis the majority of nearby galaxies, to the quiescent andhard states of black hole X-ray binaries (see Yuan &Narayan 2014 for a review of our current theoretical un-derstanding of hot accretion flow and its astrophysicalapplications). One of the most important progresses inour understanding of hot accretion flows in recent yearsis the finding of strong wind (Yuan et al. 2012b, hereafterYBW12; Narayan et al. 2012; Li et al. 2013; Yuan et al.2015; Gu 2015). This is an important topic since windis not only an important ingredient of accretion physicsbut also plays an important role in AGN feedback (e.g.,Ostriker et al. 2010). Because of the mass loss via wind,the mass inflow rate decreases inwards ˙ M in ( r ) ∝ r s with s ∼ . − Chandra obser-vation to the supermassive black hole in our Galacticcenter, Sgr A* (Wang et al. 2013).The detailed properties of winds such as the massflux and terminal velocity are calculated by Yuan et al.(2015). This is achieved by performing trajectory analy-sis of virtual test particles based on three dimensionalgeneral relativistic magneto-hydrodynamic (GRMHD)simulation data. Among others, they find that the massflux of wind can be described by ˙ M wind ≈ ˙ M BH ( r/ r s ), Key Laboratory for Research in Galaxies and Cosmol-ogy, Shanghai Astronomical Observatory, Chinese Academyof Sciences, 80 Nandan Road, Shanghai, 200030, China;[email protected] Department of Physics, Chongqing University, Chongqing400044, China with ˙ M BH is the mass accretion rate at the black holehorizon and r s is the Schwarzschild radius. From thisequation, most of the wind comes from the region of largeradius. Then a question is how large the value of r canbe. To investigate this question, Bu et al. (2016, here-after Paper I) have performed HD simulations to studythe accretion flow at very large radii. We take into ac-count the gravity of both the central black hole and thenuclear star clusters. The velocity dispersion of stars isassumed to be a constant and the gravitational poten-tial of the nuclear star cluster ∝ σ ln( r ), where σ is thevelocity dispersion of stars. We find that there is veryfew wind launched from the accretion flow in the regionwhere the gravity is dominated by the star cluster.In Paper I, we introduce an anomalous stress to trans-fer the angular momentum. In reality, it is the MHDturbulence associated with the magneto-rotational insta-bility (MRI; Balbus & Hawley 1991, 1998) that is re-sponsible for the angular momentum transfer. So MHDsimulation is more realistic. Therefore in this paper, wecontinue our study by performing MHD simulations. Asin Paper I, we define a radius R A at which the gravita-tional force due to the central black hole is equal to thatdue to the nuclear star cluster. We call R A to be theboundary of the accretion flow or active galactic nuclei(AGNs). This is partly because the value of R A is in re-ality roughly equal to that of the Bondi radius (Paper I).Hereafter, we use BHAF to refer the hot accretion flowclose to the black hole. We use CAAF (circum-AGN ac-cretion flow) to refer the hot accretion flow at larger radiibeyond R A .The structure of the paper is as follows. In §
2, wedescribe the basic equations and the set up of the initialconditions. The results of simulations will be given in § § METHOD2.1.
Equations
In a spherical coordinate ( r, θ, φ ), we solve the follow-ing magnetohydrodynamical equations describing accre-tion: dρdt + ρ ∇ · v = 0 , (1) ρ d v dt = −∇ p − ρ ∇ ψ + 14 π ( ∇ × B ) × B , (2) ρ d ( e/ρ ) dt = − p ∇ · v + η J , (3) ∂ B ∂t = ∇ × ( v × B − η J ) . (4)Here, ρ , p , v , ψ , e , B and J (= ( c/ π ) ∇ × B ) are den-sity, pressure, velocity , gravitational potential, internalenergy, magnetic field and the current density, respec-tively. d/dt ( ≡ ∂/∂t + v · ∇ ) denotes the Lagrangian timederivative. We adopt an equation of state of ideal gas p = ( γ − e , and set γ = 5 / ψ can be expressed as ψ = ψ BH + ψ star . (5)The black hole potential ψ BH = − GM BH / ( r − r s ), where G is the gravitational constant, M BH is the mass of theblack hole and r s is the Schwarzschild radius. As in PaperI, in this paper, we assume that the velocity dispersionof nuclear stars is a constant of radius. This seems tobe the case of many AGNs. So the potential of the starcluster is ψ star = σ ln( r ) + C , where σ is the velocitydispersion of stars and C is a constant. So we have R A = GM BH /σ . (6)We set σ = 10 and G = M BH = 1 to define our unitsin the present work. Under this units we have R A = 0 . σ ∼ (100 − − (e.g., Kormendy & Ho 2013), R A ∼ (10 − ) r s . Figure1 in Paper I shows the gravitational force distribution.In the above equations, the final terms in Equations(3) and (4) are the magnetic heating and dissipation ratemediated by a finite resistivity η . The exact form of η issame as that used in Stone & Pringle (2001). In the codewe adopt, the energy equation is an internal energy equa-tion, numerical reconnection inevitably results in loss ofenergy from the system. By adding the anomalous resis-tivity η , the energy loss can be captured in the form ofheating in the current sheet (Stone & Pringle 2001).In this paper, time is expressed in unit of the orbitaltime at the torus center.2.2. Initial conditions
As for the initial condition, we assume a rotating equi-librium torus embedded in a non-rotating, low-densitymedium. We assume that the torus has constant specificangular momentum L and assume a polytropic equationof state, p = Aρ γ , where A is a constant. The densitydistribution of torus is ρ = ρ c (cid:26) max[Ψ( R , π/ − ψ ( r, θ ) − L / (2( r sin θ ) ) , A [ γ/ ( γ − (cid:27) / ( γ − , (7) where R is the density maximum (center) of the torus(Nishikori et al. 2006), ρ c is the density at the toruscenter. In this paper, we assume ρ c = 1 and A = 0 . ρ and pressure ρ /r . The mass and pressureof the ambient medius are negligibly small, we choose ρ = 10 − . 2.3. Models
The initial magnetic field is generated by a vector po-tential, i.e. B = ∇ × A . In models A1 and A2, theinitial magnetic field has a dipolar configuration (sameas that in Stone & Pringle 2001). We take A to be purelyazimuthal with A φ = ρ /β , with β = 200. The only dif-ference between models A1 and A2 is that the resolutionof model A1 is two times of that of model A2. We findthat, the results for models A1 and A2 are almost same.So the resolution in model A2 is enough for our problem.In order to study the dependence of results on initialmagnetic field configuration, we carry out model A3. Inthis model, the initial magnetic field has a quadrapolarconfiguration with A φ = ρ /βr cos θ , and β = 100. Table1 summarizes the models.For the initial condition in models A1, A3 and A4,over most of the central regions of the initial torus, wehave 6 grids for one wavelength of the fastest growingmode. Therefore, the fastest growing model of MRI ismarginally resolved in our simulations.In this paper, the simulations are two-dimensional(hereafter 2D). According to the antidynamo theorem(Cowling 1933; Sadowski et al. 2015), the turbulence in-duced by MRI can not be self-sustaining. Therefore therecan not be a true steady state and the quasi-steady stateof the simulations is only transient. In this paper, we stillpreform 2D simulation because on one hand we can sim-ulate a larger radial dynamical range, and on the otherhand previous works have indicated that for many prob-lems the results from 2D and 3D simulations are oftenquite similar (e.g., see the short review in Yuan et al.2012a in the case of radial profile of inflow and outflowrates). Still, for our present study of wind from accre-tion flows, it is necessary to carefully examine whetherthe results from 2D simulation are consistent with thosefrom 3D simulation.In order to answer this question, we have carried outa 2D MHD simulation of accretion flows close to a blackhole. The domain of the simulation is 2 r s -400 r s . In thissimulation, only the black hole gravity is taken into ac-count. Using the trajectory approach as described inYuan et al. (2015), we have calculated the mass flux ofwind, which is ∼
52% of the total outflow rate calcu-lated by Equation (9). As comparison, the mass flux ofwind calculated in Yuan et al. (2015), which is based on3D MHD simulation data of accretion flow, is ∼
60% ofthe total outflow rate calculated by Equation (9). Sucha good consistency suggests that the present 2D studyshould be a good approximation.2.4.
Numerical Method
We use the ZEUS-2D code (Stone & Norman 1992a,1992b) to solve Equations (1)-(4). The polar range is 0 ≤ θ ≤ π . We adopt non-uniform grid in the radial direction( △ r ) i +1 / ( △ r ) i = 1 . θ wo-dimensional magnetohydrodynamical case 3 TABLE 1models in this paper
Models Initial magnetic field Resolution Computational domainA1 dipolar 294 ×
160 0.03-4A2 dipolar 147 ×
80 0.03-4A3 quadrapolar 294 ×
160 0.02-4A4 dipolar 294 ×
160 0.002-0.4 direction in the northern and southern hemispheres aresymmetric about the equatorial plane. The resolution at θ in the northern hemisphere is same as that at π − θ in the southern hemisphere. In order to well resolve theaccretion disk around the equatorial plane, the resolutionis increased from the north and south rotational axis tothe equatorial plane with ( △ θ ) j +1 / ( △ θ ) j = 0 . ≤ θ ≤ π/ △ θ ) j +1 / ( △ θ ) j = 1 . π/ ≤ θ ≤ π . At the poles, we use axisymmetric boundaryconditions. At the inner and outer radial boundary, weuse outflow boundary conditions. RESULTS3.1.
Mass inflow rate
Following Stone et al. (1999), we define the mass inflowand outflow rates, ˙ M in and ˙ M out as follows,˙ M in ( r ) = 2 πr Z π ρ min( v r ,
0) sin θdθ, (8)˙ M out ( r ) = 2 πr Z π ρ max( v r ,
0) sin θdθ. (9)The net mass accretion rate is,˙ M acc ( r ) = ˙ M in ( r ) + ˙ M out ( r ) . (10)Note that the above rates are obtained by time-averagingthe integrals rather than integrating the time averages.Figure 1 shows the time-averaged (from 130-136 orbits)and angle-integrated mass rates of model A1. Both themass inflow and outflow rates decrease inward. This isconsistent with that found in the HD simulations in Pa-per I. We note that the mass inflow and outflow rates arenot good power-law function of radius, just like the caseof a BHAF (Stone & Pringle 2001). This is likely be-cause the radial distribution of the strength of magneticfield is not very smooth. Because of the accumulation ofmagnetic flux during the simulation, the magnetic fieldin the inner region of the accretion flow is much strongerthan the other region. We will see in § r ∼ .
6. We note that the flow in model A1 is convec-tively unstable (see section 3.3). The accretion timescaleis roughly equal to one turnover time of local convec-tive eddies. It may require many turnover times for theconvection to reach a steady state. Thus it will be inter-
Fig. 1.—
The radial profile of the time-averaged (from t = 130 to136 orbits) and angle integrated mass inflow rate ˙ M in (solid line),outflow rate ˙ M out (dashed line), and the net rate ˙ M acc (dottedline) in model A1. They are defined in Equations (8), (9), and(10), respectively. Fig. 2.—
Snapshot of velocity vector of model A1 at t = 130orbits. Clearly, turbulent eddies occupy the whole domain and itis hard to find winds. esting to run our simulations for several times longer inthe future to check whether the results will change.3.2. Does strong wind exist in a CAAF ?
The significant mass outflow rate shown in Figure 1does not mean the existence of strong real outflow (wind)because it may be due to the turbulent motion. To studywhether wind exists, let us first directly look at the ve- Θ ( d e g r ee ) R −0.4−0.20.00.20.4 z Fig. 3.—
Trajectory of gas for model A1. The black dots locatedat r = 0 . θ angle. It is clear that the “test particle” crosses the startingradius for many times. From this figure we can see that the realwind trajectories, i.e., the trajectories which extend from r = 0 . r = 0 . locity field shown in Figure 2. We see that turbulenteddies occupy the whole domain and it is hard to findsystematic winds.To investigate this problem precisely, following Yuanet al. (2015), we use the trajectory method to studythe motion of the virtual particles. The details of thisapproach can be found in Yuan et al. (2015). To getthe trajectory, we first need to choose some virtual “testparticle” in the simulation domain. They are of coursenot real particles, but some grids representing fluid ele-ments. Their locations and velocities at a certain time t are obtained directly from the simulation data. We canthen obtain their location at time t + δt from the velocityvector and δt . Trajectory is more loyal than streamlineto reflect the motion of particles which is crucial for us toinvestigate whether real wind exists. Trajectory is onlyequivalent to the streamline for strictly steady motion,which is not the case for accretion flow since it is alwaysturbulent.Trajectory approach can easily tell us which particlesare real outflows (i.e., winds) and which are doing turbu-lent motions. Combining with the density and velocityfield information from the simulation data we can thenobtain the various properties of wind such as the massflux, angular distribution, and velocity (see Yuan et al.2015 for details). Figure 3 shows the trajectory of some gas particles starting from r = 0 . r = 0 . r = 0 . r = 0 . . Why the inflow rate decreases inward in a CAAF?
If wind is absent, what is the reason for the inwarddecrease of inflow rate? To answer this question, we needto examine the convective stability of the accretion flow.Hydrodynamical simulation of BHAFs (e.g., Stone et al.1999; Igumenshchev & Abramowicz 1999, 2000; Yuan& Bu 2010) have found that the flows are convectivelyunstable, consistent with what has been suggested bythe one-dimensional analytical study of BHAFs (Narayan& Yi 1994). The physical reason is that the entropyof the flow increases inward, which is resulted by theviscous heating and negligible radiative loss. However,in the presence of magnetic field, numerical simulationshave found that a BHAF becomes convectively stable(Narayan et al. 2012; YBW12). In the case of a CAAF,Paper I has found that the flow is convectively unstable,same as the case of a BHAF.We now study whether a CAAF is convectively stableor not in the presence of magnetic field. We use theHøiland criteria (e.g., Tassoul 1978; Begelman & Meier1982):( ∇ s · dr )( g · dr ) − γv φ R [ ∇ ( v φ R ) · dr ] dR < . (11)In Equation (11), R = r sin θ is the cylindrical radius, dr = dr ˆ r + rdθ ˆ θ is the displacement vector, s = ln( p ) − γ ln( ρ ) is ( γ −
1) times the entropy, g = −∇ ψ + ˆ Rv φ /R isthe effective gravity, and v φ is the rotational velocities.For a non-rotating flow, this condition is equivalent toan inward increase of entropy, which is the well-knownSchwarzschild criteria.Taking model A1 as an example, Figure 4 shows theresult. The result is obtained according to Equation (11)based on the simulation data at t=132 orbits at the ini-tial torus center r = 1. At t = 132 orbits, the flow hasachieved a steady state since the net accretion rate av-eraged between t = 130 to 136 orbits is a constant ofradius (see figure 1, the dotted line). The red regions areconvectively unstable. We can see from the figure that aCAAF is mostly convectively unstable. This is differentfrom the case of a BHAF with magnetic field. The reasonshould be due to the change of the gravitational potentialbut the detail remains unclear. This result strongly im-plies that the inward decrease of inflow rate is because ofconvection and it reminds us the scenario of convection-dominated accretion flow (CDAF) proposed by Narayanet al. (2000) and Quauaert & Gruzinov (2000), althoughthat model was proposed to explain the dynamics ofBHAFs rather than CAAFs.wo-dimensional magnetohydrodynamical case 5 Fig. 4.—
Convective stability analysis of Model A1. The result isobtained according to Equation (11) based on the simulation dataat t=132 orbits at the initial torus center r = 1. The red region isunstable. Fig. 5.—
Convective stability analysis of Model A4. The result isobtained according to Equation (11) based on the simulation dataat t=100 orbits at the initial torus center r = 0 .
1. The red regionis unstable.
From Figure 4, it seems that the region r < . r < .
1, itis true that the black hole gravity is bigger than that ofthe star cluster, but the gravity of the stars is not neg-ligible. In fact, in the region 0 . < r < .
1, the blackhole gravity is bigger than that of star cluster at most bya factor of 4. To further investigate this point, we haveanalyzed the convective stability of model A4. Figure 5shows the results of the region very close to the blackhole where the black hole gravity is strongly dominant.We can clearly see that the flow is convectively stable.3.4.
Varying the initial configuration of magnetic field
In order to study whether the results depend on theinitial configuration of the magnetic field, we carry outmodel A3. In this model, the initial configuration ofmagnetic field is quadrapolar. Figure 6 shows the inflow,
Fig. 6.—
Same as Figure 1, but for model A3. The time-averageis taken from t = 127 to 130 orbits. outflow, and net rates. The inflow rate is smoother thanthat in Figure 1, as we have explained in § r = 0 . . Moving the computational domain inward
From Figures 1 and 6, in the region r < .
1, the out-flow rate is very small. This region is dominated by thegravity of black hole. Previous works (Yuan et al. 2012b;Yuan et al. 2015) have shown that in this case outflow(wind) should be strong. The apparent discrepancy be-tween this work and our previous works is due to the factthat the region r < . . < r < .
4. The gravitational potential from boththe black hole and the nuclear star cluster are includedbut obviously the former dominates. The initial condi-tion of the magnetic field is dipolar. Figure 7 shows thetrajectory of some gas particles starting from r = 0 . Θ ( d e g r ee ) R −0.08−0.06−0.04−0.020.000.020.040.060.08 z Fig. 7.—
Trajectory of gas for model A4. The black dots locatedat r = 0 .
03 are starting points of the “test particle”. Differentcolors denote trajectory of “test particle” starting from different θ angle. It is clear that wind is present when the black gravitydominates. CONCLUSION AND DISCUSSIONNumerical simulations show that strong winds exist inblack hole hot accretion flow (e.g., YBW12). The massflux of wind follows ˙ M wind = ˙ M BH ( r/ r s ) (Yuan et al.2015). A question is then what the value of r can be, i.e.,whether or where the wind production stops. In orderto answer this question, in Paper I, we have performedHD simulations and take into account the gravity of boththe black hole and nuclear stars cluster. We find that themass inflow rate decreases inward. However, our trajec-tory analysis indicates that there is no wind when thepotential of star cluster dominates, i.e., beyond a certain radius R A ≡ GM BH /σ , with σ being the velocity dis-persion of stars. The inward decrease of inflow rate is notbecause of strong wind, as in the case of accretion flowin which the black hole potential dominates, but becauseof the convective instability of the accretion flow. In thispaper, we revisit the same problem by performing morerealistic MHD simulations. We find again that the con-clusion remains unchanged, i.e., there is no wind beyond R A . Our stability analysis again indicates that the MHDaccretion flow beyond R A is convectively unstable. Thisis different from the case of accretion flow when the blackhole potential dominates (YBW12; Narayan et al. 2012).So the inward decrease of inflow rate is likely because ofthe convective motion of the flow.This result indicates that the mass flux of wind foundby Yuan et al. (2015):˙ M wind = ˙ M BH ( r/ r s ) , (12)can only be applied to the region where the black holegravitational force dominates. In the star cluster poten-tial dominates region, i.e., beyond R A , no wind will beproduced. In practice, the value of R A is close to theBondi radius R B ≡ GM BH /c s (Paper I).What is the reason for the absence of wind beyond R A ?We speculate that it may be related with the change ofthe slope of the gravitational potential. Such a changewould change the shear of the accretion flow and then theturbulent stress in the accretion flow. Analytical modelsof accretion disk (e.g., Shakura & Sunyaev 1973) usuallyassume that viscous stress is proportional to the shear ofthe accretion disk, T rφ = ρνr d Ω dr (13)with ν = αc s / Ω k . c s , Ω and Ω k are sound speed, angularvelocity and Keplerian angular velocity, respectively. Re-cent 3D MHD numerical simulations show that the tur-bulent stress is not linearly proportional to the shear, butthe dependence is stronger than that predicted by Equa-tion (13) (Pessah et al. 2008; Penna et al. 2013). Thisimplies that the change of the potential changes someproperties of turbulence which then change the wind pro-duction. ACKNOWLEDGEMENTSWe thank the referee for rasing very useful questionswhich improve the paper significantly. We thank RameshNarayan for helpful discussions. This work was sup-ported in part by the National Basic Research Pro-gram of China (973 Program, grant 2014CB845800),the Strategic Priority Research Program The Emergenceof Cosmological Structures of the Chinese Academy ofSciences (grant XDB09000000) and the Natural Sci-ence Foundation of China (grants 11103061, 11133005,11121062 and 11573051). This work has made use ofthe High Performance Computing Resource in the CoreFacility for Advanced Research Computing at ShanghaiAstronomical Observatory. REFERENCESBalbus S. A., Hawley J. F., 1991, ApJ, 376, 214 Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70,1 wo-dimensional magnetohydrodynamical case 7wo-dimensional magnetohydrodynamical case 7