Sgr A* Modeling by 3D GRMHD and Polarized Transfer

Sagittarius A* Accretion Flow and Black Hole Parameters from General Relativistic Dynamical and Polarized Radiative Modeling


We obtain estimates of Sgr A* accretion flow and black hole parameters by fitting polarized sub-mm observations with spectra computed using three-dimensional (3D) general relativistic (GR) magnetohydrodynamical (MHD) (GRMHD) simulations. Observations are compiled from averages over many epochs from reports in papers for estimating the mean fluxes , linear polarization (LP) fractions, circular polarization (CP) fractions, and electric vector position angles (EVPAs). GRMHD simulations are computed with dimensionless spins over a time interval. We perform fully self-consistent GR polarized radiative transfer using our new code to explore the effects of spin , inclination angle , position angle (PA), accretion rate , and electron temperature ( is reported for radius ). By fitting the mean sub-mm fluxes and LP/CP fractions, we obtain estimates for these model parameters and determine the physical effects that could produce polarization signatures. Our best bet model has , , , , and  K at . The sub-mm CP is mainly produced by Faraday conversion as modified by Faraday rotation, and the emission region size at  GHz is consistent with the VLBI size of as. Across all spins, model parameters are in the ranges , , and K. Polarization is found both to help differentiate models and to introduce new observational constraints on the effects of the magnetic field that might not be fit by accretion models so-far considered.

Subject headings:
accretion, accretion disks –– black hole physics –– Galaxy: center –– radiative transfer –– relativistic processes — polarization

1. Introduction

The mass of the Galactic Center black hole (BH) is (Ghez et al., 2008; Reid et al., 2008; Gillessen et al., 2009) and the spin is uncertain (Huang et al., 2009b; Broderick et al., 2009; Moscibrodzka et al., 2009; Broderick et al., 2010; Dexter et al., 2010). It resides at a distance  kpc. Because of its proximity, it has been observed in many wavelengths: -rays, X-rays, IR, (sub-)mm, and radio. X-ray bremsstrahlung emission originates from hot gas at large radii where the BH’s gravity becomes important (Narayan, Yi & Mahadevan, 1995; Narayan et al., 1998; Shcherbakov & Baganoff, 2010) and Compton-scattered emission originates from close to the horizon (Moscibrodzka et al., 2009). X-rays at large radii are spatially resolved and have been used to constrain dynamical models for this region (Shcherbakov & Baganoff, 2010). The sub-mm emission is cyclo-synchrotron emission originating from close to the BH. Cyclo-synchrotron emission is polarized, and both linear and circular polarizations have been observed from Sgr A* at several sub-mm wavelengths. The accretion flow was recently resolved at  GHz (Doeleman et al., 2008; Fish et al., 2011). General relativistic (GR) effects were deemed necessary to explain the small size with full width at half maximum (FWHM) of as. Radio emission is also produced by cyclo-synchrotron at larger distances from the BH. Relativistic frame-dragging is important near the BH, so sub-mm polarized observations and the Compton-scattered X-rays might help to constrain the BH spin. The goal of the present paper is to model the sub-mm in the range of  GHz to  GHz in order to estimate the accretion flow and black hole parameters.

Sgr A* is a variable source with a variability amplitude routinely reaching in sub-mm. A popular approach is to fit simultaneous observations (e.g. Yuan, Quataert & Narayan 2004; Broderick et al. 2009), in particular, the set from Falcke et al. (1998). However, in such an approach, one would use a single simultaneous set of observations. However, simultaneous observations of fluxes, linear polarization (LP), and circular polarization (CP) fractions at several frequencies are not available. So we consider non-simultaneous statistics of all observations at all frequencies and find the mean values and standard errors of quantities at each frequency.

Numerous accretion flow models have been applied to the Galactic Center: advection-dominated accretion flow (ADAF) (Narayan & Yi, 1995), advection-dominated inflow-outflow solution (ADIOS) (Blandford & Begelman, 1999), jet-ADAF (Yuan, Markoff & Falcke, 2002), jet (Maitra, Markoff & Falcke, 2009), and viscous and magnetohydrodynamical (MHD) numerical simulations. The quasi-analytical models are useful because there is little expense in changing parameters. However, they have a large number of free parameters and also incorporate many assumptions that are not justifiable from first principles (Huang et al., 2008, 2009a), which leads to systematic uncertainties in all fits. Numerical simulations require fewer inputs and are useful for more quantitative modeling of the plasma near a rotating BH. General relativistic (GR) MHD (GRMHD) simulations (especially three-dimensional (3D) simulations), which are run over a sufficiently long duration, are still computationally expensive and involve state-of-the-art codes that are still being developed (McKinney & Blandford, 2009; Fragile et al., 2009; Noble & Krolik, 2009; Moscibrodzka et al., 2009; Penna et al., 2010). Yet, these expensive 3D simulations are required to model the turbulent disk flow, because 2D axisymmetric simulations cannot sustain turbulence as shown by generalizations of Cowling’s anti-dynamo theorem (Hide & Palmer, 1982). Given their expense, such 3D GRMHD simulations are limited to a region relatively close to the BH (Dexter et al., 2009; Moscibrodzka et al., 2009), whereas some emission and some Faraday rotation might happen far from the BH. So we analytically extend the modeled region out to , perform radiative transfer, and find the best fit to the data. The extension to large radius allows us to define the electron temperature more consistently (Sharma et al., 2007). We find a posteriori (see Appendix A) that the simulated polarized spectra are not overly sensitive to the details of the analytic extensions of density and temperature, but may depend on the extension of the magnetic field.

