Spectral and imaging properties of Sgr A* from high-resolution 3D GRMHD simulations with radiative cooling

ABSTRACT

1 INTRODUCTION

It is widely believed that most galaxies harbour supermassive black holes (SMBHs) in their galactic centres, with masses ranging from millions to billions of solar masses. Over the past few decades, the black hole (BH) candidate in the centre of the Milky Way, Sagittarius A* (hereafter Sgr A*), has proven an exceptional laboratory for studies of accretion and outflow physics of BHs due to its proximity to Earth. A significant effort has been invested in determining the BH mass and distance for Sgr A* (e.g. Reid 1993; Schödel et al. 2002; Bower et al. 2004; Ghez et al. 2008; Gillessen et al. 2009, 2017; Boehle et al. 2016; Gravity Collaboration 2018; Reid et al. 2019). We adopt the current best-fitting BH mass of $4.1\pm 0.03\times {10}^6{M}_{\odot }$⁠, which was measured by the orbital motion of stars and gas clouds (Gillessen et al. 2017; Gravity Collaboration 2018), and the distance of 8.15 ± 0.15 kpc, which was obtained from trigonometric parallaxes and proper motions of massive stars around Sgr A* (Reid et al. 2019). Given the mass and the distance, the angular size of the Schwarzschild radius, r_S = 2GM_BH/c², is $\approx 10\mu \mathrm{a}\mathrm{s}$⁠, which subtends a larger area in the sky than any other known BH, including all stellar-mass BHs.

Recently, mounting attention has been paid to the study of Sgr A* with the advent of the pioneering instruments gravity  (Gravity Collaboration 2017, 2018) and the Event Horizon Telescope (EHT ; Doeleman et al. 2009; Event Horizon Telescope Collaboration 2019), capable of probing 10–30 $\mu \mathrm{a}\mathrm{s}$ scales. These instruments allow us to profoundly improve our understanding of the physical processes associated with the accretion and relativistic jet formation in the immediate vicinity of SMBHs and thus demand equal measures of theoretical support and precise modelling of the spectrum generated by the radiation from the accretion flow around Sgr A*.

The mass accretion rate around Sgr A* (in units of $M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠) is estimated to be in the range of $\sim {10}^{-9}<\dot{M}<{10}^{-7}$⁠, as constrained by the measured Faraday rotation measure at mm/sub-mm wavelengths (Aitken et al. 2000; Bower et al. 2003; Marrone et al. 2007). Such a low accretion rate favours hot accretion flow models for the accretion disc instead of the radiatively efficient, thin disc models (Shakura & Sunyaev 1973). Many theoretical scenarios have been invoked and excluded to account for the nature of accretion and outflows in the hot accretion flow: the standard advection-dominated accretion flow (ADAF; Narayan & Yi 1994; Narayan, Yi & Mahadevan 1995; Narayan et al. 1998) and the Bondi–Hoyle models are ruled out, since these models are expected to yield an accretion rate of $\sim {10}^{-5}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠, which is two orders of magnitude higher than the measured upper limit (Bower et al. 2003). Yuan, Quataert & Narayan (2003) re-examined the radiatively inefficient accretion flow (RIAF) model for the spectrum of Sgr A* and argued that the presence of outflows within the Bondi radius plays a vital role in reducing the mass accretion rate. Alternatively, the convection-dominated accretion flow (CDAF) and the jet-dominated models are capable of reproducing the spectrum that is consistent with the observed accretion rate (Falcke & Markoff 2000; Quataert & Gruzinov 2000; Markoff, Bower & Falcke 2007). The existence of 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 confirmed that $\gtrsim 99\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}$ of the matter initially captured by the SMBH is lost before it reaches the innermost region around Sgr A*, which is consistent with the CDAF model or adiabatic inflow–outflow solution (ADIOS; see Blandford & Begelman 1999; Yuan & Narayan 2014 for the detailed model descriptions) model.

Although semi-analytic models provide an important framework for understanding the nature of the accretion flow around Sgr A*, numerical simulations are required to capture the time-dependent, turbulent evolution of the accretion flow. In particular, a self-consistent magnetohydrodynamical (MHD) description enables us to demonstrate accretion processes induced by the magnetorotational instability (MRI; Balbus & Hawley 1991) without imposing an arbitrary anomalous viscosity. In earlier numerical studies, 3D pseudo-Newtonian MHD simulations were carried out to model the synchrotron radiation from accretion flows (e.g. Goldston, Quataert & Igumenshchev 2005; Ohsuga, Kato & Mineshige 2005; Chan et al. 2009; Huang et al. 2009). However, the non-relativistic treatment in the simulations has disadvantages for modelling the synchrotron radiation, since it is mainly emitted in the immediate vicinity of the central BH, where relativistic effects cannot be ignored: shocks develop differently for the relativistic plasma when subject to intense magnetic and gravitational fields (Del Zanna, Bucciantini & Londrillo 2003), and the curvature of space–time becomes significant. Several other works made use of general–relativistic magnetohydrodynamical (GRMHD) simulations for studying the dynamics and spectral properties of Sgr A* in two dimensions (2.5D; e.g. Noble et al. 2007; Mościbrodzka et al. 2009; Hilburn et al. 2010; Dibi et al. 2012; Drappeau et al. 2013; Mościbrodzka & Falcke 2013), or in three dimensions (3D) (e.g. Dexter, Agol & Fragile 2009; Dexter et al. 2010; Dolence et al. 2012; Shcherbakov, Penna & McKinney 2012; Dexter & Fragile 2013; Mościbrodzka et al. 2014; Chael et al. 2018; Davelaar et al. 2018). In general, 2.5D simulations are a reasonable and computationally cheaper option to conduct a parameter study for reproducing the spectrum of Sgr A*, but it is known that simulations with axisymmetric coordinates cannot sustain MRI-driven turbulence, which decays over the local orbital time as a consequence of Cowling’s antidynamo theorem (Cowling 1933, see also Hide & Palmer 1982 for the generalized description). Therefore, full 3D GRMHD simulations are necessary to perform detailed studies of the nature of the accretion flows around Sgr A* and the emitted spectrum. For instance, it was confirmed that thick accretion discs are able to generate and advect large-scale poloidal magnetic flux through dynamo action when resolved properly (Liska, Tchekhovskoy & Quataert 2018a).

