Spectral and Imaging properties of Sgr A* from High-Resolution 3D GRMHD Simulations with Radiative Cooling
Doosoo Yoon, Koushik Chatterjee, Sera Markoff, David van Eijnatten, Ziri Younsi, Matthew Liska, Alexander Tchekhovskoy
MMNRAS , 1–16 (2020) Preprint 1 October 2020 Compiled using MNRAS L A TEX style file v3.0
Spectral and Imaging properties of Sgr A ∗ from High-Resolution 3DGRMHD Simulations with Radiative Cooling D. Yoon (cid:63) , K. Chatterjee , S.B. Marko ff , , D. van Eijnatten , Z. Younsi , ,M. Liska , & A. Tchekhovskoy Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands Gravitation and Astroparticle Physics Amsterdam (GRAPPA) Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, United Kingdom Institut für Theoretische Physik, Goethe-Universität Frankfurt, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany Institute for Theory and Computation, Harvard University, 60 Garden Street, Cambridge, MA 02138, USA; John Harvard Distinguished Science and ITCFellow Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The candidate supermassive black hole in the Galactic Centre, Sagittarius A* (Sgr A ∗ ), isknown to be fed by a radiatively ine ffi cient accretion flow (RIAF), inferred by its low ac-cretion rate. Consequently, radiative cooling has in general been overlooked in the study ofSgr A ∗ . However, the radiative properties of the plasma in RIAFs are poorly understood. Inthis work, using full 3D general-relativistic magneto-hydrodynamical simulations, we studythe impact of radiative cooling on the dynamical evolution of the accreting plasma, presentingspectral energy distributions and synthetic sub-millimeter images generated from the accretionflow around Sgr A ∗ . These simulations solve the approximated equations for radiative cool-ing processes self-consistently, including synchrotron, bremsstrahlung, and inverse Comptonprocesses. We find that radiative cooling plays an increasingly important role in the dynam-ics of the accretion flow as the accretion rate increases: the mid-plane density grows and theinfalling gas is less turbulent as cooling becomes stronger. The changes in the dynamical evo-lution become important when the accretion rate is larger than 10 − M (cid:12) yr − ( (cid:38) − ˙ M Edd ,where ˙ M Edd is the Eddington accretion rate). The resulting spectra in the cooled models alsodi ff er from those in the non-cooled models: the overall flux, including the peak values at thesub-mm and the far-UV, is slightly lower as a consequence of a decrease in the electron tem-perature. Our results suggest that radiative cooling should be carefully taken into account inmodelling Sgr A ∗ and other low-luminosity active galactic nuclei that have a mass accretionrate of ˙ M > − ˙ M Edd . Key words: galaxies: black hole physics – accretion, accretion disks, jets – galaxies: individ-ual (SgrA*) – magnetohydrodynamics (MHD) – methods: numerical
It is widely believed that most galaxies harbour supermassive blackholes (SMBHs) in their galactic centres, with masses ranging frommillions to billions of solar masses. Over the past few decades, theblack hole (BH) candidate in the centre of the Milky Way, Sagittar-ius A* (hereafter Sgr A ∗ ), has proven an exceptional laboratory forstudies of accretion and outflow physics of BHs due to its proximityto Earth. A significant e ff ort has been invested in determining theBH mass and distance for Sgr A ∗ (e.g., Reid 1993; Reid et al. 2019;Schödel et al. 2002; Bower et al. 2004; Ghez et al. 2008; Gillessen (cid:63) E-mail: [email protected] et al. 2009, 2017; Boehle et al. 2016; Gravity Collaboration et al.2018). We adopt the current best-fit BH mass of 4 . ± . × M (cid:12) ,which was measured by the orbital motion of stars and gas clouds(Gillessen et al. 2017; Gravity Collaboration et al. 2018), and thedistance of 8 . ± .
15 kpc, which was obtained from trigonomet-ric parallaxes and proper motions of massive stars around Sgr A ∗ (Reid et al. 2019). Given the mass and the distance, the angular sizeof the Schwarzschild radius, r S = GM BH / c , is ≈ µ as, whichsubtends a larger area in the sky than any other known BH, includ-ing all stellar-mass BHs.Recently, mounting attention has been paid to the studyof Sgr A ∗ with the advent of the pioneering instruments GRAVITY (Gravity Collaboration et al. 2017, 2018) and the
Event c (cid:13) a r X i v : . [ a s t r o - ph . H E ] S e p Yoon et al.
Horizon Telescope ( EHT ; Doeleman et al. 2009; Event Hori-zon Telescope Collaboration et al. 2019), capable of probing 10 –30 µ as scales. These instruments allow us to profoundly improveour understanding of the physical processes associated with the ac-cretion and relativistic jet formation in the immediate vicinity ofSMBHs and thus demand equal measures of theoretical supportand precise modelling of the spectrum generated by the radiationfrom the accretion flow around Sgr A ∗ .The mass accretion rate around Sgr A ∗ (in units of M (cid:12) yr − ) isestimated to be in the range of ∼ − < ˙ M < − , as constrainedby the measured Faraday rotation measure at mm / sub-mm wave-lengths (Aitken et al. 2000; Bower et al. 2003; Marrone et al. 2007).Such a low accretion rate favours hot accretion flow models for theaccretion disk instead of the radiatively e ffi cient, thin disk models(Shakura & Sunyaev 1973). Many theoretical scenarios have beeninvoked and excluded to account for the nature of accretion andoutflows in the hot accretion flow: the standard advection domi-nated accretion flow (ADAF; Narayan & Yi 1994; Narayan et al.1995, 1998) and Bondi-Hoyle models are ruled out, since thesemodels are expected to yield an accretion rate of ∼ − M (cid:12) yr − ,which is two orders of magnitude higher than the measured upperlimit (Bower et al. 2003). Yuan et al. (2003) reexamined the radia-tively ine ffi cient accretion flow (RIAF) model for the spectrum ofSgr A ∗ and argued that the presence of outflows within the Bondiradius plays a vital role in reducing the mass accretion rate. Al-ternatively, the convection dominated accretion flow (CDAF) andjet-dominated models are capable of reproducing the spectrum thatis consistent with the observed accretion rate (Quataert & Gruzinov2000; Falcke & Marko ff ff et al. 2007). The existenceof such outflows is supported by the weak hydrogen-like Fe K α line around Sgr A ∗ via the X-ray Visionary Program (Wang et al.2013). In this study, the flat density profile in the spectrum con-firmed that (cid:38)
99% of the matter initially captured by the SMBH islost before it reaches the innermost region around Sgr A ∗ , which isconsistent with the CDAF model or adiabatic inflow-outflow solu-tion (ADIOS; see Blandford & Begelman 1999; Yuan & Narayan2014 for the detailed model descriptions) model.Although semi-analytic models provide an important frame-work for understanding the nature of the accretion flow aroundSgr A ∗ , numerical simulations are required to capture the time-dependent, turbulent evolution of the accretion flow. In partic-ular, a self-consistent magneto-hydrodynamical (MHD) descrip-tion enables us to demonstrate accretion processes induced by themagneto-rotational instability (MRI; Balbus & Hawley 1991) with-out imposing an arbitrary anomalous viscosity. In earlier numer-ical studies, three-dimensional pseudo-Newtonian MHD simula-tions were carried out to model the synchrotron radiation from ac-cretion flows(e.g., Goldston et al. 2005; Ohsuga et al. 2005; Chanet al. 2009; Huang et al. 2009). However, the non-relativistic treat-ment in the simulations has disadvantages for modelling the syn-chrotron radiation, since it is mainly emitted in the immediatevicinity of the central BH, where relativistic e ff ects cannot be ig-nored: shocks develop di ff erently for the relativistic plasma whensubject to intense magnetic and gravitational fields (Del Zannaet al. 2003), and the curvature of space-time becomes signifi-cant. Several other works made use of general-relativistic magneto-hydrodynamical (GRMHD) simulations for studying the dynamicsand spectral properties of Sgr A ∗ in two dimensions (2.5D; e.g,Noble et al. 2007; Mo´scibrodzka et al. 2009; Hilburn et al. 2010;Dibi et al. 2012; Drappeau et al. 2013; Mo´scibrodzka & Falcke2013), or in three dimensions (3D) (e.g., Dexter et al. 2009, 2010;Dolence et al. 2012; Shcherbakov et al. 2012; Dexter & Fragile 2013; Mo´scibrodzka et al. 2014; Davelaar et al. 2018; Chael et al.2018). In general, 2.5D simulations are a reasonable and computa-tionally cheaper option to conduct a parameter study for reproduc-ing the spectrum of Sgr A ∗ , but it is known that simulations withaxisymmetric coordinates cannot sustain MRI-driven turbulence,which decays over the local orbital time as a consequence of Cowl-ing’s anti-dynamo theorem (Cowling 1933, see also Hide & Palmer1982 for the generalised description). Therefore, full 3D GRMHDsimulations are necessary to perform detailed studies of the natureof the accretion flows around Sgr A ∗ and the emitted spectrum.For instance, it was confirmed that thick accretion disks are ableto generate and advect large-scale poloidal magnetic flux throughdynamo action when resolved properly (Liska et al. 2018a).The bolometric luminosity of Sgr A ∗ is extremely low, L bol ∼ erg s − ≈ − L Edd , where L Edd is the Eddington luminosity.Given such a low luminosity, it has been thought that the radiativecooling losses of Sgr A ∗ are negligible, since the losses are likelynot strong enough to have a significant impact on the dynamicsof the accretion flow. Based on this argument, all previous workswith full 3D GRMHD simulations have ignored the radiative cool-ing losses for simplicity. Although this assumption may be reason-able, Dibi et al. (2012) argued based on their 2.5D simulations thatcooling losses play an increasingly important role for higher accre-tion rates and possibly alter the dynamics and resulting spectra ofSgr A ∗ , even within the allowed range of accretion rates based onpolarisation and X-ray studies. One potential impact is that the ra-diative cooling reduces the gas pressure and the disk vertical scaleheight, resulting in a decrease in turbulence and a more orderedmagnetic field (Fragile & Meier 2009). Moreover, many questionsremain unanswered: how do radiative cooling losses a ff ect the tur-bulence features of the disk, and thus the angular momentum trans-fer of the accreting plasma? How does radiative cooling togetherwith GR e ff ects result in the observed spectra from Sgr A ∗ ? Is ra-diative cooling indeed negligible for the mass accretion rate rangeof Sgr A ∗ ? Even though the e ff ects of radiative cooling can be mi-nor for the case of Sgr A ∗ , the quantitative evaluation of coolinge ff ects is highly demanded because it must play a greater role forSMBHs with higher mass accretion rates, such as M87.In this paper, we perform the first full 3D GRMHD simu-lations which include cooling losses via bremsstrahlung, thermalsynchrotron emission, and inverse Compton scattering. Due to thesignificant computational expense of the full 3D simulations, wecannot explore the full range of various parameters (e.g., spin,magnetic configuration, electron distribution function, misaligneddisk). Instead, we use parameters compatible with earlier studies,assuming a rapidly rotating BH, weak poloidal initial magneticfield, and a fixed temperature ratio between ions and electrons of T i / T e = ff erentelectron temperature prescriptions for comparison). We then exam-ine the e ff ect of radiative cooling on the dynamics of the accretionflow and the resulting spectra and images, for di ff erent accretionrates within the allowed range.This paper is structured as follows. In § 2, we give a technicaldescription of the numerical methods used, including the simula-tion setup, and the treatment of radiative cooling losses. In § 3, wedescribe the results of how cooling losses play a role in chang-ing the dynamical evolution of accreting matter. In § 4 we discussthe best-bet model for Sgr A ∗ , the e ff ects of cooling on the re-sulting spectra and sub-mm images, and the variability of multi-wavelength spectra. We also compare our 3D work to previous2.5D work. We summarise our results in § 5. MNRAS000
99% of the matter initially captured by the SMBH islost before it reaches the innermost region around Sgr A ∗ , which isconsistent with the CDAF model or adiabatic inflow-outflow solu-tion (ADIOS; see Blandford & Begelman 1999; Yuan & Narayan2014 for the detailed model descriptions) model.Although semi-analytic models provide an important frame-work for understanding the nature of the accretion flow aroundSgr A ∗ , numerical simulations are required to capture the time-dependent, turbulent evolution of the accretion flow. In partic-ular, a self-consistent magneto-hydrodynamical (MHD) descrip-tion enables us to demonstrate accretion processes induced by themagneto-rotational instability (MRI; Balbus & Hawley 1991) with-out imposing an arbitrary anomalous viscosity. In earlier numer-ical studies, three-dimensional pseudo-Newtonian MHD simula-tions were carried out to model the synchrotron radiation from ac-cretion flows(e.g., Goldston et al. 2005; Ohsuga et al. 2005; Chanet al. 2009; Huang et al. 2009). However, the non-relativistic treat-ment in the simulations has disadvantages for modelling the syn-chrotron radiation, since it is mainly emitted in the immediatevicinity of the central BH, where relativistic e ff ects cannot be ig-nored: shocks develop di ff erently for the relativistic plasma whensubject to intense magnetic and gravitational fields (Del Zannaet al. 2003), and the curvature of space-time becomes signifi-cant. Several other works made use of general-relativistic magneto-hydrodynamical (GRMHD) simulations for studying the dynamicsand spectral properties of Sgr A ∗ in two dimensions (2.5D; e.g,Noble et al. 2007; Mo´scibrodzka et al. 2009; Hilburn et al. 2010;Dibi et al. 2012; Drappeau et al. 2013; Mo´scibrodzka & Falcke2013), or in three dimensions (3D) (e.g., Dexter et al. 2009, 2010;Dolence et al. 2012; Shcherbakov et al. 2012; Dexter & Fragile 2013; Mo´scibrodzka et al. 2014; Davelaar et al. 2018; Chael et al.2018). In general, 2.5D simulations are a reasonable and computa-tionally cheaper option to conduct a parameter study for reproduc-ing the spectrum of Sgr A ∗ , but it is known that simulations withaxisymmetric coordinates cannot sustain MRI-driven turbulence,which decays over the local orbital time as a consequence of Cowl-ing’s anti-dynamo theorem (Cowling 1933, see also Hide & Palmer1982 for the generalised description). Therefore, full 3D GRMHDsimulations are necessary to perform detailed studies of the natureof the accretion flows around Sgr A ∗ and the emitted spectrum.For instance, it was confirmed that thick accretion disks are ableto generate and advect large-scale poloidal magnetic flux throughdynamo action when resolved properly (Liska et al. 2018a).The bolometric luminosity of Sgr A ∗ is extremely low, L bol ∼ erg s − ≈ − L Edd , where L Edd is the Eddington luminosity.Given such a low luminosity, it has been thought that the radiativecooling losses of Sgr A ∗ are negligible, since the losses are likelynot strong enough to have a significant impact on the dynamicsof the accretion flow. Based on this argument, all previous workswith full 3D GRMHD simulations have ignored the radiative cool-ing losses for simplicity. Although this assumption may be reason-able, Dibi et al. (2012) argued based on their 2.5D simulations thatcooling losses play an increasingly important role for higher accre-tion rates and possibly alter the dynamics and resulting spectra ofSgr A ∗ , even within the allowed range of accretion rates based onpolarisation and X-ray studies. One potential impact is that the ra-diative cooling reduces the gas pressure and the disk vertical scaleheight, resulting in a decrease in turbulence and a more orderedmagnetic field (Fragile & Meier 2009). Moreover, many questionsremain unanswered: how do radiative cooling losses a ff ect the tur-bulence features of the disk, and thus the angular momentum trans-fer of the accreting plasma? How does radiative cooling togetherwith GR e ff ects result in the observed spectra from Sgr A ∗ ? Is ra-diative cooling indeed negligible for the mass accretion rate rangeof Sgr A ∗ ? Even though the e ff ects of radiative cooling can be mi-nor for the case of Sgr A ∗ , the quantitative evaluation of coolinge ff ects is highly demanded because it must play a greater role forSMBHs with higher mass accretion rates, such as M87.In this paper, we perform the first full 3D GRMHD simu-lations which include cooling losses via bremsstrahlung, thermalsynchrotron emission, and inverse Compton scattering. Due to thesignificant computational expense of the full 3D simulations, wecannot explore the full range of various parameters (e.g., spin,magnetic configuration, electron distribution function, misaligneddisk). Instead, we use parameters compatible with earlier studies,assuming a rapidly rotating BH, weak poloidal initial magneticfield, and a fixed temperature ratio between ions and electrons of T i / T e = ff erentelectron temperature prescriptions for comparison). We then exam-ine the e ff ect of radiative cooling on the dynamics of the accretionflow and the resulting spectra and images, for di ff erent accretionrates within the allowed range.This paper is structured as follows. In § 2, we give a technicaldescription of the numerical methods used, including the simula-tion setup, and the treatment of radiative cooling losses. In § 3, wedescribe the results of how cooling losses play a role in chang-ing the dynamical evolution of accreting matter. In § 4 we discussthe best-bet model for Sgr A ∗ , the e ff ects of cooling on the re-sulting spectra and sub-mm images, and the variability of multi-wavelength spectra. We also compare our 3D work to previous2.5D work. We summarise our results in § 5. MNRAS000 , 1–16 (2020) pectra of Sgr A ∗ from 3D radiative GRMHD All simulations are performed with the
H-AMR code (Liska et al.2019a; Porth et al. 2019), which branched o ff HARM2D (Gammieet al. 2003; Noble et al. 2006) in its early days. It is accelerated byGraphical Processing Units (GPUs) and improved with a staggeredgrid for constrained transport of magnetic fields (Gardiner & Stone2005) to preserve ∇ · B =
0, more robust inversion (Newman &Hamlin 2014) adaptive mesh refinement (AMR, not utilised in thiswork), static mesh refinement (SMR), and a locally adaptive timestep (LAT; see Chatterjee et al. 2019, Appendix A). It adopts apiece-wise parabolic method (PPM; Colella & Woodward 1984)for reconstruction of cell-centred quantities at cell faces, which isthird-order accurate, for the spatial reconstruction at cell faces fromcell centres, and a second-order time-stepping.The broadband spectrum is calculated from the GRMHDoutput, using the general-relativistic Monte Carlo scheme
GRMONTY (Dolence et al. 2009), which includes synchrotron emis-sion and absorption, and inverse Compton scattering for a rela-tivistic thermal Maxwell-Jüttner distribution of electrons. Techni-cally,
GRMONTY cannot produce synthetic images but only spectra.Thus, we ray-trace the GRMHD-produced spectra by integratingthe general-relativistic radiative transfer (GRRT) equations usingthe
BHOSS code (Younsi et al. 2012, 2016, 2020b) to generate syn-thetic images at 230 GHz that can help us infer the expected imagesof Sgr A ∗ from the upcoming EHT project. In
BHOSS , the cal-culation of radiative processes includes synchrotron emission andabsorption only, which is su ffi cient for imaging at the EHT fre-quency of 230 GHz. Since the sub-mm regime of the spectra aredominated by synchrotron emission, which both codes calculate,we verify consistency in our spectral calculations from both codesby comparing the resulting spectra in the radio to NIR bands (seeAppendix A). For convenience, we adopt Heaviside-Lorentz units, which absorba factor of √ π for the magnetic field 4-vector, b µ , so that themagnetic pressure is P B ≡ b /
2. Furthermore, the typical natu-ral units are used, GM = c =
1, which sets the length unit tobe the gravitational radius, r g ≡ G M / c , and the time unit tobe the light crossing time, t g ≡ G M / c , where G , M , c are thegravitational constant, BH mass, and the speed of light, respec-tively. We use a spherical-polar axisymmetric computational grid( r , θ, φ ) extending from 0.85 r H to 250 r g , where the event hori-zon radius r H ≡ r g (cid:16) + (cid:112) − a (cid:63) (cid:17) . Here we set the dimensionlessBH spin parameter to a (cid:63) ≡ c J / GM in a Kerr-Schild foliation,where J = r H c is the angular momentum at the event horizon.The grid is uniformly spaced with respect to a set of internal coordi-nates (cid:16) x , x , x (cid:17) , which can be converted to ( r , θ, φ ), respectively . This conversion leads to a logarithmic spacing in r such that thecells have a higher resolution for smaller values of r . The spatialresolution near the event horizon is ∆ r ≈ . r g for the model withthe highest resolution. To prevent the aspect ratio of the cells frombecoming too large near the polar singularity, we reduce the reso-lution in φ − direction gradually towards both poles. We use outflowboundary conditions for both inner and outer radial boundaries, and The coordinate transformation is made using the following relations: t = x , r = exp ( x ), θ = π x , and φ = x . See Appendix B in Chatterjee et al.(2019) for the detailed coordinate conversion. reflecting boundary conditions in the θ − direction. Note that the in-ner boundary is causally disconnected from the flow, as it is locatedwithin the event horizon.It is common for GRMHD simulations to crash if either thedensity or the internal energy become very low, particularly inthe funnel region along the polar axes or near the outer radialboundaries. To avoid this, we apply numerical floors for the den-sity and the internal energy (see Appendix B3 of Ressler et al.2017 for more detailed discussions): a minimum rest mess density, ρ fl = max (cid:20) b / , u g / , − (cid:16) r / r g (cid:17) − (cid:21) , and a minimum internalenergy density, u g , fl = max (cid:20) b / , − (cid:16) r / r g (cid:17) − Γ (cid:21) , where b and u g are the co-moving magnetic field strength and the internal energydensity, respectively. We normalise the mass density such that themaximum density is ρ max = . We perform a set of GRMHD simulations, in which the magnetisedgas is accreting onto a supermassive and spinning BH. All simula-tions are initialised with a steady-state hydrostatic torus around arapidly spinning Kerr BH (Fishbone & Moncrief 1976). We set thespin parameter to a (cid:63) = . r in = r g , and the radius of thepressure maximum, r max = r g . We adopt an ideal gas equation ofstate, P g = ( Γ − u g , (1)where P g and u g are the gas pressure and internal energy, respec-tively. We set the adiabatic index to Γ = /
3, which assumes thedominance of a non-relativistic plasma in the accretion flow.As an initial magnetic configuration, we adopt a single loop ofweak magnetic field, which is computed from the magnetic vectorpotential, A φ ∝ (cid:40) ρ − . , if ρ > . , , otherwise . (2)The centre of the loop is at the density maximum, and the loop isfully contained within the initial torus. The initial magnetic fieldis normalised such that β mag = P g / P B ≥ We take into account the radiative cooling self-consistently in ourcalculation of the gas temperature, by including the e ff ects ofbremsstrahlung, synchrotron, and the inverse Compton losses. Weadopt the equations of Esin et al. (1996) for computing the radiativecooling losses. These formulae have been implemented and testedin previous numerical studies of Sgr A ∗ (Fragile & Meier 2009;Dibi et al. 2012; Straub et al. 2012; Drappeau et al. 2013).The total cooling rate for an optically thin gas is computedfrom the cooling function, q − thin = η br , C q − br + η s , C q − s , (3)where q − br and q − s are the bremsstrahlung and synchrotron coolingrates, respectively, and η br , C and η s , C are the Compton enhancementfactors, which are the average energy gain of the photon in an as-sumption of single scattering (Esin et al. 1996). We note that the MNRAS , 1–16 (2020)
Yoon et al.
Table 1.
Simulation setup parameters. All models are initialised with a ∗ = . r in = r g , and r max = r g .Model Name Cooling ρ a scale T i / T e (cid:104) ˙ M (cid:105) b (10 − M (cid:12) yr − ) Resolution C3D01RM on 1 × − . ± .
01 256 × × C3D1RL on 1 × − . ± .
07 128 × × C3D1RM on 1 × − . ± .
31 256 × × C3D1RH on 1 × − . ± .
17 648 × × C3D1RMFT20 on 1 × −
20 1 . ± .
20 256 × × C3D1RMRh20 on 1 × − R l = R h =
20 1 . ± .
11 256 × × C3D10RM on 1 × − . ± .
97 256 × × C3D100RM on 1 × − . ± .
18 256 × × NC3RM o ff – – – 256 × × NC2RH o ff – – – 648 × × ca conversion factor for the mass density from code units to c.g.s units. b mass accretion rate at the event horizon, which is averaged over 3000 – 8000 t g . c axisymmetric 2.5D run for the purpose of comparison. Compton enhancement of the bremsstrahlung is negligible as syn-chrotron is dominant at the temperature where the Comptonizationbecomes important.While the whole system is generally optically thin, we use thefollowing generalised cooling formula, from Narayan & Yi (1995)and Esin et al. (1996), to reproduce the equilibrium solution cor-responding to the optically thick disk (Shakura & Sunyaev 1973): q − = σ T T e / H T . τ + √ + τ − , (4)where σ T and T e are the Thomson cross-section and the electrontemperature, respectively, and the local temperature scale height H T is computed from H T = T e |∇ ( T e ) | . (5)The scale heights are locally calculated such that T e drops o ff bya factor of 1 / e , which was adopted in Fragile & Meier (2009) as asuitable and robust treatment in multi-dimensional simulations.The total optical depth of the disk is calculated by τ = τ es + τ abs , where τ es = n e σ T H T is the Thomson optical depth in thevertical direction and τ abs is the optical depth for absorption, whichis expressed as τ abs = H T q − thin σ T T e . (6)For a small optical depth, Eq. (4) reduces to Eq. (3), while, in theoptically thick limit ( τ (cid:29) q − = σ T T e / H T τ , which isthe appropriate black body limit. Therefore, the formula providesan approximate interpolation between the optically thin and thicklimits.At low temperatures ( T e (cid:46) × K) or the outer torus re-gions, the emission is dominated by bremsstrahlung (Straub et al.2012). The bremsstrahlung cooling rate is computed by the inter-actions of pairs among electrons ( e − ), positrons ( e + ) and ions ( i ).Since the cooling processes of e − i and e + i are identical, and thesame is true for e − e − and e + e + , the cooling rate can be written as, q − br = q − ei + q − ee + q −± , (7)where q − ei , q − ee and q −± are the radiative cooling throughelectron(positron)-ion ( e ± i ), electron(positron)-electron(positron)( e ± e ± ) and positron-electron ( e + e − ) interactions, respectively (Esinet al. 1996). However, for most regions of inner hot accretion flows, thesynchrotron emission dominates the losses as the electrons arerelativistic due to the high electron temperature. The synchrotroncooling occurs through both optically thick and thin emission: be-low some critical frequency ν c , the emission is completely self-absorbed, and thus the volume emissivity can be approximated bythe Rayleigh-Jeans black body emission. For frequencies above ν c ,the emission is optically thin. The synchrotron cooling rate can bewritten as, q − s = π k B T e H T c (cid:90) ν c ν d ν + (cid:90) ∞ ν c (cid:15) s ( ν ) d ν , (8)where k B and c are the Boltzmann constant and the speed of light,respectively, and the synchrotron emissivity (cid:15) s ( ν ) is calculated as(see Pacholczyk 1970), (cid:15) s ( ν ) = e c √ π ν ( n e + n + ) K (1 / Θ e ) I (cid:48) ( x M ) , (9)where K is the modified Bessel function of the second kind, x M ≡ ν ν Θ e , v ≡ e B π m e c , (10)and Θ e ≡ k B T e / m e c is the dimensionless electron temperature.The dimensionless spectrum I (cid:48) ( x M ), which is averaged over the an-gle between the velocity vector of the electron and the direction ofthe local magnetic field, is fitted by the function (Mahadevan et al.1996), I (cid:48) ( x M ) = . x / M + . x / M + . x / M exp (cid:16) − . x / M (cid:17) . (11)Fragile & Meier (2009) found that the Bessel function K in Eq. (9)causes errors for the low-temperature flows ( T e < K) due to themismatch of the normalisation factor between the Bessel functionand the spectrum I (cid:48) ( x M ). Following their suggested modification,we replace K (1 / Θ e ) by 2 Θ e , thereby assuming the same high-temperature limit.We numerically compute ν c in Eq. (8) by equating the opti-cally thin and thick volume emissivities at ν c , (cid:15) s ( ν c ) = e c √ π ν c ( n e + n + ) K (1 / Θ e ) I (cid:48) ( x M ) = π k B T e H T c ν c . (12) MNRAS000
18 256 × × NC3RM o ff – – – 256 × × NC2RH o ff – – – 648 × × ca conversion factor for the mass density from code units to c.g.s units. b mass accretion rate at the event horizon, which is averaged over 3000 – 8000 t g . c axisymmetric 2.5D run for the purpose of comparison. Compton enhancement of the bremsstrahlung is negligible as syn-chrotron is dominant at the temperature where the Comptonizationbecomes important.While the whole system is generally optically thin, we use thefollowing generalised cooling formula, from Narayan & Yi (1995)and Esin et al. (1996), to reproduce the equilibrium solution cor-responding to the optically thick disk (Shakura & Sunyaev 1973): q − = σ T T e / H T . τ + √ + τ − , (4)where σ T and T e are the Thomson cross-section and the electrontemperature, respectively, and the local temperature scale height H T is computed from H T = T e |∇ ( T e ) | . (5)The scale heights are locally calculated such that T e drops o ff bya factor of 1 / e , which was adopted in Fragile & Meier (2009) as asuitable and robust treatment in multi-dimensional simulations.The total optical depth of the disk is calculated by τ = τ es + τ abs , where τ es = n e σ T H T is the Thomson optical depth in thevertical direction and τ abs is the optical depth for absorption, whichis expressed as τ abs = H T q − thin σ T T e . (6)For a small optical depth, Eq. (4) reduces to Eq. (3), while, in theoptically thick limit ( τ (cid:29) q − = σ T T e / H T τ , which isthe appropriate black body limit. Therefore, the formula providesan approximate interpolation between the optically thin and thicklimits.At low temperatures ( T e (cid:46) × K) or the outer torus re-gions, the emission is dominated by bremsstrahlung (Straub et al.2012). The bremsstrahlung cooling rate is computed by the inter-actions of pairs among electrons ( e − ), positrons ( e + ) and ions ( i ).Since the cooling processes of e − i and e + i are identical, and thesame is true for e − e − and e + e + , the cooling rate can be written as, q − br = q − ei + q − ee + q −± , (7)where q − ei , q − ee and q −± are the radiative cooling throughelectron(positron)-ion ( e ± i ), electron(positron)-electron(positron)( e ± e ± ) and positron-electron ( e + e − ) interactions, respectively (Esinet al. 1996). However, for most regions of inner hot accretion flows, thesynchrotron emission dominates the losses as the electrons arerelativistic due to the high electron temperature. The synchrotroncooling occurs through both optically thick and thin emission: be-low some critical frequency ν c , the emission is completely self-absorbed, and thus the volume emissivity can be approximated bythe Rayleigh-Jeans black body emission. For frequencies above ν c ,the emission is optically thin. The synchrotron cooling rate can bewritten as, q − s = π k B T e H T c (cid:90) ν c ν d ν + (cid:90) ∞ ν c (cid:15) s ( ν ) d ν , (8)where k B and c are the Boltzmann constant and the speed of light,respectively, and the synchrotron emissivity (cid:15) s ( ν ) is calculated as(see Pacholczyk 1970), (cid:15) s ( ν ) = e c √ π ν ( n e + n + ) K (1 / Θ e ) I (cid:48) ( x M ) , (9)where K is the modified Bessel function of the second kind, x M ≡ ν ν Θ e , v ≡ e B π m e c , (10)and Θ e ≡ k B T e / m e c is the dimensionless electron temperature.The dimensionless spectrum I (cid:48) ( x M ), which is averaged over the an-gle between the velocity vector of the electron and the direction ofthe local magnetic field, is fitted by the function (Mahadevan et al.1996), I (cid:48) ( x M ) = . x / M + . x / M + . x / M exp (cid:16) − . x / M (cid:17) . (11)Fragile & Meier (2009) found that the Bessel function K in Eq. (9)causes errors for the low-temperature flows ( T e < K) due to themismatch of the normalisation factor between the Bessel functionand the spectrum I (cid:48) ( x M ). Following their suggested modification,we replace K (1 / Θ e ) by 2 Θ e , thereby assuming the same high-temperature limit.We numerically compute ν c in Eq. (8) by equating the opti-cally thin and thick volume emissivities at ν c , (cid:15) s ( ν c ) = e c √ π ν c ( n e + n + ) K (1 / Θ e ) I (cid:48) ( x M ) = π k B T e H T c ν c . (12) MNRAS000 , 1–16 (2020) pectra of Sgr A ∗ from 3D radiative GRMHD x ( r g ) z ( r g ) &