The radiation close to the BH has been modeled in Newtonian (Yuan, Quataert & Narayan, 2004) and quasi-Newtonian approximations (Goldston, Quataert & Igumenshchev, 2005; Chan et al., 2009). It has been modeled in GR assuming unpolarized (Fuerst & Wu, 2004; Dexter et al., 2009; Dolence et al., 2009) and polarized (Broderick et al., 2009; Shcherbakov & Huang, 2011) light. Fitting the total flux spectrum might not be sufficient to estimate the spin, and naturally one expects polarization to provide extra observational constraints. Spin values from (Broderick et al., 2009) to (Moscibrodzka et al., 2009) have been estimated. We neglect Comptonization (Moscibrodzka et al., 2009) and radiation from non-thermal electrons (Mahadevan, 1998; Özel, Psaltis & Narayan, 2000; Yuan, Quataert & Narayan, 2004). Emissivities are calculated in the synchrotron approximation (Legg & Westfold, 1968; Sazonov, 1969; Pacholczyk, 1970; Melrose, 1971) with an exact thermal electron distribution. Discrepancies with the exact cyclo-synchrotron emissivities (Leung, Gammie & Noble, 2011; Shcherbakov & Huang, 2011) are negligible as estimated in § 5. Exact Faraday rotation and conversion expressions are used (Shcherbakov, 2008).

We compare simulated spectra to observed ones at many frequencies simultaneously, extending an approach pioneered by Broderick et al. (2009) and Dexter et al. (2009). We compute the average observed spectra, find the deviations of the means, and then compare them to the average simulated spectra. In the search for the best fit models, we are guided by the value of , which is the normalized sum of squares of normalized deviations. Yet, we leave the exploration of the statistical meaning of to future work. We search the space of all parameters: spin , inclination , ratio of proton to electron temperatures ( is reported for radius ), and accretion rate to find the minimum models.

We summarize the radio/sub-mm observations of Sgr A* in § 2. Our 3D GRMHD simulations are described in § 3 together with the physically-motivated extension to large radii, and the electron heating prescription. We run simulations for dimensionless spins . The GR polarized radiative transfer technique is described in § 5.

The set of observations we consider consists of the spectral energy distribution (SED) within the  GHz to  GHz frequency range, linear polarization (LP) fractions at  GHz,  GHz, and  GHz, and circular polarization (CP) fractions at  GHz and  GHz. In § 6 we discuss our results: the best fit models to the observations, the importance of various physical effects in producing the observed CP and LP and electric vector position angle (EVPA), and image size estimates. We produce the simulated images of total and polarized intensities. Discussion in § 7 compares the results to previous estimates, emphasizes the significance of polarization, notes the sources of error, and outlines prospects for future work. In Appendix A we describe a number of convergence tests of our GR polarized radiative transfer code and the radial extension of the dynamical model. Throughout the paper we measure distance and time in the units of BH mass by setting the speed of light and gravitational constant to unity.

2. Observations

Sgr A* is known to be a highly variable source, yet quiescent models of Sgr A* emission are popular and useful. Unlike the drastic variations of X-ray and NIR fluxes (Baganoff et al., 2001; Genzel et al., 2003), sub-mm fluxes do not vary by more than a factor of (Zhao et al., 2003). We compile the set of observed polarized fluxes at each frequency, then we find the mean spectrum and the errors of the mean fluxes.

Previously, the observed flux spectra were compiled by Yuan, Quataert & Narayan (2004); Broderick et al. (2009). However, both papers summarize a limited set of observations and concentrate on simultaneously observed fluxes. Sub-mm flux data reported in Yuan, Quataert & Narayan (2004) consist of a short set of observations by Falcke et al. (1998) and one set of SMA observations by Zhao et al. (2003). Broderick et al. (2009) adds to these the rest of SMA total flux data (Marrone et al., 2006a, b, 2007, 2008). So out of at least papers on sub-mm observations of Sgr A* were taken into account. We compute an averaged spectrum based on papers reporting sub-mm observations of Sgr A*.

[GHz] Telescopes [Jy] LP [] CP [] EVPA []
8.45 VLA (Serabyn et al., 1997; Falcke et al., 1998; Bower et al., 1999a; An et al., 2005) 4 (Bower et al., 1999a)
14.90 VLBA, VLA (Serabyn et al., 1997; Falcke et al., 1998; Bower et al., 2002; Herrnstein et al., 2004; An et al., 2005; Yusef-Zadeh et al., 2009) 5 (Bower et al., 2002)
22.50 VLBA, VLA (Serabyn et al., 1997; Falcke et al., 1998; Bower et al., 1999b; Herrnstein et al., 2004; An et al., 2005; Lu et al., 2008; Yusef-Zadeh et al., 2007, 2009) 6 (Bower et al., 1999b; Yusef-Zadeh et al., 2007)
43 GMVA, VLBA, VLA (Falcke et al., 1998; Lo et al., 1998; Bower et al., 1999b; Herrnstein et al., 2004; An et al., 2005; Shen et al., 2005; Krichbaum et al., 2006; Yusef-Zadeh et al., 2007; Lu et al., 2008; Yusef-Zadeh et al., 2009) 7 (Bower et al., 1999b; Yusef-Zadeh et al., 2007)
88 BIMA, MPIfR, VLBA, VLA, Nobeyama, NMA, CARMA (Falcke et al., 1998; Krichbaum et al., 1998; Bower et al., 1999b; Doeleman et al., 2001; Miyazaki et al., 2004; Shen et al., 2005; Krichbaum et al., 2006; Macquart et al., 2006; Lu et al., 2008; Yusef-Zadeh et al., 2009) 8 (Bower et al., 1999b; Macquart et al., 2006) -49 (Bower et al., 1999b; Shen et al., 2005; Macquart et al., 2006)
102 OVRO, CSO-JCMT, Nobeyama, NMA, IRAM (Serabyn et al., 1997; Falcke et al., 1998; Miyazaki et al., 2004; Mauerhan et al., 2005; Yusef-Zadeh et al., 2009)
145 Nobeyama, NMA, IRAM, JCMT (Falcke et al., 1998; Aitken et al., 2000; Miyazaki et al., 2004; Yusef-Zadeh et al., 2009)
230 IRAM, JCMT, BIMA, SMA, OVRO (Serabyn et al., 1997; Falcke et al., 1998; Aitken et al., 2000; Bower et al., 2003, 2005; Zhao et al., 2003; Krichbaum et al., 2006; Marrone et al., 2006a, 2007, 2008; Doeleman et al., 2008; Yusef-Zadeh et al., 2009) (Bower et al., 2003, 2005; Marrone et al., 2007, 2008) 10 (Munoz et al. (2009, 2011)) (Bower et al., 2003, 2005; Marrone et al., 2007, 2008)
349 SMA, CSO, JCMT (Aitken et al., 2000; An et al., 2005; Marrone et al., 2006b, 2007, 2008; Yusef-Zadeh et al., 2009) (Marrone et al., 2006b, 2007) 11 (Munoz et al. (2011)) (Marrone et al., 2006b, 2007)
674 CSO, SMA (Marrone et al., 2006a, 2008; Yusef-Zadeh et al., 2009)
857 CSO (Serabyn et al., 1997; Marrone et al., 2008; Yusef-Zadeh et al., 2009)
Table 1Summary of Sgr A* radio/sub-mm observations