The bolometric luminosity of Sgr A* is extremely low, $L_{\mathrm{b}\mathrm{o}\mathrm{l}}\sim {10}^{36}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}\approx {10}^{-9}{L}_{\mathrm{E}\mathrm{d}\mathrm{d}}$⁠, where L_Edd is the Eddington luminosity. Given such a low luminosity, it has been thought that the radiative cooling losses of Sgr A* are negligible, since the losses are likely not strong enough to have a significant impact on the dynamics of the accretion flow. Based on this argument, all previous works with full 3D GRMHD simulations have ignored the radiative cooling losses for simplicity. Although this assumption may be reasonable, Dibi et al. (2012) argued based on their 2.5D simulations that cooling losses play an increasingly important role for higher accretion rates and possibly alter the dynamics and resulting spectra of Sgr A*, even within the allowed range of accretion rates based on polarization and X-ray studies. One potential impact is that the radiative cooling reduces the gas pressure and the disc vertical scale height, resulting in a decrease in turbulence and a more ordered magnetic field (Fragile & Meier 2009). Moreover, many questions remain unanswered: how do radiative cooling losses affect the turbulence features of the disc, and thus the angular momentum transfer of the accreting plasma? How does radiative cooling together with GR effects results in the observed spectra from Sgr A*? Is radiative cooling indeed negligible for the mass accretion rate range of Sgr A*? Even though the effects of radiative cooling can be minor for the case of Sgr A*, the quantitative evaluation of cooling effects is highly demanded because it must play a greater role for SMBHs with higher mass accretion rates, such as M87.

In this paper, we perform the first full 3D GRMHD simulations, which include cooling losses via bremsstrahlung, thermal synchrotron emission, and inverse Compton scattering. Due to the significant computational expense of the full 3D simulations, we cannot explore the full range of various parameters (e.g. spin, magnetic configuration, electron distribution function, misaligned disc). Instead, we use parameters compatible with earlier studies, assuming a rapidly rotating BH, weak poloidal initial magnetic field, and a fixed temperature ratio between ions and electrons of T_i/T_e = 3 (we also carry out additional simulations with different electron temperature prescriptions for comparison). We then examine the effect of radiative cooling on the dynamics of the accretion flow and the resulting spectra and images for different accretion rates within the allowed range.

This paper is structured as follows. In Section 2, we give a technical description of the numerical methods used, including the simulation setup, and the treatment of radiative cooling losses. In Section 3, we describe the results of how cooling losses play a role in changing the dynamical evolution of accreting matter. In Section 4, we discuss the best-bet model for Sgr A*, the effects of cooling on the resulting spectra and sub-mm images, and the variability of multiwavelength spectra. We also compare our 3D work to previous 2.5D work. We summarize our results in Section 5.

2 TECHNICAL DESCRIPTION OF METHOD

All simulations are performed with the H-AMR code (Liska et al. 2019a; Porth et al. 2019), which branched off HARM2D (Gammie, McKinney & Tóth 2003; Noble et al. 2006) in its early days. It is accelerated by graphical processing units and improved with a staggered grid for constrained transport of magnetic fields (Gardiner & Stone 2005) to preserve ∇ · B = 0, more robust inversion (Newman & Hamlin 2014) adaptive mesh refinement (AMR, not utilized in this work), static mesh refinement, and a locally adaptive time-step (see Chatterjee et al. 2019, Appendix  A). It adopts a piece-wise parabolic method (Colella & Woodward 1984) for reconstruction of cell-centred quantities at cell faces, which is third-order accurate, for the spatial reconstruction at cell faces from cell centres, and a second-order time-stepping.

The broad-band spectrum is calculated from the GRMHD output, using the general–relativistic Monte Carlo scheme grmonty (Dolence et al. 2009), which includes synchrotron emission and absorption, and inverse Compton scattering for a relativistic thermal Maxwell–Jüttner distribution of electrons. Technically, grmonty  cannot produce synthetic images but only spectra. Thus, we ray-trace the GRMHD-produced spectra by integrating the general–relativistic radiative transfer (GRRT) equations using the bhoss code (Younsi, Wu & Fuerst 2012; Younsi et al. 2016, 2020) to generate synthetic images at 230 GHz that can help us to infer the expected images of Sgr A* from the upcoming eht  project. In bhoss , the calculation of radiative processes includes synchrotron emission and absorption only, which is sufficient for imaging at the EHT frequency of 230 GHz. Since the sub-mm regime of the spectra is dominated by synchrotron emission, which both codes calculate, we verify consistency in our spectral calculations from both codes by comparing the resulting spectra in the radio to near-infrared (NIR) bands (see Appendix  A).

2.1 Numerical grid and floors

