Spherical Accretion with Anisotropic Thermal Conduction

# Spherical Accretion with Anisotropic Thermal Conduction

## Abstract

###### keywords:
convection, conduction – magnetic fields, MHD – methods: numerical

## 1 Introduction

In magnetized plasmas with a collisional mean free path much greater than the Larmor radius, there is appreciable heat and momentum transport along magnetic field lines but negligible transport across the field lines (Braginskii, 1965). Balbus (2000) found that anisotropic thermal conduction fundamentally alters the convective stability of dilute stratified plasmas. For a highly collisional fluid, convection sets in when the entropy increases in the direction of gravity. A low collisionality plasma, however, is unstable to convective-like motions regardless of the sign of the temperature gradient (Balbus, 2000; Quataert, 2008). For a temperature increasing in the direction of gravity, horizontal fields are the most unstable (the magnetothermal instability; MTI), while when the temperature decreases in the direction of gravity, vertical fields are the most unstable (the heat flux driven buoyancy instability; HBI). The MTI & HBI have been studied using local simulations by Parrish & Stone (2005, 2007) and Parrish & Quataert (2008), respectively.

In this paper, we study the effects of anisotropic thermal conduction and the MTI on the global dynamics of spherical accretion flows using axisymmetric magnetohydrodynamic (MHD) simulations. Thermal conduction is particularly important for hot low density plasmas in which the electron mean free path is an appreciable fraction of the size of the system, as in, e.g., the accretion flow onto Sgr A* in our Galactic Center and accretion flows onto massive black holes in elliptical galaxies (Tanaka & Menou, 2006).

There are several motivations for studying anisotropic conduction in spherical accretion flows. First, the MTI is predicted to be present at large radii in spherical accretion flows, for (where is the gravitational sphere of influence of the central object and is the isothermal sound speed of the ambient gas). At these radii, the inflow is subsonic and there is sufficient time for the instability to grow (Balbus, 2000). At smaller radii (), the MTI is unlikely to grow significantly because the inflow is supersonic (at least according to the classic Bondi 1952 solution). However, if there is resistive dissipation of the magnetic energy at small radii, conduction may nonetheless be dynamically important because it provides an efficient way to transport energy to larger radii, which can suppress the accretion rate onto the central object (Johnson & Quataert, 2007). As we discuss in more detail later in the paper, reduction of accretion rate does not specifically require anisotropic heat transport, but is in fact analogous to the suppression of accretion by magnetic dissipation seen in the resistive MHD simulations of Igumenshchev & Narayan (2002).

A direct observational motivation for studying the MTI in spherical accretion flows is provided by mm observations of Sgr A*. These find a rotation measure (RM) of rad m that has been roughly constant since it was first measured over years ago (Aitken et al., 2000; Bower et al., 2005; Marrone et al., 2007). The stability of the observed RM is difficult to explain if the RM is produced by the accretion flow at small radii (Sharma, Quataert, & Stone, 2007). In local simulations of the MTI, the field is amplified by factors of and becomes aligned with the temperature gradient (Parrish & Stone, 2007). Thus the MTI at large radii in spherical accretion flows is expected to produce a primarily radial magnetic field that could help explain the stable measured RM towards Sgr A* (a possibility we speculated on in Sharma, Quataert, & Stone 2007).

The remainder of this paper is organized as follows. We describe our numerical set up in §2. Because it is not computationally feasible to directly simulate both the small radii close to the central object where most of the gravitationally potential energy is released and the radii where the MTI develops, we carry out two different sets of simulations in this paper, covering different radial scales in the accretion flow. In §3 we present the results of simulations at that study the development of the MTI for weak initial magnetic fields (see Table 1). In §4 we describe simulations at that study the effect of thermal conduction and magnetic dissipation on spherical accretion for reasonably strong initial magnetic fields (see Table 2). We then conclude and discuss the astrophysical implications of our results (§5). Throughout this paper we consider the simplified problem of purely spherical accretion; in future studies, we will examine the effects of anisotropic thermal conduction and the MTI in rotating accretion flows.

## 2 Numerical Simulations

We use the axisymmetric ZEUS-2D MHD code in spherical geometry (Stone & Norman, 1992a, b). We solve the MHD equations with anisotropic thermal conduction:

 ∂ρ∂t+∇⋅(ρV)=0, (1) ∂V∂t+V⋅∇V=−∇(p+Π)ρ+J×Bρ−∇Φ, (2) ∂B∂t=∇×(V×B−ηJ), (3) ∂e∂t+∇⋅(eV)=−(p+Π)∇⋅V−∇⋅Q+ηJ2, (4) Q=−χ^b(^b⋅∇)T, (5)

where is the mass density, is the fluid velocity, is the pressure, is an explicit artificial viscosity (e.g., see Stone & Norman, 1992a), is the current density, is the gravitational potential due to a central mass (we use the pseudo-Newtonian potential of Paczynski & Wiita 1980), , is the magnetic field strength, is the explicit resistivity used in some of our simulations, is the internal energy density ( is the adiabatic index), is the heat flux along field lines, is the thermal diffusivity, is the unit vector along , and is the temperature. For convenience we often use which has the dimensions of a diffusion coefficient (cm s). Anisotropic thermal conduction is added in an operator split fashion with subcycling (e.g., see Parrish & Stone 2005) and is implemented using the method based on limiters that guarantees that heat flows from hot to cold regions (see Sharma & Hammett, 2007).

We initialize our simulations using the hydrodynamic Bondi solution with a constant vertical magnetic field . Note that the initial density and temperature gradients are in the radial direction; thus, the initial field is not everywhere perpendicular to the initial temperature gradient, unlike in previous local studies of the MTI. We consider both weak and strong fields to test the effect of the strength of the field on the flow dynamics. The field strength is specified by the dimensionless number , defined by , where is the outer radius of the computational domain and . This is equivalent to , where is the ratio of gas pressure to magnetic pressure.