The reported observations vary in covered period from several hours (An et al., 2005) to several years (Zhao et al., 2003; Krichbaum et al., 2006). We know that variations of a factor of may happen within several hours (Yusef-Zadeh et al., 2009), whereas variations by more than a factor of several are never observed in the sub-mm. So, fluxes observed more than a day apart are weakly correlated. The issue of autocorrelation in timescales will be addressed in future work. We consider the following averaging technique to sample the distributions of fluxes. First, we define groups of close frequencies, the frequencies in each group being different by no more than several percent from the mean. There are such groups (see Table 1). We exclude papers reporting single frequencies far from the mean of each group. In particular, the  GHz and  GHz observations of Li et al. (2008); Falcke et al. (1998) and the  GHz observations of Bower et al. (2001) are excluded. A mean frequency is ascribed to represent each group. Then, we take all reported observations of each polarization type (total flux, LP and CP fraction, EVPA) for each group and draw the largest sample of fluxes/polarization fractions, taking observations separated by at least hours. When several fluxes are reported over a period of several hours (Yusef-Zadeh et al., 2009), we draw one data point from the very beginning of the observation, unless a flare is reported to occur at that time. Some of the published observations have large error bars. Often such data are produced by observing in sub-mm with large beam size, but light from Sgr A* is blended with dust and other sources. In particular, SMT data (Yusef-Zadeh et al., 2009), early CSO measurements (Serabyn et al., 1997), and early JCMT measurements (Aitken et al., 2000) may have such issues, so we exclude these data from the sample. The interferometric observations, especially with VLBI, help to reduce the error from otherwise unreliable observations, e.g. with BIMA array (Bower et al., 2001). However, some inconsistencies still exist for simultaneous observations at the same frequency with different instruments (Yusef-Zadeh et al., 2009).

After the sample of fluxes, polarization fractions, and EVPAs are found for each frequency group, we compute the mean and the standard error. The summary of results is presented in Table 1. CP fractions of at  GHz and at  GHz are based on SMA observations by Munoz et al. (2011) with the reported instrumental error. Note that standard errors in our total flux samples are smaller than the error bars of prior observations (Falcke et al., 1998; Yuan, Quataert & Narayan, 2004; Broderick et al., 2009), but still larger compared to contemporary single-observation instrumental errors (Marrone et al., 2007). That is, we do not incorporate instrumental error in our estimates of standard error of the mean fluxes or and at  GHz and  GHz (even though the instrumental error of at  GHz is large). We do not incorporate the source size measurements (Doeleman et al., 2008) in calculating , but we check that the best bet model is not inconsistent with those observations. Figure 1 shows a compilation of the mean quantities and their Gaussian standard errors. The data are represented by both error bars and the interpolated shaded area. A red dashed curve on the plot represents the analytic approximation , where flux is in Jy and frequency is in GHz.

Figure 1.— Mean observed SEDs of specific flux , linear polarization (LP) fraction, electric vector position angle (EVPA), and circular polarization (CP) fraction. The error bars show standard error of the mean. The dashed line on the plot represents the analytic approximation for frequency in GHz (not the simulated SED). As noted in Table 1, the error is instrumental for CP at high frequencies and LP at  GHz, whereas it is computed from a sample of observed quantities for flux, EVPA at all frequencies, and LP at high frequencies.

3. Three-Dimensional Grmhd Simulations

Our radiative transfer calculations take the results of simulations of accretion flows onto BHs as input. These simulations are similar to those in Penna et al. (2010). Below, we review the methodology.

3.1. Governing Equations

We simulate radiatively inefficient accretion flows (RIAFs) onto rotating BHs using a three-dimensional fully general relativistic code (see § 3.3). The BH is described by the Kerr metric. We work with Heaviside-Lorentz units. Our five simulations correspond to different choices of the dimensionless BH spin parameter: , and . The self-gravity of the RIAF is ignored.

The RIAF is a magnetized fluid, so we solve the GRMHD equations of motion (Gammie et al., 2003). Mass conservation gives:


where is the fluid frame rest-mass density, is the contravariant 4-velocity, and is the covariant derivative. Energy-momentum conservation gives


where the stress energy tensor includes both matter and electromagnetic terms,