For convenience, we adopt Heaviside–Lorentz units, which absorb a factor of $\sqrt{4\pi }$ for the magnetic field 4-vector, b^μ, so that the magnetic pressure is P_B ≡ b²/2. Furthermore, the typical natural units are used, GM = c = 1, which sets the length unit to be the gravitational radius, $r_{\mathrm{g}}\equiv \mathrm{G}\mathrm{M}/c^2$⁠, and the time unit to be the light crossing time, $t_{\mathrm{g}}\equiv \mathrm{G}\mathrm{M}/c^3$⁠, where $G,M,c$ are the gravitational constant, BH mass, and the speed of light, respectively. We use a spherical–polar axisymmetric computational grid (⁠$r,\theta ,\varphi )$ extending from 0.85 r_H to 250 r_g, where the event horizon radius $r_{\mathrm{H}}\equiv r_{\mathrm{g}}(1+\sqrt{1-a_{\star }^2})$⁠. Here, we set the dimensionless BH spin parameter to $a_{\star }\equiv cJ/{\mathrm{G}\mathrm{M}}^2$ in a Kerr–Schild foliation, where $J=r_{\mathrm{H}}{c}^2$ is the angular momentum at the event horizon. The grid is uniformly spaced with respect to a set of internal coordinates $(x^1,x^2,x^3)$⁠, which can be converted to (⁠$r,\theta ,\varphi )$⁠, respectively¹. This conversion leads to a logarithmic spacing in r such that the cells have a higher resolution for smaller values of r. The spatial resolution near the event horizon is $\mathrm{\Delta }r\approx 0.01r_{\mathrm{g}}$ for the model with the highest resolution. To prevent the aspect ratio of the cells from becoming too large near the polar singularity, we reduce the resolution in ϕ direction gradually towards both poles. We use outflow boundary conditions for both inner and outer radial boundaries and reflecting boundary conditions in the θ direction. Note that the inner boundary is causally disconnected from the flow, as it is located within the event horizon.

It is common for GRMHD simulations to crash if either the density or the internal energy becomes very low, particularly in the funnel region along the polar axes or near the outer radial boundaries. To avoid this, we apply numerical floors for the density and the internal energy (see Appendix B3 of Ressler et al. 2017 for more detailed discussions): a minimum rest mess density, ${\rho }_{\mathrm{f}\mathrm{l}}=\mathrm{m}\mathrm{a}\mathrm{x}[b^2/20,u_g/150,{10}^{-6}(r/r_{\mathrm{g}}{)}^{-2}]$⁠, and a minimum internal energy density, $u_{g,\mathrm{f}\mathrm{l}}=\mathrm{m}\mathrm{a}\mathrm{x}[b^2/750,{10}^{-7}(r/r_{\mathrm{g}}{)}^{-2\mathrm{\Gamma }}]$⁠, where b and u_g are the co-moving magnetic field strength and the internal energy density, respectively. We normalize the mass density such that the maximum density is ρ_max = 1.

2.2 Simulation models

2.3 Radiative cooling

We take into account the radiative cooling self-consistently in our calculation of the gas temperature by including the effects of bremsstrahlung, synchrotron, and the inverse Compton losses. We adopt the equations of Esin et al. (1996) for computing the radiative cooling losses. These formulae have been implemented and tested in previous numerical studies of Sgr A* (Fragile & Meier 2009; Dibi et al. 2012; Straub et al. 2012; Drappeau et al. 2013).

3 RESULTS

We initialize our fiducial model by following the ‘best-bet’ model, that is widely agreed by previous 2.5D parameter surveys (e.g. Mościbrodzka et al. 2009; Dibi et al. 2012; Drappeau et al. 2013), which have T_i/T_e = 3 and a_⋆ = 0.9375 (see Table 1).

3.1 General evolution

Our simulations start with the initial torus in hydrostatic equilibrium (Fishbone & Moncrief 1976). As the turbulence triggered by tangled magnetic fields transports angular momentum outward, the gas flows towards the central BH generating a thick disc, akin to an RIAF. The cooled models require a specific density unit to be pre-set in order to achieve the designated mass accretion rate in the simulations (see Table 1). However, the non-cooled models are scale-free and therefore scaled by the corresponding density unit (in GRRT post-processing) to enable comparison with the cooled models.