We use an adiabatic index so that the initial solution has a sonic point that lies within the computational domain for our simulations at (Table 1).1 For our simulations at large radii, we use the Bondi radius, , to normalize length scales in the problem, while for the simulations at small radii, we use the gravitational radius of the central object, . We express time in units of inverse rotation frequency at the outer radius of the domain, .

The simulations at described in §3 (Table 1) go from to . For these simulations, we keep the temperature at the outer radial boundary fixed at its initial value of , i.e., K. The simulations described in §4 (Table 2) go from to . The temperature is fixed at the virial temperature at the outer boundary (except for one simulation with a very weak initial magnetic field [S2]; in this simulation there is no significant change in the flow properties relative to the hydrodynamic Bondi solution, so fixing the temperature to be that of the initial Bondi solution is more physical). For the simulations at (§3), fixing the temperature at the outer boundary is reasonable because of the large thermal inertia of the plasma at large radii. For the simulations from to described in §4, however, the choice of the outer temperature boundary condition is more subtle and has a nontrivial effect on the solution. Because we initialize the flow with a solution, the initial temperature is significantly less than the virial temperature. Fixing to its initial value is inappropriate because dissipation leads to a temperature throughout most of the domain and thus setting equal to its initial value leads to an unphysically large temperature gradient at . Setting is unphysical because this sets the heat flux to 0. An alternative boundary condition that we considered, setting a constant heat flux across the outer boundary (corresponding to constant), led to a roughly constant temperature throughout the domain with after several ; this drove all the mass out of the simulation domain. We thus settled on the physically motivated outer boundary condition of . Factor of few variations in the exact value of did not significantly change the results of the simulations. In both sets of simulations, inflow/outflow boundary conditions are applied at the radial boundaries.

The radial grid is logarithmically spaced (grid spacing ) so that the number of grid points from to is roughly the same as the number of grid points from to . A logarithmic grid gives reasonable resolution at all radii; e.g., for our simulations at small radii (), the cell size at the smallest radii is . A uniform polar grid extends from to . Reflective boundary conditions are applied at the poles. Our standard resolution is , but we have also carried out simulations at higher resolution to assess convergence. Our calculations use a version of ZEUS that is not parallel, so that even these modest 2D calculations are time intensive.

We do not use an explicit resistivity for the simulations in §3 (Table 1). These simulations are primarily aimed at studying the MTI for weak fields in which case the magnetic energy is so small that using an explicit resistivity is not required to correctly capture the dynamics (in addition, we have verified that using the explicit resistivity below does not change the results). However, for (§4; see Table 2), the gravitational potential energy of the inflowing plasma is converted into magnetic energy by flux freezing. The inflow brings together oppositely directed field lines, which leads to magnetic dissipation and plasma heating. This energy is then transported to large radii by thermal conduction. To capture this dynamics, it is important to conserve energy as accurately as possible. Although we are not using a conservative code, we can capture a reasonable fraction of the dissipated magnetic energy using an explicit resistivity of the form (Stone & Pringle, 2001; Igumenshchev & Narayan, 2002)

 η=η0dr2|∇×B|√4πρ (6)

with (with this choice for , we find that total energy is conserved to better than 30 % in our simulations at ). The resistive terms in the induction and internal energy equations are included using the method of Fleming, Stone, & Hawley (2000).

Although our simulations are two-dimensional, the turbulent magnetic and velocity fields do not decay at late times, because the anti-dynamo theorem does not apply when a net magnetic flux threads the simulations domain. In this case the source current for the magnetic field is outside of the simulation domain (Balbus & Hawley, 1998).

### 2.1 Thermal Conductivity

In a collisional plasma, the thermal conductivity is given by cm s (Spitzer, 1962). However, the inner regions of hot accretion flows are, in many cases, collisionless with the electron mean free path due to Coulomb collisions larger than the radius (e.g., Tanaka & Menou 2006; Sharma et al. 2007); the collisional result is thus inapplicable. In a collisionless plasma, electrons can in principle carry a free streaming heat flux as large as ( is the thermal speed; Cowie & McKee 1977; Snyder, Hammett, & Dorland 1997). The diffusion coefficient in this case is , assuming a virial temperature proportional to . Motivated by the application to collisionless systems such as Sgr A*, we take as the thermal conductivity for our simulations and define the dimensionless conductivity by

 κ≡αc(GMr)1/2. (7)

Although the scaling is well-motivated, the precise value of is difficult to calculate. It depends on the electron mean free path in the collisionless limit, which is determined by the rate of pitch angle scattering by small-scale (high frequency) turbulent fluctuations. Motivated by the results of Sharma et al. (2006) and Sharma et al. (2007) on pitch angle scattering in collisionless accretion flows, we take , i.e., , as our fiducial diffusion coefficient. This value is much smaller than the maximum free-streaming thermal conductivity that could be obtained if the electrons are virial. Our choice for is such that the large-scale thermal conduction timescale () is comparable to the free-fall timescale (). Thermal conduction can thus lead to a significant modification of the thermal structure of the plasma (as we shall see). Since both and scale as , this is true at all radii. In addition, with on large scales, small-scale perturbations are effectively isothermal and the MTI grows at its maximal rate (Balbus, 2000). For the MTI simulations at , we carried out simulations with factors of few smaller and larger than our fiducial choice of , and found results essentially identical to those described here. For the simulations at , resistive heating and thermal conduction significantly modify the dynamical properties of the accretion flow; in this case, we show solutions for a variety of values of to explicitly show the dependence on the conductivity (which is still relatively modest).