where is the internal energy density and is the ideal gas pressure with . Models with show minor differences compared to models with (McKinney & Gammie, 2004; Mignone & McKinney, 2007). The contravariant fluid-frame magnetic 4-field is given by and is related to the lab-frame 3-field via , where is a projection tensor and is the Kronecker delta function (Gammie et al., 2003). We often employ below, which is the orthonormal magnetic field vector in a comoving locally flat reference frame (Penna et al., 2010). The magnetic energy density () and magnetic pressure () are then given by . Note that the angular velocity of the gas is .

Magnetic flux conservation is given by the induction equation


where , and is the determinant of the metric. No explicit resistivity or viscosity is included, but we use a shock-capturing Godunov method that fully conserves energy. So, all dissipation from shocks and numerical diffusivity (e.g. in shear flows or current sheets) is fully captured, as required to study RIAFs.

In Penna et al. (2010), we studied both RIAFs and geometrically thin radiatively efficient disks. For the later case, a cooling term was added to the energy-momentum equation (2) to describe radiative losses and keep the disk thin. The current set of models are all RIAFs, so no cooling term is used. Entropy generated by viscous or resistive dissipation is advected along with the inflow or transported out via convection or in a wind.

3.2. Physical Models

The initial mass distribution is an isentropic equilibrium torus (Chakrabarti, 1985a, b; De Villiers, Hawley & Krolik, 2003) with pressure for . The torus inner edge is at and the maximum density and pressure are at . We initialize the solution so that at the pressure maximum. As in Chakrabarti (1985a), the angular velocity distribution of the initial torus is a power law, where for the Chakrabarti (1985a) -parameter we choose (At large radii .). The thickness of the torus at the pressure maximum is then , where


where is an area element in the plane, and the integral over is a time average over the period when the disk is in a steady state (see §3.6). A tenuous atmosphere fills the space outside the torus. It has the same polytropic equation of state as the torus, , with , and an initial rest-mass density of , corresponding to a Bondi-like atmosphere. The torus is threaded with three loops of weak, poloidal magnetic field: the initial gas-to-magnetic pressure ratio is , where and are the maximum values of the gas and magnetic pressure in the torus. This approach to normalizing the initial field is used in many other studies (Gammie et al., 2003; McKinney & Gammie, 2004; McKinney, 2006a; McKinney & Narayan, 2007b; Komissarov & McKinney, 2007; Penna et al., 2010).

Recent GRMHD simulations of thick disks indicate that the results for the disk (but not the wind-jet, which for us is less important) are roughly independent of the initial field geometry (McKinney & Narayan 2007a, b; Beckwith et al. 2008a, but see also McKinney et al. 2012). The magnetic vector potential we use is given by


with all other initially zero. This is the same as used in Penna et al. (2010). We use , and set if either or . Here is the maximum value of the internal energy density in the torus. We choose and , which gives initial poloidal loops that are roughly isotropic such that they have roughly 1:1 aspect ratio in the poloidal plane. The form of the potential in equation (6) ensures that each additional field loop bundle has opposite polarity. Perturbations are introduced to excite the magneto-rotational instability (MRI). The second term on the right-hand-side (RHS) of equation 6 is a random perturbation: is a random real number generator for the domain to . Random perturbations are introduced in the initial internal energy density in the same way, with an amplitude of . In Penna et al. (2010), it was found that similar simulations with perturbations of and became turbulent at about the same time, the magnetic field energy at that time was negligibly different, and there was no evidence for significant differences in any quantities during inflow equilibrium.

3.3. Numerical Methods

We perform simulations using a fully 3D version of HARM that uses a conservative shock-capturing Godunov scheme (Gammie et al., 2003; Shafee et al., 2008; McKinney, 2006b; Noble et al., 2006; Mignone & McKinney, 2007; Tchekhovskoy, McKinney & Narayan, 2007; McKinney & Blandford, 2009). We use horizon-penetrating Kerr-Schild coordinates for the Kerr metric (Gammie et al., 2003; McKinney & Gammie, 2004), which avoids any issues with the coordinate singularity in Boyer-Lindquist coordinates. The code uses uniform internal coordinates mapped to the physical coordinates . The radial grid mapping is


which spans from to , where is the radius of the outer event horizon. This just ensures the grid never extends inside the inner horizon, in which case the equations of motion would no longer be hyperbolic. The parameter controls the resolution near the horizon. For the outer radial boundary of the box, absorbing (outflow, no inflow allowed) boundary conditions are used.

The -grid mapping is


where ranges from to (i.e. no cut-out at the poles) and is chosen to concentrate grid zones toward the equator. Reflecting boundary conditions are used at the polar axes. The -grid mapping is given by , such that varies from to for a box with . Periodic boundary conditions are used in the -direction. Penna et al. (2010) considered various for thin disks and found little difference in the results. In all of their tests, and we remain above this limit as well. In what follows, spatial integrals are renormalized to refer to the full range in , even if our computational box size is limited in the -direction. For the purpose of radiative transfer, we combine two identical regions of size preserving the orientation to obtain the span of full .

3.4. Resolution and Spatial Convergence

The resolution of the simulations is . This is the fiducial resolution of Penna et al. (2010). Shafee et al. (2008) found this resolution to be sufficient to obtain convergence compared to a similar model. In the vertical direction, we have about 7 grid cells per density scale height. Turbulence is powered by the MRI, which is seeded by the vertical component of the magnetic field (Balbus & Hawley, 1998). The characteristic length scale of the MRI is the wavelength of the fastest growing mode:


where is the vertical component of the Alfvén speed. We find that the MRI is well-resolved in the midplane of disk both initially and in the saturated state.