Fig. 1 shows the density contour map overlaid with the magnetic field structure, at which the data are averaged for 5000–6000 t_g.² The overall evolution of the accreting hot accretion flow is similar between the non-cooled and cooled runs when the mass accretion rate is smaller than ${10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠. However, the effect of cooling becomes increasingly important and shows visible differences in model C3D10RM, where the target mass accretion rate is ${10}^{-7}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ (the estimated value is $8.22\times {10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠): when radiative cooling is on, the density increases significantly in the mid-plane and the magnetic field within the disc is less turbulent. This is because cooling reduces the gas pressure and the corresponding scale height of the accretion flow, thus increasing the dominance of magnetic fields: the plasma beta, β_mag ≡ P_g/P_B, decreases due to the reduced gas pressure and the compressed volume. Such highly magnetized plasma tends to be stable against the MRI, and thus the magnetic field within the disc becomes less turbulent.

Figure 1.

Comparison of the time-averaged density contour map between the NC3RM (left, non-cooled) and C3D10RM (right, cooled) simulations. The non-cooled model is scaled by multiplying the same density unit (see ρ_scale in Table 1) for the corresponding cooled model. Given the density unit, the average mass accretion rates are $8.24\times {10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ and $8.22\times {10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ for the (scaled) non-cooled and cooled models, respectively. The light black contours represent the magnetic field, and the dashed magenta lines denote the jet boundary, defined as where the magnetization parameter is σ = 1. All variables are averaged over 5000 t_g–6000 t_g.

Open in new tabDownload slide

Figure 2.

The averaged density profile along the disc overs the time interval 6000 t_g–8000 t_g (see equation 13). The solid lines represent the result from models C3D1RM, C3D10RM, and C3D100RM, and the dashed lines are the profile from the non-cooled model NC3RM with re-scaling to each of the cooled models. The vertical dotted line indicates the location of the event horizon, $r_{\mathrm{H}}=r_{\mathrm{g}}(1+\sqrt{1-a_{\star }^2})=1.35r_{\mathrm{g}}$⁠.

Open in new tabDownload slide

Radiation processes occur predominantly in hot accretion flows. To obtain insight into the physical properties of the flows, we show contour maps slices of density, B², and electron temperature at a single time-step (⁠$t=5000t_{\mathrm{g}}$⁠) for the highest resolution model (Fig. 3). For the hot accreting plasma, it is evident that the typical density is $\sim {10}^{-16}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}$⁠, which corresponds to $N_e\approx 6\times {10}^7\mathrm{c}{\mathrm{m}}^{-3}$ in fully ionized plasma, where N_e is the electron number density. The typical strength of magnetic fields is B$\simeq 30\mathrm{G}$ near the event horizon, which decreases with increasing distance from the BH. The electron temperature is maximal within the ‘funnel’ over the pole, which is up to ∼10¹³ K, and is ∼10¹² K in the mid-plane. Despite the high temperature in the funnel, it typically produces negligible emission as a consequence of the extremely low densities in this region. Note that since we assume a relativistic thermal Maxwell–Jüttner distribution for the radiative processes, the question still remains as to which fraction of non-thermal electrons can be generated within the plasma by steepening of MHD waves or inducing shocks/turbulence through mechanisms akin to magnetic reconnection, and thus how it contributes to the radio emission (Yuan et al. 2003; Davelaar et al. 2018).

Figure 3.

Colour contour maps of the mass density, the magnetic field strength, and the electron temperature for the model C3D1RH. The figure presents a single time slice at 5000 t_g.

Open in new tabDownload slide

Given the electron temperature range of 10¹¹–10¹² K in the accretion flow, as seen in Fig. 4, bremsstrahlung cooling is relatively weak. Note that in our simulations, the Comptonization of bremsstrahlung is neglected as it is never of importance compared to synchrotron emission over the temperature range of interest. On the other hand, synchrotron cooling with Compton enhancement is dominant in the mid-plane near the BH: optically thin synchrotron radiation at $r\lesssim 5r_{\mathrm{g}}$ is responsible for the sub-mm peak in the spectral energy distributions (SEDs), which lies within 10¹¹ Hz < ν < 10¹⁴ Hz (see S4.1.1). Inverse-Compton scattering of synchrotron photons is active at $r\lesssim 8r_{\mathrm{g}}$⁠. The emission region and the relative strength of radiative processes are consistent with the previous study (Moscibrodzka et al. 2007), which computed the radiative properties around a hot accretion disc considering the balance of cooling and heating. The mean electron temperature in the accretion flow at $r\lesssim 8r_{\mathrm{g}}$ is 〈T_e〉 ≈ 3 × 10¹¹ K, and therefore the average increase of energy in a single scattering can be approximated as $A=1+4{\mathrm{\Theta }}_{\mathrm{e}}+16{\mathrm{\Theta }}_{\mathrm{e}}^2\approx 4\times {10}^4$ (see Esin et al. 1996). As a result, the frequencies of these scattered photons are shifted to the range of 0.02–20 keV, which is consistent with the observed emission in X-rays (0.5–8 keV) (Baganoff et al. 2001, 2003). The quiescent X-ray emission of Sgr A* is extended, with an intrinsic size of ∼1.4^″ (Baganoff et al. 2003), which is coincident with the Bondi accretion radius calculated from the measured BH mass and ambient temperature (Yuan et al. 2003). It is known that ∼90 per cent of the total X-ray emission originates from the outer part of the disc (Neilsen et al. 2013) and is dominated by bremsstrahlung, which is beyond the scope of this work: we calculate the emission within r = 20 r_g, where the X-ray emission is predominantly produced by the synchrotron self-Compton process. We will further discuss the spectral properties of the X-ray emission in Section 4.1.

Figure 4.

Colour contour maps of the radiative cooling rates: bremsstrahlung (left), synchrotron (middle), and inverse Compton scattering (right) for the model C3D1RH. The figure presents a single time slice at 5000 t_g. Radiative cooling is calculated only when B²/ρ < 1 (see the white area, where the cooling process is off).

Open in new tabDownload slide

3.2 Mass accretion rates

The mass accretion rate is a critical factor in determining the radiation fluxes. In previous works that omit radiative cooling, the simulations are scale-free and must therefore be scaled with an arbitrary density unit during the GRRT post-processing to achieve the desired accretion rate. The scaled variables are used to calculate synthetic spectra, which match with observations (e.g. Dexter et al. 2010; Shcherbakov et al. 2012; Mościbrodzka et al. 2014). Our radiatively cooled simulations, however, are not scale-free: the calculation of cooling rates requires specifying variables in physical units.

Figure 5.

The mass accretion rate at the event horizon as a function of time. Upper panel: The solid and dashed lines represent the results from models with radiative cooling and models without radiative cooling (NC3RM), respectively. The accretion rate from the NC3RM model is scaled by the density unit that corresponds to the each cooled counterpart. Lower panel: The accretion rate from the simulations at the different resolutions: C3D1RH (high, dot–dashed), C3D1RM (intermediate, solid), and C3D1RL (low, dotted).

Open in new tabDownload slide

In the lower panel of Fig. 5, we show how the angular momentum transport through the turbulence of the accretion flow can be affected by resolution effects. To examine if the MRI is resolved properly, we calculate the ‘MRI quality factors’ (i.e. Q-factors), which are defined as the number of cells available for resolving the fastest-growing MRI mode in each direction. The Q-factor of the lowest resolution case (C3D1RL) is ∼3, which is below the nominal Q value of 10–20 for capturing the saturation level of the MRI (Hawley, Guan & Krolik 2011). Thus, it is obvious that the mass accretion rate in the model C3D1RL drops significantly after 2000 r_g since the low-resolution run fails to resolve MRI-driven turbulence. However, for the intermediate (C3D1RM) and highest (C3D1RH) resolution cases, the Q-factors are 12 and 20, respectively, which are large enough to sustain the MRI-driven turbulence. These Q-factors are indicative of a criterion above which the simulations satisfactorily resolve the MRI but cannot be used for the analysis of turbulent features in the flow. We will further discuss the disc properties in Section 3.3.

We carry out multiple simulations with different density units to target Sgr A* mass accretion rates (in units of $M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠) of ${10}^{-9},{10}^{-8},{10}^{-7},\mathrm{a}\mathrm{n}\mathrm{d}{10}^{-6}$⁠. The largest accretion rate among these target values is beyond the observed range around Sgr A* (⁠$2\times {10}^{-9}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}<\dot{M}<2\times {10}^{-7}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠), but this model is included to compare the results from other simulations with the case of an extremely high accretion rate. In the upper panel of Fig. 5, the solid and dashed lines represent the mass accretion rates, which are calculated from the cooled and non-cooled models, respectively. The non-cooled model is re-scaled by the same density unit for each cooled model. For the model with strong cooling (C3D100RM), it is clear that the overall accretion rate is smaller by a factor of 2, compared to the non-cooled case with the same density unit. However, the models with weak cooling (C3D01RM and C3D1RM) show no significant differences in the accretion rate between the cooled and non-cooled models, which is surprising since the cooled model is expected to lose less angular momentum compared to the non-cooled model. For model C3D10RM, the accretion rate is even slightly higher than in the non-cooled model. The reason for this is the enhanced density in the vicinity of the event horizon playing a role in increasing the accretion rate (see Fig. 2), which compensates for the weak angular momentum transport.

3.3 Disc properties

Figure 6.

Scale height profiles (H/R) according to equation (15) for models with different density scales and thus different accretion rates. The solid thick lines represent the mean value over a time interval 6000 t_g–8000 t_g, and the shaded regions represent the variations during this interval.

Open in new tabDownload slide

Figure 7.

Barycentric radii of the disc [i.e. density-weighted radii, see equation (16)]. The solid curves represent models with medium resolution but differing density scales, and the dot–dashed line and the dashed line represent models that are the same as the model C3D1RM but with lower (C3D1RL) and higher (C3D1RH) resolutions.

Open in new tabDownload slide

4 DISCUSSION

4.1 Radiative properties of Sgr A*

4.1.1 Spectral energy distribution

Once we consider radiative cooling, we can no longer scale the GRMHD data to fit the observed flux. Hence, we choose the best-fitting data that produce the flux lying within the observed ranges at 230 GHz (3 Jy–4.2 Jy; see observations compiled within Connors et al. 2017), but note that this is not a statistical fit. To compare the results between the cooled data and the non-cooled data, we scale the latter with the same mass density unit for each of their cooled counterparts. As seen in Fig. 8, with our fixed value of a_⋆ = 9375 and T_i/T_e = 3, the model C3D01RM, at which the target mass accretion rate is $\dot{M}={10}^{-9}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠, is reasonably consistent with the observations. The overall shapes of the light curves between the cooled (C3D01RM) and non-cooled (NC3RM) models are similar to each other; however, the average fluxes in the non-cooled model are slightly higher than in the cooled model. The average fluxes at 230 GHz, which are calculated for the time after 5000 t_g, are 3.63 ± 0.41 Jy and 3.8 ± 0.35 Jy for the cooled and non-cooled models, respectively.

Figure 8.

The synchrotron light curves at 230 GHz for the cooled model (C3D01RM; blue colour) and the non-cooled model (NC3RM; orange colour), which are calculated using bhoss (Younsi et al. 2020; Younsi et al., in preparation). The grey shaded region represents the flux range consistent with observations compiled by Connors et al. (2017). The synthetic spectra (Fig. 9) are calculated after 5500 t_g (vertical dotted line), where the synchrotron flux lies within the observed range for both the cooled and non-cooled models.

Open in new tabDownload slide

We calculate the spectra from the GRMHD simulation data using the Monte Carlo radiative transport code grmonty (Dolence et al. 2009), which computes synchrotron emission and absorption, and inverse Compton scattering in full general relativity. Fig. 9 shows the SEDs for the cooled model (C3D01RM) and the non-cooled model (NC3RM) with the same density scale. The SEDs have two peaks: the sub-mm peak and the far-UV peak. Thermal synchrotron emission from mildly relativistic electrons is responsible for the sub-mm peak and these same photons are then Compton upscattered to produce the far-UV peak. To check the dependency of the viewing angle, we set the number of θ bins to 6, within which the fluxes are averaged to represent the values for the range of the inclination angle between the BH spin axis and the observer line of sight. In general, the fluxes slightly increase with increasing inclination angle. This is mainly due to the orbiting plasma that is approaching the observer and the emission being more strongly Doppler boosted at higher inclination angles (i.e. close to edge-on).

Figure 9.

Spectral energy distributions of Sgr A*, which are calculated from the simulation results of C3D01RM (cooled, blue) and NC3RM (non-cooled, red). The GRMHD data are averaged over the time interval t_g 6000–10000 t_g (see Fig. 8 for the light curve within this interval). The solid lines represent the mean value of the spectrum, and the shaded regions represent the variation of the spectrum during the time period. Observational data points are taken from Melia & Falcke (2001) and Schödel et al. (2011) at the upper limit of near-to-mid IR band, Connors et al. (2017) in the sub-mm band, Bower et al. (2019) at terahertz frequencies (233, 678, and 870 GHz), and Baganoff et al. (2001) and Baganoff et al. (2003) for X-rays (2–10 keV). The X-ray flux in the ray-traced GRMHD data should be below ∼10 per cent of the observed quiescent (lower) X-ray flux since most of the X-rays should be emitted from the outer disc via bremsstrahlung (Neilsen et al. 2013; Wang et al. 2013), which is not included in this calculation. The different panels represent the results with different inclination angles: 45^○–60^○ (left), 60^○–75^○ (middle), and 75^○–90^○ (right).

Open in new tabDownload slide

The SED in the cooled model (C3D01RM) differs slightly from the non-cooled model (NC3RM): for the cooled model, the overall flux, including the peak value at the sub-mm bump, is slightly lower than for the non-cooled model. The NIR emission originates from the innermost regions (⁠$2r_{\mathrm{g}}<r<6r_{\mathrm{g}}$⁠; see Mościbrodzka et al. 2009), where the gas temperature and the magnetic field intensity are high. The relatively weak NIR emission in the cooled model is indicative of the lower gas temperature due to the inclusion of radiative cooling. The peak of the far-UV flux is also slightly higher in the non-cooled model than in the cooled model, mainly due to the higher flux of the seed photons over the NIR band, and the peak frequency in the non-cooled model is ∼3.5 times larger than in the cooled model. This is because the average increase of energy in the scattering is formulated to $A=1+4{\mathrm{\Theta }}_{\mathrm{e}}+16{\mathrm{\Theta }}_{\mathrm{e}}^2$ (Esin et al. 1996), implying that the higher temperature in the non-cooled model leads to upscattering of photons into the higher energy range. We note that these differences in the SED arise from our adoption of identical density unit values in the GRRT post-processing for both the cooled and non-cooled models. While the resulting SEDs lie within the observational constraint at 230 GHz (see Fig. 9), if the non-cooled model is normalized in the GRRT calculation by decreasing the density unit to match the flux of the cooled model at 230 GHz, the differences become less significant. The adjusted density unit to achieve this matching of the 230-GHz fluxes is  0.84 times (i.e. smaller than) the cooled model’s density unit value, and thus the estimated accretion rate is also smaller in the non-cooled model.

Figure 10.

Spectral energy distribution of Sgr A*, calculated using GRMHD results with different electron temperature prescriptions: T_i/T_e = 3 (C3D01RM, green line), T_i/T_e = 20 (C3D1RMFT20, red line), and T_i/T_e depending on the plasma magnetization (C3D1RMRh20, blue line). The inclination angle is 60^○–75^○.

Open in new tabDownload slide

4.1.2 Synthetic images at 230 GHz

To compute the synthetic mm (230 GHz) images, we use the GRRT code bhoss  (Younsi et al. 2012, 2016; Younsi et al. 2020; Younsi et al., in preparation), in which the radiation processes include synchrotron emission and absorption. The inclination angle of both the BH spin axis and the disc normal to the line of sight is currently poorly constrained by observations. For example, in kinematic studies of the star S2 near the Galactic Centre, the inclination angle is best fitted to 134°, which is moderate Gravity Collaboration (2018). One may also expect the system to be nearly edge-on, given that highly inclined sources often produce linear polarization in compact radio sources around Sgr A* (Bower et al. 2003). In this work, we set inclination angles of ${50}^{\circ },{70}^{\circ },\mathrm{a}\mathrm{n}\mathrm{d}{90}^{\circ }$ to clearly see how the inclination angle changes the shape of spectra and images. Note that we assume an aligned jet with the angular momentum of accreting gas. In fact, misalignment is likely to occur in Sgr A* since infalling gas cannot be quickly torqued into alignment with the BH spin, given the geometrically thick disc (Dexter & Fragile 2013; Liska et al. 2018b; White et al. 2020).³ We will present the effects of tilted accretion discs in a different work (Chatterjee et al. 2020). Fig. 11 shows the synthetic images at 230 GHz for a single time snapshot at $8000t_{\mathrm{g}}$⁠, which shows detailed turbulent substructure. The ring-like structure is produced by the gravitational lensing effect that magnifies the emissions from accretion flow, and the bright patches in the left side are the result of Doppler-boosted emissions from the approaching side of the disc approaching the observer. These patches are brighter in the non-cooled model than in the cooled model due to the relatively higher gas temperature near the BH.

Figure 11.

GRRT synthetic images at 230 GHz for the cooled model (C3D01RM, top row) and the non-cooled model (NC3RM, bottom row). They are taken from a single snapshot at 8000 t_g. The columns represent the result with different inclination angles: 50^○ (left), 70^○ (middle), and 90^○ (edge-on, right). The images are post-processed with the GRRT code, bhoss (Younsi et al., in preparation).

Open in new tabDownload slide

In Fig. 12, we take the blurred images, which are averaged over the time interval between 8000 and $9500t_{\mathrm{g}}$⁠, as an appropriate proxy for the eht  image. The time interval of $1500t_{\mathrm{g}}$ corresponds to ∼8.4 h, given the BH mass is $\sim 4.1\times {10}^6{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ (Gravity Collaboration 2018). The blurred images were obtained through convolution with a Gaussian filter, at which the full width at half-maximum (FWHM) is 20 μas. As seen in the figure, the emission is dominated by the left side of the disc, which has a symmetric crescent shape for all models. Such a crescent shape is the blurred region of the aforementioned hot patches, which are produced by Doppler beaming. The intensity contrast within the crescent in the non-cooled model is larger compared to the cooled model, as expected from the higher temperature in the regions corresponding to the hot patches. In general, the BH shadow is clearly visible for lower inclination angles (θ₀ < 70°) but becomes less visible for the edge-on images.

Figure 12.

Time-averaged and blurred images at 230 GHz for the cooled model (C3D01RM, top row) and the non-cooled model (NC3RM, bottom row). The synthetic images are averaged over 8000 t_g–9500 t_g and blurred by the convolution with a 2D Gaussian filter, at which the full width at half-maximum (FWHM) is 20 μas (white circle in bottom right of each panel). The size of the images extends to ±65 μas, which corresponds to ∼15 r_g. The columns represent the result with different inclination angles: 50° (left), 70° (centre), and 90° (edge-on, right).

Open in new tabDownload slide

4.2 Variability

The dynamical environment around Sgr A* drives flares through various mechanisms: sudden electron heating by magnetic re-connection, star–disc interactions, stochastic acceleration, gravitational lensing of ‘hotspots’ in the accretion flow, and sudden increases in the mass accretion rate due to the infall of clumps of material (e.g. Markoff et al. 2001; Nayakshin, Cuadra & Sunyaev 2004; Yuan, Quataert & Narayan 2004; Trippe et al. 2007; Hamaus et al. 2009; Gravity Collaboration 2018).

As seen in Fig. 13, both the cooled (C3D1RMRh20) and non-cooled (NC3RM) models are highly variable in the different multiwavelength bands. The NIR light curve behaves similarly to the X-ray light curve, and the eruption events in both wavebands are roughly correlated with a pronounced rise in the mass accretion rate. Fig. 14 shows the 230-GHz images for the cooled model C3D1RMRh20 during the NIR quiescent state (orange dotted line in Fig. 13) and a flaring state (blue dotted line). In this figure, the optical depth of the accretion flow increases during the flaring event, in accordance with the increase in the mass accretion rate. Fig. 15 shows the corresponding SEDs and clearly illustrates the overall rise in flux across all frequencies. However, the X-ray flux is lower than the quiescent limit even during the flaring event, while the NIR flux exceeds the quiescent flux by a factor of ∼4. Hence, an increase in the accretion rate can probably trigger an NIR flare without a clearly detectable X-ray flare, which can account for why some NIR flaring events do not exhibit simultaneous X-ray flares (Hornstein et al. 2007). This can be seen for both the cooled and non-cooled models, implying that radiative cooling perhaps plays little role in producing flares. However, it is expected that such cooling shortens the duration of the flaring events, which its importance can be significant, in tandem with electron heating (Leichtnam et al., in preparation).

Figure 13.

Light curves at three different frequency bands, from top to bottom: sub-mm (9.4 × 10¹¹ Hz), near-IR (4.5 × 10¹⁴ Hz), and X-ray (5 × 10¹⁷ Hz), for the cooled model (C3D1RMRh20, solid curve) and the non-cooled model (NC3RM, dashed curve). The bottom panel presents the mass accretion rate for the cooled and non-cooled models. The light curves are calculated by using grmonty  with the assumption of D_Sgr = 8.2 kpc, where D_Sgr is the distance to Sgr A*. The fluxes are averaged over the theta bin of 60°–75°. The bottom panel shows the mass accretion rate at the event horizon. The cadence of the simulation is 50 t_g, which corresponds to ∼16 min in Sgr A*. The vertical dotted lines indicate the selected time for a flaring-like event (blue) and a quiescent state (orange).

Open in new tabDownload slide

Figure 14.

Comparison of GRRT synthetic images at 230 GHz for the cooled model (C3D1RMRh20) between two different states: NIR-quiescent state at $t=53h$ and NIR-flaring state at $t=37h$ (see Fig. 13). The inclination angle is 70^○.

Open in new tabDownload slide

Figure 15.

Comparison of SEDs between two different states: NIR-quiescent state at $t=53h$ and NIR-flaring state at $t=37h$ (see Fig. 13). The inclination angle is 70°.

Open in new tabDownload slide

However, the maximum peak of NIR emission is ∼0.5 mJy, which is an order of magnitude lower than the observed flaring flux (Dodds-Eden et al. 2011). Moreover, the X-ray luminosity lies below the value of $2\times {10}^{34}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}$ for the entire time period, which is still identified as an X-ray ’quiescent’ state: the observed luminosities of the bright X-ray flares are $>{10}^{35}\mathrm{e}\mathrm{r}\mathrm{g}{\mathrm{s}}^{-1}$ (Haggard et al. 2019). One possible reason is that our assumption of purely thermal electrons is not sufficient to produce NIR flares, since non-thermal electrons can be produced near the BH via relativistic magnetic re-connection (Werner et al. 2016). The contribution of non-thermal electrons to flaring will be discussed in an upcoming paper (Chatterjee et al., in preparation). As an alternative solution, Dexter et al. (2020) suggested that saturation of magnetic flux can trigger the flaring events in a magnetically arrested disc.

4.3 Comparison with previous works: 2.5D versus 3D

Recent axisymmetric 2.5D GRMHD simulations have explored the effects of radiative cooling on the dynamical evolution of hot accretion flows around Sgr A* (e.g. Fragile & Meier 2009; Dibi et al. 2012; Straub et al. 2012; Drappeau et al. 2013). However, it is known that MRI-driven turbulence is not sustainable in the axisymmetric 2.5D simulations, which decays over the local orbital time as a consequence of the Cowling’s antidynamo theorem (Hide & Palmer 1982). As seen in Fig. 16, the mass accretion rate in the 2.5D run never reaches a steady state and significantly decreases after 1000 t_g due to the lack of angular momentum transport via the MRI, while the mass accretion rate in the 3D run reaches a quasi-stationary state. Therefore, the 2.5D, axisymmetric approximation does not allow running over a long simulation time, and instead it requires choosing the data at a certain period of time before the MRI decays dramatically, or including an artificial magnetic dynamo term in the induction equation (Sądowski et al. 2015).

Figure 16.

Dimensionless mass accretion rate as a function of time for the non-cooled high-resolution 2.5D model (NC2RH, orange) and the non-cooled medium-resolution 3D model (NC3RM, blue).

Open in new tabDownload slide

Nevertheless, in many respects our 3D results agree with the previous 2.5D results (Dibi et al. 2012) in that radiative cooling plays an increasingly significant role with increasing mass accretion rate and its impact becomes important above a mass accretion rate of $\dot{M}>{10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠. The best-fitting Sgr A* model with the constant temperature ratio of T_i/T_e = 3 in our work requires the mass accretion rate of $\sim {10}^{-9}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠, which is similar to the results from the previous 2.5D results (Mościbrodzka et al. 2009; Dibi et al. 2012).

5 CONCLUSIONS

By means of GRMHD simulations and GRRT post-processing, we study the effects of radiative cooling on the dynamics of accretion flows and their resulting spectra. It is generally assumed that radiative cooling is negligible for RIAF discs, which occur at low accretion rates (i.e. $\dot{M}\lesssim {10}^{-7}{\dot{M}}_{\mathrm{E}\mathrm{d}\mathrm{d}}$⁠). However, the importance of radiative cooling increases with the increasing accretion rate. It is poorly understood what the critical value of the BH mass accretion rate beyond which radiative cooling becomes effective actually is, particularly in 3D.

For the calculation of radiative cooling, we adopt the approximate solution for the advection-dominated accretion disc, which includes bremsstrahlung, synchrotron, and inverse Compton scattering (Esin et al. 1996). We assume that the temperature ratio between ions and electrons, T_i/T_e, is constant or depends on the plasma beta. However, recent studies with particle-in-cell simulations show that the temperature ratio increases over time due to the weak ion–electron thermal coupling (Zhdankin, Uzdensky & Kunz 2020). While it remains unclear if there are other mechanisms for the effective energy transfer from ions to electrons (e.g. the ion cyclotron instability; Sironi & Narayan 2015), the subject of changing temperature ratio over time is beyond the scope of the current paper. In this work, full 3D GRMHD simulations with radiative cooling are extended from previous 2.5D simulations studied by Dibi et al. (2012) and Drappeau et al. (2013).

In general, radiative cooling enhances the disc mid-plane density as the gas pressure decreases due to energy loss, which is radiated away. The disc with reduced pressure and compressed volume increases the dominance of magnetic fields, which reduce the angular momentum transport outwards via the MRI. As a result, when radiative cooling is on, the disc structure is different from when cooling is neglected: a density peak appears near the central BH, and the distance of the peak from the BH increases with the strength of radiative cooling (i.e. mass accretion rate). This difference is negligible when the accretion rate is small; however, when the accretion rate is larger than ${10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠, it becomes apparent (see Fig. 2). Since this rate lies within the range of mass accretion rates for Sgr A*, we argue that cooling losses can affect the dynamical evolution to an appreciable degree.

The effects of radiative cooling on the spectra are visible even for the low accretion rate of $\dot{M}<{10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$⁠: cooling reduces the peak flux in the sub-mm bumps due to the decreased gas temperature. The decreased seed photon by synchrotron at the sub-mm bumps results in the decrease of the flux in the X-ray bumps for which inverse Compton is responsible. The synthetic images at 230 GHz, which is calculated by GRRT post-processing, show similar crescent shapes between the cooled and non-cooled GRMHD data but slightly dimmer in the cooled data due to the decreased temperature adjacent to the BH.

Recent studies by Ressler et al. (2020) indicate that the inner regions near the BH could be strongly magnetized, as magnetic fields get advected from stellar winds (see also Ressler, Quataert & Stone 2018), which can lead to the formation of MADs. While it is thought to be inevitable for MADs to produce strong outflows, which are absent in Sgr A* (e.g. Markoff et al. 2007), we plan to conduct a study of MADs with radiative cooling in future work, so as to investigate how cooling would affect the disc and the outflow in such a situation. Another notable caveat of our simulations is the absence of non-thermal electron acceleration, which is deemed to be responsible for X-ray (and perhaps, NIR) flaring in Sgr A* (Neilsen et al. 2013; Ball et al. 2016; Connors et al. 2017). We will discuss the contribution of the non-thermal electron in a different work (Chatterjee et al., in preparation).

ACKNOWLEDGEMENTS

This research was enabled in part by support provided by Oak Ridge Leadership Computing Facility, which is a DOE office of Science User Facility supported under contract DE-AC05-00OR22725, and Calcul Quebec (http://www.calculquebec.ca) and Compute Canada (http://www.computecanada.ca). DY, KC, and SM are supported by the Netherlands Organization for Scientific Research (NWO) VICI grant (no. 639.043.513), ZY is supported by a Leverhulme Trust Early Career Research Fellowship, ML was supported by the NWO Spinoza Prize (PI M.B.M. van der Klis), and AT is supported by Northwestern University and by National Science Foundation grants AST-1815304 and AST-1911080.

DATA AVAILABILITY

The data from the GRMHD simulations and GRRT calculations used in this work are publicly available at https://doi.org/10.5281/zenodo.3988208.

Footnotes

REFERENCES