As we discuss in §3, it turns out that a conductivity has the property that for the initial Bondi solution near . Since the initial condition is already approximately a steady state solution even in the presence of conduction (modulo the MTI), including conductivity does not affect the structure of the plasma. We thus carried out simulations at with other choices for (namely, and ). Although we do not show detailed results for these simulations, in all cases we found that the MTI rapidly creates radial magnetic field lines and that conduction along the radial magnetic field lines forces the density and temperature structure of the plasma to readjust in order to satisfy . This can lead to modest modifications in the density and temperature profiles of the flow, but in all cases the accretion rate was within a factor of a few of the Bondi rate.

## 3 Simulations at r∼rB

Figure 1 shows several different timescales characterizing the initial Bondi solution. For the infall timescale () is much longer than the timescale for the MTI to grow () and the instability grows exponentially. The MTI can thus saturate before the plasma flows in. At smaller radii (), however, where the inflow time is comparable to the MTI growth time, the growth is algebraic and the MTI is less likely to have a significant effect on the dynamics. Instead, at , the magnetic field will be primarily amplified by flux freezing. Here we discuss our simulations of the MTI at . In the next section we present the results of simulations at .

### 3.1 The Fiducial MTI Simulation

Our fiducial simulation (labeled “R1”) uses and has an initial magnetic field strength of , corresponding to at . This simulation describes the evolution of the MTI for weak magnetic fields. We choose the lower resolution simulation R1 as our fiducial case because it finishes relatively quickly (in a week); this enables us to readily study the effects of different boundary conditions, different choices for the thermal conductivity, etc. We also present higher resolution calculations (e.g., D1 & Q1, which take several months to finish) that are similar to R1, although with a modest increase in the field amplification by the MTI (see Fig. [9]).

Figure 2 shows the and components of the magnetic energy density as a function of time for the fiducial simulation and for an MHD simulation without conduction, but with the same initial conditions. The energy is averaged from to , to isolate the radii where (see Fig. 1) and the instability is unaffected by inflow. The top panel in Figure 2 shows the early linear growth of the MTI, while the bottom panel shows the longer time evolution. Because we are averaging over a range of radii, there is no unique growth rate for the MTI, but the maximal MTI growth rate evaluated at provides a reasonable fit to the early time growth (long-dashed line). After a few free-fall times (several MTI growth times; see Fig. 1), the magnetic field ceases to grow exponentially and the growth becomes algebraic. This represents the saturation of the MTI; the later growth of the field is due to flux freezing. This can be seen explicitly by noting that the MHD simulation (short-dashed line) shows the same algebraic growth in the magnetic field with time, only at a reduced magnitude of because there is no early MTI growth. Note that the net increase in magnetic energy because of the MTI is quite modest: even in the saturated state.

Figure 3 compares the plasma at the end () of the fiducial run and the corresponding MHD simulation. The magnetic field unit vectors are also shown. The plasma with conduction is a factor of smaller than without conduction (as is also seen in Fig. 2), indicating that the MTI alone amplifies the field by a factor of . Since the simulations do not run for even one infall timescale at , the magnetic field lines at large radii do not change significantly in the MHD simulation and remain primarily vertical. By contrast, the magnetic field is completely restructured in the presence of anisotropic conduction because of the MTI. The magnetic field at the end of the simulation is primarily radial. This can also be seen in Figure 4, which shows the volume averaged angle of the magnetic field with respect to the radial direction as a function of time for several runs. In the linear stage of the MTI, the radial field is amplified exponentially; after saturation, the field direction, although turbulent, continues to become radial because of flux freezing (as also occurs in the pure MHD simulation shown in Fig. 4), although this occurs on the slower inflow timescale. Our conclusion that the MTI generates primarily radial field lines in the saturated state is fully consistent with the local MTI simulations of Parrish & Stone (2005, 2007).

To understand the saturation of the MTI, it is useful to consider the linear dynamics of the instability. In the WKB limit, and assuming rapid conduction and weak magnetic fields, the linear growth rate is given by (Quataert, 2008)

 γ2≈1ρdpdr(dlnTdr)[(1−2b2r)k2θk2+2bθbrkθkrk2], (8)

where and are the and wavenumbers, and and are the and magnetic field unit vectors. Our initial field is vertical with a dependent and . In the saturated state, where and , the fastest growth rate of the MTI is , a factor of smaller than the growth rate in the initial configuration. Growth occurs only for modes satisfying , i.e., for small radial wavelengths. The significant decrease in the growth-rate of the MTI for nearly radial fields qualitatively explains why magnetic field reorientation leads to saturation of the instability. In addition, even if there is growth at late times, it is on small radial scales and so is unlikely to modify the large-scale field structure or magnetic energy. The MTI can also saturate by making , i.e., by making the plasma isothermal. Whether the plasma becomes isothermal or not depends on (and on the boundary conditions), which influences how the density and temperature adjust to make in steady state. Figure 5 shows that for our fiducial simulation, there is very little change in the temperature and density profiles during the simulation. For our choice of , it turns out that the initial Bondi flow nearly satisfies near . Thus when the MTI reorients the field lines to be primarily radial, little evolution in and is required to make . As noted in §2.1, we carried out simulations with different in order to assess its effects on the flow dynamics (although we believe that the fiducial choice considered here is the most physical for collisionless systems). In all cases, we found that the MTI reorients the field lines to be primarily radial and that and adjust so that .

Figure 6 shows the time and angle averaged mass accretion rate in the saturated state as a function of radius. In MHD the velocity is always radial and inwards and the mass accretion rate is equal to the Bondi rate (for this simulation). The solid line in Figure 6 shows the net mass accretion rate for the fiducial run, ; also shown are the mass fluxes of material with and . Although there are MTI-driven convective motions, they are modest and the net accretion rate is only slightly smaller than the Bondi accretion rate. This is partially a consequence of the fact that the density and temperature profiles change very little in the course of the simulation (Fig. 5). However, even for different choices of and , we still found for MTI simulations with initially weak fields. The MTI by itself thus does not significantly change the accretion rate onto the central object. This is in contrast to the very strong suppression of found in the next section by the combined action of magnetic dissipation, plasma heating, and thermal conduction at .