Penna et al. (2010) studied convergence in , , and and found that models with or , or , and or behaved similar for disks with similar resolution across the disk. Our resolution of the MRI and prior convergence testing by Penna et al. (2010) for similarly-resolved disks justify our choice of grid resolution. It is currently not computationally feasible to perform a similar spin parameter study at much higher resolutions, and future studies will continue to explore whether such simulations are fully converged (Hawley et al., 2011; McKinney et al., 2012).

3.5. Ceiling Constraints

During the simulation, the rest-mass density and internal energy densities can become low beyond the corona, but the code remains accurate and stable for a finite value of , , and for any given resolution. We enforce , , and by injecting a sufficient amount of mass or internal energy into a fixed zero angular momentum observer (ZAMO) frame with 4-velocity , where is the lapse.

We have checked that the ceilings are rarely activated in the regions of interest of the flow. Figure 2 shows the constrained ratios, , , and , as a function of at six radii (, and ) for the model. The data has been time-averaged over the steady state period from to . The ceiling constraints are shown as dashed red lines. The solution stays well away from the ceilings. Thus, the ceilings are sufficiently high.

Figure 2.— Ratios of , , and versus . Black curves correspond to different radii in the flow; from top to bottom, , and . The data is time-averaged over the steady state period of the flow, from to . Numerical ceilings constrain the solution to lie below the dashed red lines, but we see that the solution does not approach these limits.

3.6. Approach to Steady State

We run the simulations from to . The accretion rate, the height- and averaged plasma , and other disk parameters, fluctuate turbulently about their mean values. The simulation reaches a quasi-steady state, when the mean parameter value are time-independent. Figure 3 shows the accretion rate and height- and averaged at the event horizon as a function of time for all five models. We take the period from to to define steady state.

As shown in Penna et al. (2010), for disk models like the one considered, the disk outside the innermost stable circular orbit (ISCO) behaves like the -disk model with across disk thicknesses of . This allows one to accurately infer the timescale for reaching “inflow equilibrium,” corresponding to a quasi-steady flow across all quantities, at a given radius. For by - (the simulation runs till , but the initial are transients not necessarily associated with achieving inflow equilibrium for a simple viscous disk), we use the results in Appendix B of Penna et al. (2010) and find that inflow equilibrium is achieved within a radius of - for models with and for models with . Even for a doubling of the viscous timescale, inflow equilibrium is achieved by - depending upon the BH spin. This motivates using an analytical extension of the simulation solution for radii beyond as described later in § 4.1.

Figure 3.— Accretion rate and height- and averaged versus time at the event horizon for all five models: (solid light cyan), (solid dark red), (long-dashed green), (short-dashed brown), and (dot-dashed orange).

3.7. Evolved Disk Structure

Figure 4 shows matter stream lines as vectors and number density as greyscale map. The large scale vortices existing on a single time shot (panel a) almost disappear when averaged over the duration (panel b) from times to . The density is highest in the equatorial plane on average, but deviations are present on the instantaneous map. The ISCO does not have any special significance: density and internal energy density increase through ISCO towards the BH horizon.

Figure 5 shows magnetic field lines as vectors and comoving electromagnetic energy density as a greyscale map. The structure of magnetic field at early times remembers the initial multi-loop field geometry (Penna et al., 2010), but switches at late times to a helical magnetic field structure resembling a split-monopole in meridional projection. Such switching of magnetic field structure suggests that the final helix with projected split-monopole is a natural outcome of any vertical flux being dragged into the BH (although the amount of magnetic flux threading the hole and disk may be chosen by initial conditions as described in McKinney et al. 2012). The magnetic field structure of a single snapshot (panel a) looks similar to the structure of the linear average between and (panel b). The polar region of the flow has the strongest magnetic field. The magnetic field lines on Figure 5 illustrate only the direction of the field’s poloidal component. The toroidal magnetic field is stronger above and below the midplane of the disk outside of ISCO. The toroidal field strength is comparable to the poloidal field strength inside the ISCO and near the disk midplane.

Figure 4.— Stream lines of velocity (red vectors) and number density (greyscale map) for spin at in the meridional plane( as cylindrical radius): single time snapshot at on the upper (a) panel and time average between and on the lower (b) panel. The corresponding calibration bars of are shown on the right. Number density is normalized by its maximum value, and the vectors show the poloidal velocity direction.
Figure 5.— Magnetic field lines (red vectors) and comoving electromagnetic energy density (greyscale map) for spin at in the meridional plane( as cylindrical radius): single time snapshot at on the upper (a) panel and time average between and on the lower (b) panel. The corresponding calibration bars of comoving are shown on the right. Magnetic field energy density is normalized by its maximum value. The magnetic field lines illustrate only the direction of the field’s poloidal component.

4. Dynamical Model Based on Simulations

We now discuss extensions of the numerical simulations, which we need to perform radiative transfer computations. We extend the simulations to large radii and define the electron temperature.

4.1. Extension to Large Radius

The flow is evolved in a quasi-steady state for from until , which corresponds to orbits at . The flow is not sufficiently settled at larger radii. However, outside , some Faraday rotation and emission might occur. So, we extend the dynamical model to larger radii (i.e. ) in a continuous way and check (see Appendix A) how variations of our large radius prescriptions change the results of radiative transfer. The outer radial boundary of radiative transfer is situated at . The profiles of number density , internal energy density , magnetic field , and velocity are extended as power-laws until radius . The power-law index for number density is obtained by matching the known value at about (Baganoff et al., 2003) and the average value at in the equatorial plane for each model. The value of may be different for different models. The radial flow velocity is then obtained from the continuity relation in the equatorial plane . The power-law of internal energy density is obtained in a similar way by matching the values  K and at distance (Baganoff et al., 2003; Shcherbakov & Baganoff, 2010). The meridional physical velocity is extended as and toroidal velocity as to approximately match the power law between and , where the relationship is used to connect the 4-velocity components with physical velocity components. All components of comoving magnetic field are extended as , which appears valid across a diverse set of GRMHD models (McKinney et al., 2012). This power-law slope corresponds roughly to equipartition of magnetic field energy density, since constant fraction magnetic field is for . Exploration of various extensions of the magnetic field will be the topic of future studies.