Figure 7 shows time and angle averaged luminosities ( energy fluxes) for energy transport by conduction and MTI-initiated convection. The luminosities are normalized to the field free (“Spitzer”) conductive luminosity at . The radial conductive luminosity () is roughly constant as a function of radius for , as would be expected in steady state. The conductive heat flux is similar to the local field-free value at all radii because the field lines are primarily radial. The non-conductive energy transport has contributions from the kinetic energy flux (), the enthalpy flux (), and the Poynting flux. The latter is negligible for the weak magnetic field values considered here. Both the kinetic energy flux and the enthalpy flux can be split up into terms proportional to the mass accretion rate (the advective flux of energy due to the bulk motion of the fluid) and terms which represent the turbulent/convective transport of energy (). Mathematically, this can be seen by decomposing the fluid variables into mean and fluctuating components; e.g., the radial velocity is taken to be , where represents an appropriate average and represents deviation from the average. For our purposes, the average is taken over at each time at a given . With this decomposition, the average kinetic energy flux is given by

 12⟨ρV2Vr⟩=QK,a+QK,t,

with

 QK,a = 12⟨ρVr⟩[⟨Vr⟩2+⟨Vθ⟩2+⟨δV2r⟩+⟨δV2θ⟩], (9) QK,t = 12⟨δ(ρVr)(δV2r+δV2θ)⟩+⟨δ(ρVr)δVr⟩⟨Vr⟩ (10) + ⟨δ(ρVr)δVθ⟩⟨Vθ⟩,

where is the advective flux of kinetic energy and is the turbulent/convective kinetic energy flux. Similarly, the average enthalpy flux can be written as

 γγ−1⟨pVr⟩=Qe,a+Qe,t,

with

 Qe,a = γγ−1⟨ρVr⟩⟨c2s⟩, (11) Qe,t = γγ−1[⟨Vr⟩⟨δρδc2s⟩+⟨ρ⟩⟨δVrδc2s⟩ (12) + ⟨δρδVrδc2s⟩],

where . Figure 7 shows the total turbulent energy flux in the MTI simulations at large radii, ; the turbulent energy transport – which is dominated by the term in equation (12) – is much smaller than the conductive energy transport at all radii.

### 3.2 Effect of B0