After defining the extension power-laws for quantities in the equatorial plane, we extend the quantities radially at arbitrary and in a continuous way. For example, for density at arbitrary and and at we have


where is taken from the simulations. We similarly extend other quantities. As shown in Appendix A, small variations in power-law indices of number density and temperature have little influence on radiation intensities and linear/circular polarization fluxes, but variations of magnetic field slope can make a substantial difference.

4.2. Electron Temperature

Neither the proton nor the electron temperature is given directly by the simulations. However, it is crucial to know the electron temperature to determine the emission. Our solution is to split the total internal energy density , given by the simulations and their power-law extension, between the proton energy and the electron energy. The energy balance states


where and are the respective heat capacities, is the rest-mass density, and is Boltzmann’s constant. The difference of temperatures is influenced by three effects: equilibration by Coulomb collisions at large radii, the difference in heating rates and of protons and electrons operating at intermediate radii, and the difference in heat capacities operating close to the BH. Radiative cooling is ignored since, according to Sharma et al. (2007), the radiative efficiency of the flow is negligible for realistic . The relevant effects can be summarized by the equation:




is the non-relativistic temperature equilibration rate by collisions (Shkarofsky et al., 1966), all quantities being measured in CGS units. We consider protons to always have non-relativistic heat capacity and collisions to always obey the non-relativistic formula. The magnitudes of errors introduced by these simplifications are negligible. The exact expressions for total electron heat capacity and differential heat capacity are approximated as


correspondingly, with the error , where


is the dimensionless electron temperature. It was recently shown (Sharma et al., 2007) that the ratio of heating rates in the non-relativistic regime in a disk can be approximated as


with coefficient . This formula is adopted in the relativistic regime as well, since no better prescription is available. Sharma et al. (2007) found the value in simulations, whereas we find for the best fit models (see Table 2 and § 6).

The proton and electron temperatures are determined at each point in the following way. We first take a single snapshot of a simulation with spin and extend the flow quantities to (see § 4.1). Then we compute azimuthal averages of radial velocity , number density , and at the equatorial plane, extend them as power laws to , and solve the equations (11,12) from down to the inner grid cell point. Temperatures are set to  K at (Baganoff et al., 2003; Shcherbakov & Baganoff, 2010). On the next step we compare the values of to the calculated and and determine the functional dependence and . At each point of the simulation (including off the equator), we draw temperatures from this correspondence. That is, GRMHD simulation directly provides and at the equatorial plane, so the function gives at each point in space. Typical profiles of proton and electron temperatures are shown on Figure 6. Temperatures stay equal until due to collisions, despite different heating prescriptions. Within the timescale of collisional equilibration becomes relatively long and electrons become relativistic, thus deviates down from . The electron and proton temperature profiles in the region are used to conduct the radiative transfer.

Figure 6.— Temperatures of protons (upper red line) and electrons (lower blue line) for the dynamical model with spin giving the best fit to polarization observations (see Table 2 and § 6).

For a given accretion rate, we find that there exists a unique dependence of the ratio of temperatures (measured at at the equator) upon the heating coefficient , so that we can use and interchangeably.

5. General Relativistic Polarized Radiative Transfer

5.1. Description of Radiative Transfer

Now we convert the dynamical model of the accretion flow into a set of observable quantities using polarized radiative transfer (Broderick et al., 2009; Shcherbakov & Huang, 2011). We closely follow Shcherbakov & Huang (2011) for the transfer technique. Similar to Huang et al. (2009a), we define the polarized basis in the picture plane, where one vector points North, another vector points East, and the wavevector points towards the observer. We parallel transport this basis in the direction of the BH and do the radiative transfer along the ray in the opposite direction towards the observer. At each point along the ray we go into the locally-flat comoving frame, calculate the angles between the magnetic field and basis vectors, and compute the Faraday conversion, Faraday rotation, emissivities, and absorptivities.

Radiative transfer involves shooting a uniform grid of geodesics from the picture plane down to the black hole. The total polarized fluxes are computed by integration of intensities along each ray backwards to the picture plane. We found that is good enough to compute the spectrum (Dexter et al. 2009 used ). For radiative transfer we employ all 3D data in each numerical simulation snapshot and, following Moscibrodzka et al. (2009), perform multilinear interpolation in three dimensions for the quantities in between the grid points. We make no approximations in the use of spatial 3D data. We self-consistently take into account the evolution of the numerical simulation as the light geodesics travel around the BH. Since it is too time-consuming to look up simulation data over a long period of time, we only evolve the simulation between and to get a spectrum at time . The offset appears, since the picture plane is located away from the BH center. The extension to the large radius outside , however, is not evolved with time. It is taken to be that of a single snapshot at time . The snapshot at times and are taken to represent the numerical simulations at earlier and later times, respectively. We find that is large enough to achieve accurate simulated spectra. The total fluxes are found at regular time intervals within period of quasi-steady accretion from till , e.g. for . We compute spectra over the quasi-steady accretion phase and average them to find the mean simulated spectra. To compute the polarized fluxes we take the integration domain in the picture plane to be a square with a side


in the units of , where frequency is in GHz. This square is centered at the BH. The size based on Equation (18) is larger than the photon orbit visible diameter and follows the intrinsic size dependence on frequency (Shen et al., 2005; Doeleman et al., 2008) at low frequencies. An important radiative transfer parameter is the distance from the BH, where intensity integration starts. The dependence of synchrotron emissivity on temperature and magnetic field strength is so strong that it overwhelms the sole effect of gravitational redshift close to the BH. We obtain accurate results in the sub-mm for computation out from , where is the horizon radius. To quantify the needed accuracy of computations, we define a quantity in Appendix A. We conduct multiple tests of radiative transfer convergence for best fit models at each spin. In Appendix A, we justify the chosen values of radiative transfer parameters , , , , etc.

Our calculation of plasma response is different from Shcherbakov & Huang (2011). They offered a way to find exact emissivities, absorptivities, Faraday rotation, and conversion coefficients for thermal and other isotropic particle distributions. Here, for simplicity, we employ fitting formulas for Faraday rotation and Faraday conversion and synchrotron approximation for emissivities for a thermal plasma. We define


where is - angle, is electron gamma factor, and is the cyclotron frequency. Then following Legg & Westfold (1968); Melrose (1971), we write down emissivities in the , , and modes as


Here is the Bessel function of the 2nd kind of order . We employed IEEE/IAU definitions of Stokes , , and (Hamaker & Bregman, 1996), and we define counter-clockwise rotation of the electric field as seen by the observer as corresponding to positive – as also chosen in Shcherbakov & Huang (2011). So, the sign of the emissivity (Eq. 5.1) is opposite to the sign in Rybicki & Lightman (1967). A variation of emissivity formulas (19,5.1) exists: Sazonov (1969); Pacholczyk (1970) define , integrating over particle energy instead of . This approximation appears to give significantly larger errors at low particle energies.

Next, one needs to identify the accurate thermal particle distribution . Various correspond to various synchrotron approximations. The ultrarelativistic thermal approximation (Pacholczyk, 1970; Huang et al., 2009a) has the simplest distribution . However, the exact thermal distribution of particles


allows for more precise computation of radiation. Synchrotron emissivities based on the equations (19,5.1) with the exact thermal distribution (21) agree with the exact cyclo-synchrotron emissivities , , and (Leung, Gammie & Noble, 2011; Shcherbakov & Huang, 2011) to within for typical dynamical models and frequencies  GHz. Emissivities integrated over the ultrarelativistic thermal distribution typically have error.

Thermal absorptivities are found from emissivities (Eq. 5.1) via Kirchhoff’s law


where is the source function for low photon energies (). Faraday rotation and Faraday conversion coefficients are taken from Shcherbakov (2008):






are the fitting formulas for deviations of and from analytic results for finite ratios of . The deviation of from is significant for the set of observed frequencies , temperatures , and magnetic fields found in the typical models of Sgr A*. These formulas constitute a good fit to the exact result for the typical parameters of the dynamical model (Shcherbakov, 2008).

Polarized radiative transfer can take much longer to perform compared to non-polarized radiative transfer when using an explicit integration scheme to evolve the Stokes occupation numbers , , and . Large Faraday rotation measure and Faraday conversion measure lead to oscillations between occupation numbers. One of the solutions is to use an implicit integration scheme, while another solution is to perform a substitution of variables. In the simple case of Faraday rotation leading to interchange of and , our choice of variables is the amplitude of oscillations and the phase. Thus, the cylindrical polarized coordinates arise as follows:


Then, the amplitude slowly changes along the ray and the angle changes linearly, and this translates into a speed improvement. In the presence of substantial Faraday conversion, the polarization vector precesses along some axis on a Poincaré sphere, adding an interchange of circularly and linearly polarized light. So, polar polarized coordinates are more suitable in this case:


where is the total polarized intensity, angle changes are mainly due to Faraday rotation, and angle changes are mainly due to Faraday conversion. The application of this technique speeds up the code enormously at low frequencies of  GHz.

Model Inclination angle , deg Spin position angle PA, deg Heating constant Ratio at Electron at , K Accretion rate ,
spin 42.0 171.0 0.42107 15.98
spin 74.5 115.3 0.37012 20.14
spin 64.5 84.7 0.37239 20.16
spin 53.5 123.4 0.39849 18.16
spin 57.2 120.3 0.41343 17.00
spin short period 1 70.0 79.3 0.38934 18.50
spin short period 2 72.8 113.1 0.40507 17.31
spin short period 3 73.4 57.4 0.37302 19.87
spin short period 4 74.4 115.4 0.36147 20.95
spin short period 5 71.9 95.7 0.37420 19.79
spin short period 6 76.4 116.7 0.38853 18.59
spin fast light 41.4 187.5 0.41929 16.09
spin fast light 72.7 105.9 0.39804 17.83
spin fast light 59.4 131.8 0.35708 21.62
spin fast light 53.3 123.3 0.40215 17.86
spin fast light 57.7 119.6 0.41720 16.73

Note. – Mean values are shown for ratio , electron temperature , and the accretion rate . These are the simple means over all simulation snapshots, which were employed for radiative transfer in a particular model. The values of ranges from across all models.

Table 2Properties of the best fit models with different spins.
Figure 7.— Fits to the observed fluxes, LP and CP fractions by best models for each spin. The inclination angle , accretion rate , ratio of temperatures were adjusted for each spin to minimize . Fits to total flux are in the upper left panel, LP fraction in the lower left, and CP fraction in the lower right. Shown are the best fit models with spin (short-dashed brown), spin (solid dark red), spin (long-dashed green), spin (solid light cyan), and spin (dot-dashed orange). The upper right panel shows the dependence of EVPA on frequency for the best models. Note, that EVPAs are not included into our fitting procedure. The thick blue curve represents observations. Simulated EVPA curves are arbitrarily shifted to approximate EVPA at  GHz. The addition of an external (to the emitting region) Faraday rotation screen helps to fit .

5.2. Search for the Best Fits

We define quantities to discriminate between models. We define for fitting total fluxes as


for the set of frequencies , and  GHz, where are the errors of the means. We incorporate LP fractions at , and  GHz and CP fractions at and  GHz to obtain


Then we define (as degrees of freedom) to be for flux fitting and for fitting all polarized data. The quantity would be drawn from statistics if -s were the true observational errors and if the observed fluxes were drawn from a Gaussian distribution. However, for the purpose of the present work, we only employ as a measure of fitting the data. That is, lower indicates better agreement with the data. We do not attempt to ascribe any statistical meaning to the quantity .

We explore models with parameters: spin , inclination angle , accretion rate , and the ratio of proton to electron temperature ( is reported for radius ). For the radiative transfer calculations, the density from the simulations is scaled to give the desired accretion rate.

6. Results

In previous sections, we described our compiled observations, GRMHD numerical simulations of the flow structure, our method for obtaining the electron temperature, and our method for polarized radiative transfer. In this section, we discuss our results for accretion flow and BH parameters, as guided by a minimization of for our model applied to the observations.

Figure 7 shows best fits to observations by models with five different spins. Inclination angle , accretion rate , and heating coefficient were adjusted to reach the lowest . Fits to fluxes (upper left) are not substantially different, although models with higher spins fit better at high frequencies. Larger deviations can be seen on (lower left) and (lower right) plots. Models with high spins require lower accretion rate (i.e. density) to fit the flux spectrum. As a consequence, they are not subject to Faraday depolarization, which leads to a decrease of LP at low , and the models end up having larger linear polarization fractions at  GHz. Not all models reproduce the observed decrease of mean LP fraction between  GHz and  GHz groups. The discrepancies in fitting the CP fraction are also large: all the lowest models give at  GHz. The best bet model with spin reproduces and fractions well, but fails in fitting the total flux. Most solutions predict the wrong sign of the difference, which could be fixed with stronger magnetic field (e.g. as seen in models by McKinney et al. 2012) to yield stronger Faraday rotation. In sum, crude agreement of simulated polarized spectra to the observed ones was achieved, but the improved dynamical models may be needed for better fits.

We now isolate the physical effects responsible for the observed polarized quantities for our best bet model with spin that has the lowest (see subsection 6.1).

Figure 8.— Contributions of different effects to the CP fraction as a function of frequency for our best bet model with BH spin . Shown are observations (blue error bars), the best bet model (solid red line), the same dynamical model computed with zero V emissivity () in radiative transfer so that CP is produced by Faraday conversion (dot-dashed orange), the same model with zero Faraday conversion () (short-dashed brown), and the same model with zero Faraday rotation () (long-dashed green). Emissivity in circular V mode contributes little to the observed CP, which is mainly due to Faraday conversion.
Figure 9.— Contributions of different effects to the LP fraction (on the left) and EVPA (on the right) as functions of frequency for the best bet model with spin . Shown are observations (blue error bars and thick blue line), the best bet model (solid red line), and the same dynamical model computed with zero Faraday rotation () in radiative transfer (long-dashed green). Beam depolarization is weak: if Faraday rotation is absent, then LP stays high at low frequencies. Even when the Faraday rotation is set to zero, the EVPA depends on frequency due to varying intrinsic emission EVPA. Faraday rotation in the best bet model is too weak to reproduce EVPA observations, so stronger magnetic fields or more magnetic flux near the black hole than in the simulations may be required.

There are several radiative transfer effects that contribute similarly to the polarized fluxes. Let us consider the production of circular polarization in the flow. Figure 8 shows the consequence of switching off each physical effect for our best bet model with spin . The solid red curve is the result with all physics on. The dot-dashed orange line below is for zero circular emissivity having . The brown dashed line corresponds to zero Faraday conversion (). Switching off emissivity leads to a minor correction, whereas setting Faraday conversion to zero results in CP of the opposite sign with several times smaller absolute value. Most of the CP in this model is produced by Faraday conversion. It would be incorrect, however, to think that the simple linear to circular conversion explains the observed CP. The dashed green line in Figure 8 shows the CP fraction, when Faraday rotation is switched off (). The effect of Faraday rotation is insignificant at  GHz, but the rotation of the plane of linear polarization simultaneous with conversion between linear and circular polarizations produces a unique effect at lower . This is the so-called “rotation-induced conversion” (Homan et al., 2009). Sign oscillations of with frequency do not happen when the Faraday rotation is on, but they do happen when . For the best fit model it is the rotation-induced Faraday rotation, which is responsible for the most of circularly polarized light.

In Figure 9 we illustrate the influence of Faraday rotation on LP fraction (left panel) and EVPA (right panel). The solid curves are produced with all physics on for our best bet model with spin . The green dashed lines are computed when switching off Faraday rotation (). The Faraday rotation is small at high frequencies and curves look similar at  GHz. As the rotation of polarization plane is much stronger at low , a significant phase shift accumulates between different rays at the low end of the spectrum and cancellations of LP become strong at  GHz. This illustrates the effect of Faraday depolarization (Bower et al., 1999a). In the absence of Faraday rotation, the dependence of EVPA on frequency is not constant: the variations of intrinsic emitted EVPA are significant. Thus, the change of EVPA with should not always be ascribed to the effect of Faraday rotation. The positive observed slope of EVPA with at high , acquired due to negative Faraday rotation measure (), is comparable to the slope of intrinsic emitted EVPA.

Figure 10.— Correlated flux as a function of baseline at  GHz normalized to the averaged observed flux at  Jy for the best fit models with spin (darker red lines) and (lighter orange lines). The upper solid lines show the smallest size (largest correlated flux) over all position angles of BH spin axis, and the lower dashed lines show the large