If the MTI saturates by making the magnetic field lines primarily radial, as we showed in the previous section, it is reasonable to expect that the saturation magnetic field strength will scale with the initial magnetic field strength (since a smaller radial field is required to make an initially weaker vertical magnetic field radial). To test this hypothesis we carried out the simulation R2 with a weaker vertical magnetic field (, corresponding to at . All other parameters are same as the fiducial run.

Figure 8 shows the results from run R2 compared to those of R1. The most interesting result is that even nonlinearly for R1 and R2 are similar. That is, the late time magnetic energy is roughly proportional to the initial magnetic energy in the domain. This indicates that the saturation of the MTI is quasilinear: for weak magnetic fields, the MTI saturates when the field lines become radial, which requires a .

We also carried out simulations with much stronger initial magnetic fields, , corresponding to at (run R3 with conduction and run MHD2 without conduction). The magnetic field is sufficiently strong that the MTI is suppressed by magnetic tension. Indeed, we find that there are only small differences between the simulations with and without conduction. In both cases the magnetic field at small radii is amplified by the converging inflow. Magnetic tension prevents plasma infall in the equatorial region; mass accretion is dominated by the polar regions. There is a strong equator-pole anisotropy because of the strong magnetic field. The polar regions are cold, low density, and magnetically dominated (), with large radial infall velocities. We find a modest reduction in the accretion rate for these strong-field simulations at , with the mass accretion rate through being of the Bondi rate. Finally, we find that including an explicit resistivity does not have a significant effect on the flow dynamics at large radii even for energetically significant magnetic fields. This is in contrast to the simulations described in §4 in which there is significant magnetic dissipation when at . The difference is probably that the supersonic inflow at is much more effective than the slow inflow velocities at at forcing reconnection/dissipation of field lines.

### 3.3 Effect of Resolution

To assess the effect of resolution we carried out several higher resolution runs: run D1, which is analogous to the fiducial simulation R1 but with twice the resolution, run D2 which is analogous to R2 but with twice the resolution, and run Q1 which has four times the resolution of R1. While the overall dynamics is very similar to the lower resolution runs, the magnetic energy is a factor of few larger in the higher resolution simulations. Figure 9 compares the fiducial run R1 with D1 and Q1, and shows that both and are larger in higher resolution runs than they are in R1, by a factor of . The same is true for D2 compared with R2. In particular, Figure 9 shows that the phase of linear exponential magnetic field amplification occurs for a longer period of time in the higher resolution simulations. Since the late-time magnetic field growth is simply flux freezing of the early-time field, this leads to an enhanced field at all times. The fastest growing MTI modes for approximately radial magnetic fields are those with , which is a resolution dependent criterion. Higher resolution simulations can thus resolve growth for longer times and to larger radial fields, i.e., to smaller .

Although the fastest growth of the MTI occurs at small scales, nonlinearly one expects the energy contained in small scales to be subdominant for sufficiently high resolution simulations. Indeed, Figure 9 shows that both the radial and polar magnetic energies are similar for the higher resolution runs D1 and Q1. This indicates that nonlinear saturation of the MTI is reasonably independent of resolution for sufficiently high resolution simulations (e.g., runs D1 and Q1).

## 4 Simulations at r≪rB

As shown in the previous section, the MTI generates a radial field at in dilute spherical accretion flows. As the plasma flows in to yet smaller radii, this field will be amplified by flux freezing. It is well known that for inflow according to the Bondi solution (Bondi, 1952), the magnetic energy increases more rapidly with decreasing radius than the gravitational energy of the gas, and eventually the field becomes energetically important (e.g., Shvartsman 1971). Igumenshchev & Narayan (2002) showed that when magnetic dissipation of the amplified field becomes significant and the dynamics of the inflow changes appreciably. Resistive heating increases the temperature of the plasma and may drive thermal convection. The accretion rate onto the central object decreases by a few orders of magnitude. In this section, we examine the modification to these results caused by anisotropic thermal conduction. Our general conclusion that mass accretion is significantly suppressed is consistent with that of Igumenshchev & Narayan (2002), but we find that the conductive transport of energy is more important than the convective transport.

Our simulations at small radii have and ; their properties are summarized in Table 2. The initial conditions are chosen from the same Bondi solution used in the previous section. Run S1 takes and . This value of corresponds to at . Igumenshchev & Narayan (2002) have shown that a weaker field at means that significant dissipation does not begin until smaller radii (and later times), but that otherwise the results do not depend sensitively on the choice of . In addition to the fiducial value of , we also performed simulations with and ; the latter corresponds roughly to the maximum rate of heat transport by free-streaming electrons. We present results for a variety of but focus on the simulation (S3) because it most clearly demonstrates the relevant physics.

S2 and MHD3 represent control simulations. S2 is a simulation with anisotropic conduction but with a very weak field. As expected, it is effectively a laminar hydrodynamic Bondi solution (see Fig. 14), with no evidence for the MTI or dynamically significant resistive heating. MHD3 is a resistive MHD simulation with at , analogous to the simulations performed by Igumenshchev & Narayan (2002).

### 4.1 Results

Figure 10 shows the mass accretion rate at the inner radius as a function of time. The weak-field simulation S2 shows no significant differences relative to the Bondi solution; the factor of decrease in for S2 in Figure 10 is due to the fact that thermal conduction leads to small changes in the density and temperature profiles of the flow. In addition, there is no indication that the MTI is present, unlike in the simulations at large radii discussed in the previous section. By contrast, all of the simulations with at – both with and without anisotropic conduction – show a significant decrease in after a few free-fall times, with at late times. We focus our interpretation on the simulations with conduction, since our results for simulations without conduction are similar to those presented in Igumenshchev & Narayan (2002). Note that the simulations with anisotropic conduction show episodic “flares” in and other dynamical quantities that are not as prominent in the MHD simulations; these are discussed more below.

Figure 11 shows the overall energetics for the simulation with (S3). A significant fraction of the gravitational potential energy released by the inflowing matter (solid line) is converted into magnetic energy by flux freezing. Resistive dissipation of oppositely directed field lines leads to plasma heating. Although our numerical routine is not conservative, we do have an explicit resistivity and Figure 11 shows that the total resistive heating in the computational domain (dotted line) is similar to the gravitational potential energy released by accretion.2 In turn, much of this energy is transported out of the computational domain by conduction, as shown by the short-dashed line in Figure 11, which is the conductive luminosity at the outer radius . The results for runs S1 & S4 (which have different conductivities ) look similar.

Figure 13 shows how the time and angle averaged density and temperature profiles of the accretion flow change in the presence of resistive heating and conduction. The initial temperature of the flow is much less than the virial temperature because we initialize a Bondi solution with . By contrast, both the MHD simulation and the simulations with anisotropic conduction have time-averaged temperatures of because much of the accretion energy is thermalized as heat via resistive dissipation. The resulting high temperatures lead to a flow in approximate hydrostatic equilibrium, with a mean radial velocity that is highly subsonic at (see Fig. 14). Figure 13 also shows that the density of the flow in the inner region is much lower in the simulations with heating than in the initial solution. This is consistent with accretion rate being (Fig. 10).

In Figures 12-14 we have presented results averaging over the full range of , from to . However, the initially vertical magnetic field lines impose an asymmetry on the flow so that the flow profiles are not spherically symmetric. As seen in Figure 16 (discussed below), the temperature is larger in the equatorial region compared to the pole, by roughly a factor of 10 at small radii; the same is true for the density. The equatorial region is primarily gas pressure dominated while the poles are magnetically dominated. The relative importance of the different energy transport mechanisms (shown in Fig. 12) also varies with angle; in particular, conduction dominates at the poles (and averaged over all angles), but is less important in the equatorial region where the field lines are less ordered. Overall, the angular asymmetry we find in these nonrotating simulations is not as large as in simulations with significant angular momentum (e.g., Hawley & Balbus, 2002) which show “funnels” with radial outflows. In our spherical simulations, the time averaged flow is always inwards at all angles (although there can be outflow at certain times; see Fig. 16).

Although there are modest quantitative differences between the simulations with anisotropic conduction and those without (Figs. 10-14), our results demonstrate that many of the properties of magnetized spherical accretion are in fact rather similar in the two cases. An important difference, however, is that for even modest thermal conductivities (; eq. [7]), conduction is the dominant mechanism of energy transport and transports a significant fraction of the energy released by magnetic dissipation to large radii (Fig. 11 & 12). By contrast, in the MHD simulations without conduction, we find that the net energy transport is inwards. Although convective motions are present, the convective luminosity is much smaller than at all radii.

An additional interesting difference between the simulations with and without conduction is the “flaring” in the simulations with anisotropic conduction (Fig. 10-11). The accretion rate, resistive heating, and conduction luminosity all show regular rapid increases (“flares”) by a factor of on a timescale orbital periods at (see Fig. 15). The characteristic timescale between such flares is a fraction of an orbital period at . The flare timescales do not depend sensitively on the value of the thermal conductivity, as long as . As shown in Fig. 10, flares are most prominent for run S2, the run with the smallest mass accretion rate among simulations (see Table 2).

To assess the robustness of this flaring, we performed several checks. Similar flaring is much less pronounced in the MHD simulations and in simulations with an isotropic thermal conductivity, even though the reduction in relative to is similar; in addition, the flaring is more pronounced for larger values of . To assess whether the flaring is a numerical artifact, we carried out simulations implementing anisotropic conduction using the method of Parrish & Stone (2005) and Sharma et al. (2006), which uses simple finite differencing (our method, based on Sharma & Hammett 2007, limits the heat fluxes so that the temperature is always positive). We find quite similar flaring with the two methods.3 We also changed the magnitude of the artificial viscosity and the explicit resistivity, and the form of the explicit resistivity (e.g., instead of eq. 6, where and are the local isothermal sound speed and Alfvén speed, respectively), but the flaring in run S3 (Fig. 10) persisted. These checks argue against the flaring being a numerical artifact. The physical origin of this flaring is, however, not clear. The MTI could in principle be present because the inflow is subsonic, but it should be stabilized by magnetic tension. The flaring we find is qualitatively similar to the “relaxation oscillations” seen in some 1D time-dependent models of spherical accretion with radiative heating (e.g., Cowie et al., 1978) and could be analogous, with spatially extended heating and redistribution of energy by anisotropic conduction precluding the existence of a steady state solution.

Figure 15 shows a number of quantities characterizing one of the accretion/heating flares over a small interval in time. The gravitational potential energy released by accretion () peaks slightly earlier in time than the volume integrated resistive heating rate (); this is consistent with the field being amplified by flux freezing and then dissipated by resistive diffusion. Figure 15 also shows the maximum values of the local resistive heating rate (), the internal energy per unit volume (), and the temperature, and the minimum value of the density (where the maximum/minimum are over all grid cells in the computational domain). The maximum resistive heating rate and maximum internal energy have a time-dependence similar to that of the accretion energy. The peak temperature and minimum density, on the other hand, occur after the majority of the resistive heating has taken place, and their evolution in time is much more rapid. In fact, the increase in temperature (decrease in density) occurs on a timescale comparable to the free-fall timescale at small radii.

Figure 16 shows contour plots of temperature relative to the virial temperature and local resistive heating rate (in arbitrary units) in the inner at the time of the peak in max[T] in Figure 15 (). Also shown are the velocity unit vectors. The resistive heating is strongly concentrated in the equatorial region at (near the sonic point; see Fig. [14]) where rapid inflow brings together oppositely directed field lines. Strong heating leads to modestly supervirial temperatures at small radii. This energy diffuses to large radii along a finite set of field lines because of the anisotropic nature of conduction. The supervirial temperatures also reverse the inflow along particular field lines, leading to a weak mass outflow (as shown by the velocity unit vectors in Fig. 16).

## 5 Discussion

We have studied the effects of anisotropic thermal conduction on the dynamics of spherical accretion flows using global axisymmetric simulations. Our simulations address two different physical effects that operate on different radial scales in spherical accretion flows.

### 5.1 The MTI in Spherical Accretion Flows

We have calculated the nonlinear saturation of the magnetothermal instability (MTI) at radii in spherical accretion flows (§3), where is the gravitational sphere of influence of the central object (of mass ) and is the isothermal sound speed of the ambient plasma. The MTI does not grow appreciably at radii for the canonical Bondi spherical accretion solution because the inflow is supersonic. Our simulations are the first global simulations of the MTI, but our results are largely consistent with the previous local simulations of Parrish & Stone (2005, 2007). We have focused on studying the saturation of the MTI for initially weak vertical magnetic fields ().

We find that the MTI saturates by aligning the magnetic field lines with the temperature gradient, i.e., by making the field lines primarily radial (Figs. 2 - 4). This occurs on approximately the local dynamical time (). Since the field lines end up largely radial, the rearrangement of the magnetic field results in a radial heat flux similar to the field-free (Spitzer) value; in the saturated state, energy transport by conduction dominates transport by the convective motions associated with the MTI (Fig. 7). Given the significant radial heat flux, the density and temperature profiles of the accretion flow adjust to satisfy , where is the conductive flux of energy. In all cases of that we considered (see §2.1), the mass accretion rate in the presence of the MTI differs only modestly from the canonical Bondi rate (Fig. 6).

Our results indicate that the saturation of the MTI is quasilinear in the weak field limit: the magnetic and convective energy densities in the saturated state are roughly proportional to the initial magnetic energy density (Fig. 8).4 The net result is that the MTI amplifies the magnetic field by a factor of independent of the initial magnetic field strength (for weak fields). For stronger initial fields, magnetic tension can saturate the MTI before the field has become largely radial (Parrish & Stone, 2005, 2007). In this case, the amplification of the magnetic field will be less than we have found here.

The quasilinear saturation of the MTI can be understood qualitatively using linear theory (Quataert, 2008); as a result, we believe that this is a general property of the MTI, rather than being specific to the spherical accretion problem considered in this paper. The linear growth rate of the MTI is reduced by a factor of for nearly radial fields, and the growth occurs on progressively smaller scales for small (eq. [8]). Thus, as the field lines become radial, the MTI effectively shuts off its own growth. This is in stark contrast to the magnetorotational instability in differentially rotating plasmas, which is an exact nonlinear solution of incompressible MHD, and for which growth continues unimpeded even when (Goodman & Xu, 1994).

One simplification in the current calculations is that we have not included anisotropic viscosity; the diffusion coefficient for anisotropic viscosity in a collisional plasma is smaller than the diffusion coefficient for anisotropic conduction by a factor of ; in a collisionless plasma, the relative importance of these two effects is more complex and depends on small-scale kinetic microinstabilities (e.g., Sharma et al. 2007). Local simulations of the MRI find that the saturation of the instability is sensitive to the ratio of the resistivity to the viscosity (e.g., Fromang et al. 2007). In future calculations, the effects of anisotropic viscosity on the evolution of the MTI should thus be studied, but this is probably first best done in a local calculation, rather than in the global geometry considered here.

#### Application to the Galactic Center

One potential application of the MTI at large radii in spherical accretion flows is to the observed Faraday rotation from the Galactic center. Millimeter observations of Sgr A* find a rotation measure (RM) of rad m that has been roughly constant over the past years (Aitken et al., 2000; Bower et al., 2005; Marrone et al., 2007). Propagation through the interstellar medium cannot account for this large RM. In addition, X-ray observations of the thermal plasma in the central parsec of the Galactic center find that the electron number density is cm and the temperature is keV at a distance of from the black hole (Baganoff et al., 2003; Quataert, 2004). If the magnetic field on this scale were uniform, radial, and in rough equipartition with the thermal pressure, mG and rad m, two orders of magnitude smaller than what is observed. This shows that the observed RM must be produced at distances from Sgr A*.

Our calculations show that given a sufficiently coherent seed magnetic field at large radii, the MTI will generate a relatively coherent radial field at distances of from Sgr A*. Given the measured density and temperature of the plasma, the electron mean free path is at radii and thus the MTI grows on a dynamical time. Numerical simulations of accretion onto Sgr A* from stellar winds located at find that the circularization radius of the flow is (Cuadra et al., 2006). If the Bondi solution is applicable from the scale of the observed plasma at to , then at fixed and the RM produced by the roughly spherical inflowing plasma at large radii around Sgr A* is

It is difficult to determine with confidence whether the stellar winds feeding Sgr A* – or the shocks between such winds – produce a sufficiently large and coherent ’seed’ magnetic field at to account for the RM towards Sgr A*. However, for an initial value of of (assumed ) at , the combined action of the MTI and flux freezing will decrease to at . Thus even a relatively weak magnetic field (; ) on scales of pc will be amplified by the MTI and flux freezing to a value capable of explaining the observed RM from Sgr A*.5 Our explanation for the observed RM also naturally accounts for its stability, since the dynamical time at is years. By contrast, if the RM is produced in the rotating accretion flow even closer to the black hole, it is harder to account for the observed stability (Sharma, Quataert, & Stone, 2007).

### 5.2 The Effects of Thermal Conduction on the Dynamics of Spherical Accretion

The second set of simulations we have carried out address the effects of anisotropic conduction and resistive dissipation of magnetic energy on the dynamics of spherical accretion flows at small radii (§4). Our results are in reasonable agreement with those of Igumenshchev & Narayan (2002), who studied the same physical problem but without anisotropic conduction. Although the MTI does not dramatically modify the dynamics of spherical accretion at , it does generate a significant radial magnetic field that can be further amplified by flux freezing. For , supersonic inflow in (magneto)hydrodynamic Bondi accretion leads to magnetic field amplification by flux freezing, converging field lines, and resistive heating of the plasma. When the magnetic energy density is comparable to the gravitational energy density, resistive heating can heat the plasma up to the virial temperature and significantly modify the dynamics of spherical accretion (Igumenshchev & Narayan, 2002). This is only likely to occur for sub-Eddington accretion rates onto compact objects, since only under these conditions is radiative cooling negligible.

All of our calculations with resistivity and strong magnetic fields at – both with and without anisotropic conduction – show a significant decrease in after a few free-fall times, with in the quasi steady state at late times (Fig. 10). The high temperatures generated by resistive heating stifle the inflow, reducing the accretion rate (Fig. 10) and generating a subsonic (Fig. 14) nearly hydrostatic atmosphere around the central object. For realistic values of the thermal conductivity for collisionless plasmas (see §2.1), we find that conduction transports the majority of the energy released by accretion to large radii (Figs. 11 & 12). By contrast, although turbulent motions are present, there is no net convective energy flux from small to large radii, even in our simulations without conduction (Fig. 12).6 This indicates that significant radial transport of energy (be it conductive or convective) is not necessary to suppress accretion onto the central object; instead, heating that raises the temperature to the virial temperature appears sufficient. Indeed, it is well-known that there is no supersonic Bondi solution for adiabatic indices . Heating via resistive dissipation and/or artificial/numerical viscosity effectively increases the adiabatic index above this value, thus stifling accretion onto the central object. Thermal conduction can further reduce the accretion rate by conducting the dissipated energy to larger radii.

Overall, our results are in reasonable agreement with Johnson & Quataert (2007)’s simplified one-dimensional models for the structure of spherical accretion in the presence of conduction and plasma heating. Our simulations confirm their hypothesis that a significant fraction of the energy released by accretion is transported to large radii by thermal conduction. In addition, their prediction that should decrease by orders of magnitude – with a factor of suppression due to heating alone (their Fig. 8) – is consistent with our numerical results.

The reduction in that we find is similar to that needed to explain the low luminosity of Sgr A* (e.g., Sharma, Quataert, & Stone 2007). It is also reasonably consistent with the limits on from Faraday rotation measurements (e.g., Marrone et al. 2007). However, it is unlikely that the present calculations are quantitatively applicable to Sgr A* given the estimated circularization radius of the accretion flow in Sgr A* of (Cuadra et al., 2006). The gravitational potential energy released accreting to is unlikely to significantly modify the accretion flow dynamics. Thus the accretion rate and outward conductive energy flux will be largely determined in the rotationally supported flow at radii . Nonetheless, our results demonstrate that even in the absence of angular momentum, the Bondi accretion solution does not apply to dilute magnetized accretion flows. Instead, the accretion flow is heated to nearly the virial temperature by resistive dissipation of magnetic energy, the accretion rate is , and most of the gravitational potential energy released by accretion is transported to large radii by thermal conduction. In addition, in our simulations with anisotropic conduction, the accretion and heating are time-dependent, with episodic “flaring” (Figs. 10-11) followed by mass outflow along a subset of the field lines (Fig. 16). It is important to determine if analogous effects are present in rotating accretion flows (which is by no means clear), given the possible application to the broad-band electromagnetic flares observed from Sgr A* (e.g., Marrone et al., 2008) and other systems.

## Acknowledgments

We thank Bryan Johnson and Ian Parrish for useful discussions. We thank the referee Steve Balbus for detailed comments. PS and EQ were supported in part by NASA grant NNG06GI68G and the David and Lucile Packard Foundation. PS was partly supported by DOE award DE-FC02-06ER41453 to Jon Arons. JS was supported by DOE grant DE-FG52-06NA26217. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Computational resources were also provided by the Princeton Plasma Physics Laboratory Scientific Computing Cluster. Part of the code was tested with the Tungsten cluster at the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign.

### Footnotes

1. For the sonic point lies very close to the origin. It is, however, numerically desirable to have the sonic point within the computational domain so that the inner boundary conditions do not influence the solution. We have verified that we get similar results for a run with , but with rest of the parameters same as the fiducial run R1.
2. In addition to explicit resistivity – which is the dominant source of heating in our calculations – we also include artificial viscosity for shock capturing. Explicit resistivity, artificial viscosity, and inherent numerical diffusion all lead to increase in entropy in the simulations.
3. At late times in the simulations with conduction implemented using simple finite differencing, the temperature can become negative, which demonstrates the importance of using a method that guarantees that heat flows from hot to cold, such as that of Sharma & Hammett (2007).
4. Ian Parrish has found this same result in local three-dimensional Cartesian simulations (private communication).
5. If most O stars have surface magnetic fields comparable to the kG field measured in the young O star Orinis C (Donati et al., 2002), stellar winds alone will fill the central pc with a magnetic field mG.
6. Our interpretation of the MHD simulations of spherical accretion without conduction thus differs from that of Igumenshchev & Narayan (2002), who emphasized the importance of radial energy transport by convection. This needs to be investigated further.

### References

1. Aitken, D. K. et al. 2000, ApJ, 534, L173
2. Baganoff, F. K. et al. 2003, ApJ, 591, 891
3. Balbus, S. A. & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1
4. Balbus, S. A. 2000, ApJ, 534, 420
5. Bondi, H., 1952, MNRAS, 112, 195
6. Bower, G. C. et. al. 2005, ApJ, 618, 29
7. Braginskii, S. I. 1965, Reviews of Plasma Physics, v. 1, p. 205, Consultants Bureau
8. Chandran, B. D. G. & Dennis, T. J. 2006, ApJ, 642, 140
9. Chandran, B. D. G. & Rasera, Y. 2007, arXiv:0707.2351
10. Cowie, L. L. & McKee, C. F. 1977, ApJ, 211, 135
11. Cowie, L. L., Ostriker, J. P., & Stark, A. A. 1978, ApJ, 226, 1041
12. Cuadra, J., Nayakshin, S., Springel, V., & di Matteo, T. 2006, MNRAS, 366, 358
13. Donati, J.-F., Babel, J., Harries, T. J., Howarth, I. D., Petit, P., & Semel, M. 2002, MNRAS, 333, 55
14. Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464
15. Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123
16. Goodman, J., & Xu, G. 1994, ApJ, 432, 213
17. Hawley, J. F. & Balbus, S. A. 2002, ApJ, 573, 738
18. Igumenshchev, I. V. & Narayan, R. 2002, ApJ, 566, 137
19. Johnson, B. M. & Quataert, E. 2007, ApJ, 660, 1273
20. Marrone, D. P. et al. 2007, ApJ, 654, L57
21. Marrone, D. P. et al. 2008, ApJsubmitted, arXiv:0712.2877
22. Paczynski, B. & Wiita, J. 1980, A&A, 88, 23
23. Parrish, I. J. & Stone, J. M. 2005, ApJ, 633, 334
24. Parrish, I. J. & Stone, J. M. 2007, ApJ, 664, 135
25. Parrish, I. J. & Quataert, E. 2008, ApJ, in press
26. Piffaretti, R., Jetzer, P., Kaastra, J. S., & Tamura, T. 2005, A&A, 433, 101
27. Quataert, E. 2003, Astron. Nachr., 324, 435
28. Quataert, E. 2004, ApJ, 613, 322
29. Quataert, E. 2008, ApJ, 673, 758
30. Rees, M. J. et al. 1982, Nature, 295, 17
31. Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1
32. Shvartsman, V. F., 1971, Soviet Astron., 15, 377
33. Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. M., 2006, ApJ, 637, 952
34. Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. M. 2007, ApJ, 667, 714
35. Sharma, P., Quataert, E., & Stone, J. M. 2007, ApJ, 671, 1696
36. Sharma, P. & Hammett, G. W. 2007, J. Comp. Phys., 227, 123
37. Snyder, P. B., Hammett, G. W., & Dorland, W. 1997, Phys. Plasmas, 4, 3974
38. Spitzer, L. 1962, Physics of Fully Ionized Gases (New York: Wiley Interscience)
39. Stone, J. M., Norman, M. L. 1992a, ApJS, 80, 753
40. Stone, J. M., Norman, M. L. 1992b, ApJS, 80, 819
41. Stone, J. M. & Pringle J. E. 2001, MNRAS, 322, 461
42. Tanaka, T. & Menou, K. 2006, ApJ, 649, 345
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters