Inverse Compton Cooling in the Coronae of Simulated Black Hole Accretion Flows
Brooks E. Kinch, Scott C. Noble, Jeremy D. Schnittman, Julian H. Krolik
DDraft version September 7, 2020
Typeset using L A TEX preprint style in AASTeX63
Inverse Compton Cooling in the Coronae of Simulated Black Hole Accretion Flows
Brooks E. Kinch , Scott C. Noble , Jeremy D. Schnittman , and Julian H. Krolik CCS-2: Computational Physics and MethodsLos Alamos National LaboratoryNM 87545, USA Gravitational Astrophysics LaboratoryNASA Goddard Space Flight CenterGreenbelt, MD 20771, USA Department of Physics and AstronomyJohns Hopkins UniversityBaltimore, MD 21218, USA (Revised September 7, 2020)
Submitted to ApJABSTRACTWe present a formulation for a local cooling function to be employed in the diffuse,hot corona region of 3D GRMHD simulations of accreting black holes. This new coolingfunction calculates the cooling rate due to inverse Compton scattering by consideringthe relevant microphysics in each cell in the corona and approximating the radiationenergy density and Compton temperature there by integrating over the thermal seedphoton flux from the disk surface. The method either assumes ion and electron tem-peratures are equal (1T), or calculates them separately (2T) using an instantaneousequilibrium approach predicated on the actual relevant rate equations (Coulomb andCompton). The method is shown to be consistent with a more detailed ray-tracingcalculation where the bulk of the cooling occurs, but is substantially less costly to per-form. As an example, we apply these methods to a harm3d simulation of a 10 M (cid:12) ,non-spinning black hole, accreting at nominally 1% the Eddington value. Both 1Tand 2T approaches lead to increased radiative efficiency and a larger fraction of totalcooling in the corona as compared to the original target-temperature cooling functionused by harm3d , especially in the 1T case. Time-averaged post-processing reveals thatthe continuum spectral observations predicted from these simulations are qualitativelysimilar to actual X-ray binary data, especially so for the 1T approach which yields aharder power-law component (Γ = 2 .
25) compared to the 2T version (Γ = 2 . Corresponding author: Brooks E. [email protected] a r X i v : . [ a s t r o - ph . H E ] S e p Kinch et al.
Keywords:
Magnetohydrodynamical simulations (1966) — General relativity (641) —Accretion (14) — X-ray binary stars (1811) INTRODUCTIONAstrophysical black hole accretion is a profoundly difficult problem to simulate realistically. The richand varied observational signatures from both supermassive and stellar-mass black holes arise fromthe interplay of General Relativity (GR), magnetohydrodynamics (MHD), and radiative processes. Aproper treatment must be at least finely enough resolved in space and time to capture the evolutionof the magneto-rotational instability (Balbus & Hawley 1991; Hawley & Balbus 1991)—the coreunderlying accretion mechanism—on the appropriate background metric (generally, Kerr). The firstcode to achieve this in a global (albeit 2D) context was harm (Gammie et al. 2003), which has formedthe basis of many subsequent codes (Porth et al. 2019) including harm3d (Noble et al. 2009), theglobal 3D GRMHD code we extend and apply in this work. The final component of the physics,radiation, has proved the most difficult to fully incorporate into any code.Simultaneously solving the MHD equations and the global angle- and energy-dependent radiationtransport equation, in General Relativity, is both computationally expensive (typically prohibitivelyso) and technically challenging. Even so, significant progress has been made in the last decade,though the problem is usually made tractable by introducing at least one of the following simplifyingassumptions: abandoning General Relativity in favor of a pseudo-Newtownian description of thegravitational potential, while performing realistic, multi-angle group radiation transport (Jiang et al.2014a,b, 2019b,a); limiting the possible angular-dependence of the radiation field by invoking eitherflux-limited diffusion (Zanotti et al. 2011; Roedig et al. 2012) or, more recently, the “M1 closure”relation, in either axisymmetric (2D) (S¸adowski et al. 2014) or 3D simulations (Fragile et al. 2012,2014; McKinney et al. 2014; S¸adowski et al. 2016); or Monte Carlo (Ryan et al. 2015)/hybrid MCtechniques Ryan & Dolence (2019). Most attempts have eschewed energy-dependent transfer in favorof a “grey” atmosphere—the radiation field is treated as monochromatic, coupled to the fluid onlythrough the Rosseland mean opacity (Rybicki & Lightman 1986). The first of these approximations,the pseudo-Newtonian potential, is especially problematic in regions close to the black hole whereGeneral Relativistic effects play a critical role in determining both the structure of the accretionflow and photon trajectories. The others are essentially variants of a diffusion approximation, andare best suited to the cooler, denser, and optically thick body of the accretion disk, where theenvironment is similar to those found in stellar atmospheres—the field from which these methods,and grey transfer, originate (Chandrasekhar 1960). With the exception of Monte Carlo methods(Ryan et al. 2015; Ryan & Dolence 2019), these are all especially poorly-suited to the diffuse, hot,optically thin corona, especially at small radii near the black hole. harm3d has employed a local cooling function approach to emulate the effects of radiative coolingin moderately accreting (0.01–0.3 Eddington) systems: gravitationally-bound gas hotter than a spec-ified target temperature is cooled to the target temperature over one orbital timescale; the targettemperature is chosen so as to achieve a desired small disk aspect ratio, i.e., a geometrically thindisk. In this paper, we employ a cooling function which instead considers the relevant microphysicswithin each fluid element and from there calculates the rate at which the internal energy of the gasis converted to photons. This calculation requires, in general, knowledge of the energy-dependent nverse Compton Cooling substantial overhead of full transport. We achieve this by employing a series of simplifying butenabling assumptions, which we check against a Monte Carlo radiation transport calculation. Thenew cooling function is “switched on” at a time after the simulation has evolved long enough withthe original, target-temperature cooling function that it has achieved a statistically steady-state, andapplies only in the corona. Within the disk body, the original, target-temperature cooling functionremains in place. In the second part of this paper, we also explore the consequences of weak ion-electron coupling in the corona, in the context of our more physically-motivated cooling function.Throughout, we examine the effects of the new cooling function on the dynamical, thermodynamic,and spectral properties of a fiducial 10 M (cid:12) , non-spinning ( a = 0) black hole simulation accreting atapproximately 1% Eddington. INVERSE COMPTON COOLING FUNCTIONThe target-temperature cooling function was introduced in Noble et al. (2009) and further developedin Noble et al. (2010) and Noble et al. (2011). A radius-dependent target temperature T ∗ (expressedin harm3d ’s dimensionless code units) is set according to T ∗ = π R z ( r ) r (cid:20) H ( r ) r (cid:21) , (1)where R z is the relativistic correction to the vertical component of gravity at radius r and H isthe density-weighted scale height of the disk [equation from Noble et al. (2009), corrected fromAbramowicz et al. (1997)]. The desired disk aspect ratio, H/r , is chosen a priori to achieve ageometrically thin disk—for the starting point ThinHR simulation series used here, the desired diskaspect ratio is
H/r = 0 .
05. When gas on a bound orbit exceeds the target temperature, it is cooledback to the target temperature by introducing a nonzero sink term on the right-hand side of the localstress-energy conservation equation solved by harm3d : ∇ ν T νµ = −L u µ , (2)here T νµ is the stress-energy tensor and u µ is the specific four-momentum; L is chosen so that the gascools to the target temperature over one circular orbital period at its radius. If the gas is at or below T ∗ , or is gravitationally unbound, L = 0. As the gas is cooled, its pressure and therefore supportagainst gravity decreases, settling back toward the midplane and thereby achieving a geometricallythin disk.The target-temperature approach has several key benefits: as intended, it gives rise to a geometri-cally thin, optically thick, relatively cool, dense disk, sandwiched between a hotter, diffuse corona—aconfiguration with considerable observational support (Haardt & Maraschi 1991); the implementa-tion is independent of the central black hole mass scale M and the nominal accretion rate ˙ m (inEddington units) and therefore, in principle, the results of a single simulation can be scaled to bothstellar-mass X-ray binary systems and supermassive active galactic nuclei; and it is easy to evaluateas it depends only on local properties of a given fluid element.On the other hand, of course, it is unphysical: the choice of target temperature is motivated notby the relevant microphysics, but by the desire to achieve a configuration-by-design that agrees well Kinch et al. with observational evidence. There are other concerns as well. By virtue of its implementation as alocal sink term, it is everywhere “optically thin”; that is, the dissipation rate is an entirely intensive quantity—even deep within the disk, energy lost (nominally to photons) simply vanishes, while inreality these photons would diffuse through the optically thick material, scattering and undergoingabsorption/re-emission along the way. In the corona, gravitationally unbound matter does not coolat all, while gas that is cooled does so on a circular orbital timescale which may not relate to itsactual cooling time.An optically thin cooling function is, however, a good approximation in the truly optically thincorona. Thus we seek a more physical cooling function there while still retaining the implementationof Equation 2. In addition, we understand the actual physical mechanism behind coronal cooling:the inverse Compton (IC) scattering of thermal seed photons from the disk surface off of very hotelectrons—this is exactly the physics treated with great care by pandurata (Schnittman & Krolik2013). Below we detail the development of a new cooling function L to replace the target-temperaturecooling function in the corona; to emphasize the physical origin of the new cooling function—andto distinguish it from the target-temperature cooling function which will remain in use in the diskbody—we refer to it simply as the IC cooling function.At first, we will require the ion and electron populations to be at the same temperature in eachsimulation cell, T e = T i , i.e., a one temperature fluid. In section 6, we extend the method totreat the ion and electron temperatures separately, T e (cid:54) = T i , i.e., a two temperature fluid, withthe assumption that electrons are heated only through Coulomb collisions with ions. These twoprescriptions represent, essentially, the limiting cases of maximally- and minimally-efficient radiativecooling, respectively. The one temperature fluid assumption requires no special description of theion-electron coupling mechanism—we simply assume that some strong coupling mechanism exists,or that turbulent dissipation is shared nearly equally between ions and electrons, or a combinationthereof. Later, we will posit a specific coupling mechanism, namely Coulomb collisions, in additionto the assumption that all turbulent energy is initially injected into the ions only.2.1. Inverse Compton Power
To distinguish the corona volume from the disk body, two disk photosphere surfaces, Θ top ( r, φ )and Θ bot ( r, φ ), are defined by integrating the electron scattering opacity from the z -axis toward themidplane: (cid:90) Θ top ( r,φ )0 κρ ( r, φ ) dθ √ g θθ = 1 , (3) − (cid:90) Θ bot ( r,φ ) π κρ ( r, φ ) dθ √ g θθ = 1 , (4)where κ is the electron scattering opacity, 0 . g − . Note that if no solution exists for the aboveequations for a given ( r, φ ) such that Θ top < Θ bot , then simply no disk body exists there. Fora given point ( r, θ, φ ) such that Θ top and Θ bot exist, the point is considered in the disk body ifΘ top < θ < Θ bot ; otherwise, the point is considered in the corona. The location of the photospheresurfaces depends on the choice of accretion rate through the overall scale of the density; as discussedin section 2.4 below, a larger nominal ˙ M / ˙ M Edd implies a larger overall ρ which, all else equal, results nverse Compton Cooling P IC = 43 σ T cγ β u rad , (5)where σ T is the Thomson scattering cross section, c is the speed of light, γ = 1 / (cid:112) − β with β = v/c —where v refers to the electron velocity—and u rad is the radiation energy density. Inthe nonrelativistic limit, i.e., when the dimensionless electron temperature Θ e ≡ k B T e /m e c (cid:28) (cid:104) γ β (cid:105) (cid:39) e (with the averaging indicated by angle brackets performed over a thermal electronvelocity distribution), and so P IC , non − relativistic = 43 σ T c (3Θ e ) u rad . (6)In the relativistic limit, i.e., Θ e (cid:29) (cid:104) γ β (cid:105) (cid:39) e , therefore P IC , relativistic = 43 σ T c (12Θ e ) u rad . (7)Because each expression is much larger than the other in their appropriate regimes, we can representthe IC power in either limit by their sum. In addition, we multiply by the electron density n e = χ ( ρ/m i ), where χ is the free electron fraction (number of free electrons per ion, equal to 1.21 for afully-ionized plasma with solar elemental abundances; a variable χ might also be used to account forthe presence of electron-positron pairs due to pair production) and m i is the mean ion mass ( (cid:39) m p ),to translate from energy exchanged per scatter per unit time to a volumetric cooling rate. The finalexpression is: L IC = 4 σ T cχm i ρu rad Θ e (1 + 4Θ e ) . (8)This is the expression for the IC cooling rate which enters into harm3d ’s stress-energy equation 2in place of the target-temperature cooling rate for those fluid elements in the corona. It requires asinput: the density, the electron temperature, and the radiation energy density. As discussed, we willassume for now that some strong coupling mechanism forces T e = T i . From standard thermodynamicsand the ideal gas law, the gas pressure p gas , the internal energy density u , and the electron and iontemperatures are related by p gas = ( c P /c V − u = n e k B T e + n i k B T i , (9)where c P /c V is the ratio of specific heats (the adiabatic index) [equal to 5/3 for a monatomic gas;though we assume c P /c V = 5 / e = m i m e c P /c V −
11 + χ uρc . (10) Kinch et al.
This expression depends only on the (mass) density and the internal energy density—which, like thedensity, is part of harm3d ’s fluid solution. Again, this equation holds so long as there is some strongcoupling mechanism forcing T e = T i .In order to estimate the radiation energy density u rad at each point in the corona, we first makeseveral key simplifications:1. We ignore general and special relativistic effects. The thermal seed photons launched from thedisk surface are assumed to travel in straight rays, undergoing neither red/blue-shifting norbeaming due to the bulk fluid flow of the rotating accretion disk, nor gravitational redshiftingdue to their origin in a deep gravitational potential.2. We do not account for scattering or obscuration of the disk emission by intervening coronamaterial between the disk surface and a given point in the corona.3. We adopt the “fast light” approximation, i.e., we do not account for the light travel timebetween a point on the disk photosphere and a point in the corona; rather, the radiationenergy density in the corona each time step is computed from the thermal flux from the disksurface at the same time step.With these assumptions in place, we derive an expression for u rad by integrating the thermal seedphoton flux over the disk surface with an appropriate geometric weight. Let r indicate the locationof the coronal cell in question, and let r (cid:48) locate a surface cell on the photosphere. Then: du rad ( r ) = 1 c F disk ( r (cid:48) ) cos ϑdA (cid:48) R , (11)where the factor of 1 /c translates flux into energy density, R = r − r (cid:48) , ϑ is the angle between R and the disk surface normal vector ˆ n , dA (cid:48) is the (infinitesimal) element of the disk surface area, and F disk is the assumed blackbody seed photon flux with effective temperature T eff , set by integratingthe (target-temperature) cooling function within the disk at the specified ( r (cid:48) , φ (cid:48) ): (cid:90) Θ bot Θ top L dθ √ g θθ = 2 σ SB T . (12)We express cos ϑ in terms of the spherical coordinates of the coronal and photosphere cells:cos ϑ = ± rR [sin θ cos θ (cid:48) (cos φ cos φ (cid:48) + sin φ sin φ (cid:48) ) − cos θ sin θ (cid:48) ] = rR G ± ( r , r (cid:48) ) , (13)where + is used for the lower half of the corona and − for the upper half. In addition, it can beshown that the infinitesimal solid angle d Ω (cid:48) subtended by the disk surface element at r (cid:48) from thecoronal cell at r is: d Ω (cid:48) = 2 π (cid:34) − (cid:18) dA (cid:48) πR (cid:19) − / (cid:35) . (14)Substituting equations 13 and 14 into equation 11, and integrating over the disk surface, we arriveat nverse Compton Cooling u rad ( r ) = 2 π (cid:90) ∂A (cid:48) rR F disk ( r (cid:48) ) G ± ( r , r (cid:48) ) (cid:34) − (cid:18) dA (cid:48) πR (cid:19) − / (cid:35) . (15)This expression is computed in each coronal cell. To ease the computational burden of the numericalintegration of the right-hand side of the above equation, a coarsened sampling of the photosphere grid(e.g., including only every eighth φ grid-point and every sixth r grid-point) is used without significantloss in accuracy. Furthermore, u rad does not need to be computed every time step: harm3d fluidupdate time steps are guaranteed to be sufficiently short compared to the thermal time scale thatthe integral in equation 15 varies very little from one time step to the next. In practice, we haveverified that evaluating u rad only once every 20 time steps introduces <
1% error.2.2.
The Compton Temperature
It is also useful to estimate the Compton temperature, T C , in each coronal cell. The Comptontemperature—the temperature at which Compton heating is balanced by Compton cooling—is, inthe non-relativistic limit: k B T C = 14 (cid:82) ∞ hνJ ν dν (cid:82) ∞ J ν dν = 14 (cid:104) ε (cid:105) , (16)where J ν is the mean intensity at frequency ν . In other words, the Compton temperature is equalto one quarter the mean photon energy (so defined). For a pure blackbody, it is easy to show that (cid:104) ε (cid:105) = 3 . kT eff . Therefore, T C in a given coronal cell is found by averaging T eff over the disksurface, weighted by the contribution of each particular photosphere surface element’s flux to thetotal radiation energy density. That is: T C = 3 . (cid:80) n u rad ,n T eff ,n (cid:80) n u rad ,n , (17)where T eff ,n is the effective temperature of the n th disk surface element and u rad ,n is the evaluationof the right-hand side of equation 15 for a particular photosphere element (performed in the courseof numerical integration). The expression for IC power derived above, equation 8, is valid only if T e (cid:29) T C ; otherwise, an additional term for Compton heating is required:Compton heating = σ T m e c n e u rad (cid:104) ε (cid:105) . (18)As we show below, the T e (cid:29) T C condition is always met using the 1T assumption (and is still fairlywell satisfied under 2T), which indicates that Compton heating is negligible compared to Comptoncooling. 2.3. The IC Cooling Time
Equation 8 for the IC cooling power is the time-rate change in the internal energy of the gas. Thatis, in the fluid rest frame: dudt = −L IC . (19) Kinch et al.
Substitute the expression for Θ e in equation 10 (derived assuming T e = T i ) into equation 8 to solvethe differential equation above: u ( t ) = u (1 + bu ) e at − bu , (20)in which a ≡ σ T m e c χ χ ( c P /c V − u rad , (21) b ≡ m i m e ( c P /c V − χ ρc . (22)From equation 20, we calculate the cooling time t cool , or the time over which—assuming u rad and ρ are constant—the internal energy decreases from u → u /e : t cool = 1 a ln (cid:18) e + bu bu (cid:19) . (23)By inspection of the above equations, we see that the conditions for a short cooling time are eithera high radiation energy density or a high initial electron temperature. Because of the quadratic termin the expression for the IC cooling rate, using the instantaneous rate of equation 8 in harm3d can overestimate the cooling in very hot cells. We instead define the cooling over one time step byappropriately time-averaging the cooling rate during the time step:¯ L IC = u − u (∆ τ )∆ τ , (24)where ∆ τ is the proper time interval in the given coronal cell corresponding to the global simulation(coordinate) time step ∆ t . In addition, if t cool < ∆ τ in any cell, the global time step is resetaccordingly to match the shortest t cool . As mentioned above, the value of ¯ L IC is set assuming u rad and ρ are constant. While u rad is a function of (an average over) the disk structure—and will therefore varymore slowly—a given coronal cell’s density can of course change rapidly. To maintain the integrityof the numerical fluid dynamics solution, we must be sure that no cells cool too substantially eachtime step.In practice, ¯ L IC differs from L IC as given in equation 8, and t cool < ∆ τ , only briefly right after thenew corona cooling function is first “switched on.” Because the target-temperature cooling functiononly cools bound gas (and does so less efficiently, as we see below), the corona cools rapidly underthe new regime; the usual harm3d time step determination procedure is generally sufficient afterseveral-to-ten M of simulation time have elapsed.2.4. A Note on Units and Scaling
The derivations in the previous sections used physical, cgs units. To implement these equationsin harm3d , however, we must translate to code units. Using notation such that a quantity x isconverted from code to cgs units by x cgs = [ x ] x code , we rewrite equation 8 as:[ L ] L IC , code = 4 σ T cχm i [ ρ ] ρ code [ u rad ] u rad , code Θ e (1 + 4Θ e ) . (25) nverse Compton Cooling G = c = 1, Θ e (already a dimensionless quantity) is trivially re-expressedin code units by setting c = 1, u → u code and ρ → ρ code in equation 10. Consulting the conversionfactors for L and ρ from Schnittman et al. (2013), we have:[ L ] = 4 πc κG M ˙ m/η ˙ M code , (26)[ ρ ] = 4 πc κGM ˙ m/η ˙ M code , (27)[ u rad ] = [ u ] = c [ ρ ] . (28)Substituting these expressions into equation 25, we find: L code = (cid:18) πσ T χm i κ (cid:19) ˙ m/η ˙ M code ρ code u rad , code Θ e (1 + 4Θ e ) . (29)The term in parentheses is dimensionless and, assuming a fixed free electron fraction (totally-ionizedplasma is essentially guaranteed in the corona), constant.The other dimensionless term, ( ˙ m/η ) / ˙ M code , does not appear in the “code units” form of theexpression for the target-temperature cooling rate, but it does appear in the code units form ofthe IC cooling rate expression; it appears also in the expressions 26, 27 for translating the coolingrate and density from code to cgs units. This term serves to set the overall mass (and thereforedissipation) scale of the accretion flow. The radiative efficiency η is defined such that, in cgs units,the bolometric luminosity is related to the mass accretion rate by L = η ˙ M c . The value of η usedin the above expressions, however, must be chosen a priori . For ease of comparison with analyticaccretion disk theory, we choose the Novikov & Thorne (1973) values for the nominal radiativeefficiency; for a = 0, η NT = 0 . η NT , accreting at ˙ m = 0 .
01, has a luminosity equal to 0.01 L Edd . We describethese choices as “nominal” because their purpose is not to specify a resulting luminosity, but simplyto set a scale—because the radiative efficiency may not correspond exactly to the ISCO bindingenergy, but is computed by the simulation, the actual luminosity should be on order of, but is notexpected to be exactly equal to, 0.01 L Edd .Whereas the scale-free nature of the target-temperature cooling function required choosing in ad-vance only the dimensionless spin a/M , this is no longer the case when using the IC corona coolingfunction. The nominal accretion rate ˙ m appears in the expression for L code , as described above; also,the very first step when using this new cooling function is to divide the disk and corona by calculatingthe location of the photospheres—because ρ ∝ ˙ m , the upper and lower photosphere surfaces movefurther from the midplane with increasing ˙ m , decreasing the volume of the simulation space governedby the IC corona cooling function. The structure of the accretion geometry therefore depends on thechoice of ˙ m .For a chosen ˙ m , however, the simulation results are still scalable with M . To see why this remainspossible, consider how each term in equation 8 scales with M : ρ ∝ M − and u rad ∝ F disk ∝ L r ∝ M − M ; therefore L IC ∝ M − . Thus the IC corona cooling function has the same scaling with M Kinch et al. as does the target-temperature cooling function. The Compton temperature scales with M like so: T C ∝ T eff , disk ∝ ( L r ) / ∝ M − / . The scaling is weak, and such that T C decreases for more massiveblack holes. The condition for the validity of the IC cooling rate, equation 8, is T e = ( m e c /k B )Θ e (cid:29) T C . As discussed earlier, Θ e is a dimensionless quantity that does not scale with M or ˙ m (it isproportional to the ratio of two quantities with identical scaling relationships, ρ and u ). Therefore,because T e (cid:29) T C is satisfied for stellar-mass black holes (as we show below), it is necessarily satisifedfor supermassive black holes. 2.5. Uncooled Material harm3d ’s inversion routines are susceptible to a close subtraction error when recovering the gaspressure from the total pressure in regions which are magnetically-dominated [because p B ≈ p total , the“positive pressure problem” (Balsara & Spicer 1999)]. Numerical errors can result in large, artificialpressure gradients across adjacent cells, rapidly accelerating material in gross violation of energyconservation. As a remedy, harm3d instead solves an entropy conservation equation where solutionof the stress-energy conservation equation (2) fails; that is: ∇ µ ( S u µ ) = 0 , (30)where S ≡ p/ρ c P /c V − . Cells subject to evolution using the entropy equation do not obey the energyconservation equation, and are therefore unaffected by the specified cooling function; the effect issmall (Noble et al. 2009), however, and the region of the simulation volume to which it is appliedultimately contributes little to the overall simulation dynamics or X-ray observables. The entropyconservation equation is employed if either B /ρ (the magnetization) or B /u (approximately thereciprocal of the plasma β ) exceed certain critical values. For the simulations we show below, thethresholds chosen are B /u > , B /ρ > APPLICATION OF THE INVERSE COMPTON COOLING FUNCTIONTo demonstrate the IC cooling function with strongly coupled ions and electrons ( T e = T i ), weapply it to the zero spin ThinHR simulation, after the system has already evolved for 10 , M with the original, target-temperature cooling function applied everywhere. For our example case,the (scale free) ThinHR simulation is scaled to a 10 M (cid:12) central black hole with a nominal accretionrate of 1% Eddington ( ˙ m = 0 . M in two versions: with theIC cooling function in use in the corona, and with the original, target-temperature cooling functionremaining in use everywhere. We refer to the time at which at the IC cooling function switches onas t = 0 (even when referencing the target-temperature everywhere simulation).One of the chief assumptions of analytic accretion disk theory (Shakura & Sunyaev 1973) is that theflow is time-steady; this implies mass inflow equilibrium, i.e., the rate of inward mass flow throughshells at all radii is the same and does not vary with time. In a real or simulated system the accretionflow is turbulent and therefore highly variable in time and space; however, inflow equilibrium canstill be defined in a time-averaged sense. The shell-integrated mass inflow rate as a function of theradial coordinate is (Noble et al. 2012)˙ M ( r ) = − (cid:90) π dθ (cid:90) π/ dφ ρu r √− g ; (31) nverse Compton Cooling r/M − . − . . . . . . ˙ M ( r ) / ˙ M E dd − M M M Figure 1.
The shell-integrated mass inflow rate, ˙ M , as a function of radial coordinate, averaged in timeover three 1000 M windows, expressed in ratio to the nominal Eddington mass accretion rate. The dashedcurves represent the continuation of the target-temperature everywhere simulation, while the solid red andblue curves are for the run where the IC cooling function is switched on at t = 0. the factor of 4 and the azimuthal integration bounds are necessary as these simulations are per-formed over only one quadrant. In Figure 1, we show the radial-dependence of the mass inflowrate—expressed in ratio to the nominal Eddington mass accretion rate—time-averaged over threewindows: the last 1000 M of the ThinHR simulation, and the first and second 1000 M intervals after.For comparison, the same data is shown for the continuation past t = 0 for the original, target-temperature everywhere version.From Figure 1 it is apparent that these simulations achieve time-averaged inflow equilibrium outto r ∼ M ; though the value of the mass accretion rate (most sensibly measured at the eventhorizon at r = 2 M ) does vary. The ThinHR simulations are initialized with a finite amount ofmatter available to accrete (they are not “fed” by a companion as a real X-ray binary would be),some of which must be pushed out to larger radii as other material sheds its angular momentumand moves inward. In practice, it is challenging to design a simulation that achieves total radius-and time-independent inflow equilibrium while remaining computationally feasible. Nevertheless, theregion interior to r ∼ M accounts for about half of the total cooling in the simulation volume, andis the origin of the most important observable X-ray diagnostics.The volume-integrated total cooling rate is L tot = 4 (cid:90) R max R h dr (cid:90) π dθ (cid:90) π/ dφ L u t √− g, (32)2 Kinch et al. − −
500 0 500 1000 1500 2000 t/M − − L / L E dd totalray-traced coronadisk Figure 2.
The total cooling rate, the contributions from the disk and corona separately, and the luminos-ity measured by the distant observer via the ray-tracing code pandurata , as fractions of the Eddingtonluminosity, as functions of time. The IC cooling function is switched on (in the corona only) at t = 0. where L is the fluid rest frame value of the IC cooling function (in the corona) or the target-temperature cooling function (in the disk). R h = 2 M is the radius of the event horizon, and R max = 70 M is the outer radial boundary of the simulation volume. Figure 2 shows the totalcooling rate, as well as the contribution from the disk and corona separately, expressed in ratio to theEddington luminosity, as functions of time—including the last 1000 M of the “starter” simulation,with the division between disk and corona superimposed. With the old cooling function, the diskaccounted for nearly exactly half of the total cooling; with the new cooling function, the mean diskfraction is 0.38 and nearly constant in time—even though the overall luminosity decreases over thelength of the simulation.Note that L is defined in the fluid rest frame, therefore L as defined does not account for all specialand General Relativistic effects, nor for the capture of photons by the black hole. We thereforeuse the Monte Carlo ray-tracing code pandurata (Schnittman & Krolik 2013; Schnittman et al.2013), applied to successive snapshots of the harm3d simulation, to calculate the (post-processed)bolometric luminosity that reaches infinity. As is apparent from the figure, equation 32 provides aconsistent overestimation of the more careful ray-tracing calculation.As discussed in section 2.4, the conversion from code units to physical cgs units requires specificationof the radiative efficiency η . We choose the analytic accretion disk theory (Novikov & Thorne 1973)values for convenience. However, we can also compute the simulation’s radiative efficiency post hoc bycomparing the time-averaged accretion rate (as measured at the event horizon) to the time-averagedluminosity as measured by an observer at infinity. For the 0–1000 M window, this yields η = 0 . nverse Compton Cooling Figure 3.
Azimuthally-averaged values of the density (left) and the fluid frame cooling rate (right), for thesnapshot at t = 1000 M . The white lines indicate the ( φ -averaged) photosphere surfaces. Cells with zerocooling are shown in white. for the 1000–2000 M window, η = 0 . a = 0 black hole, η NT = 0 . M of the input simulation, η = 0 . L and ˙ M scale from code to cgs units proportional to the factor ( ˙ m/η NT ) / ˙ M code , and so theinferred radiative efficiency is nearly independent of the nominal choice. It is not entirely independent,however, as the location of the photosphere surfaces which divide the corona from the disk depends onthe physical density scale (equations 3 and 4). As Figure 2 indicates, the increase in inferred radiativeefficiency is due to the increased magnitude of the IC cooling function compared to the target-temperature version. Measured values of the radiative efficiency are generally not observationally-accessible, though long time-averages for L/L
Edd are. Regardless of our choice for the overall scalingfactor ( ˙ m/η ) / ˙ M code , our simulations can be just as well characterized (and thereby compared toobservable systems) by L/L
Edd . When the IC cooling function is turned on, the coronal luminositytriples, doubling
L/L
Edd for this system. Over time, the total luminosity trends downward, returningto its nominal value but with a greater share due to the corona. While not shown in Figure 2, the longterm trend of the “starter” target-temperature everywhere simulation continued for t > t = 1000 M for the IC simulation. Figure 3 shows the density and cooling rate for thissnapshot: note the rapid decrease in density away from the midplane (at 90 ◦ ). Azimuthal averagesof the dimensionless ratios B /ρ and B /u are shown in Figure 4. The critical values which triggerevolution via the entropy equation (see section 2.5) are met only in the relatively small region nearthe z -axis at r (cid:46) M .The IC simulation’s corona is nearly everywhere more luminous than is the target-temperatureversion’s. Compare the radial distributions of coronal cooling as shown in Figure 5 (shell-integrated4 Kinch et al.
Figure 4.
Azimuthally-averaged values of the ratio of magnetic energy density to mass-energy density (left)and magnetic energy density to internal energy density (right), for the snapshot at t = 1000 M . The whitelines indicate the ( φ -averaged) photosphere surfaces. and time-averaged). Except for a deficit at very small radii, the IC-cooled gas radiates more energy.The consequence is that the internal energy of the gas decreases after the IC cooling function isswitched on, resulting in an overall lower temperature corona (see Figure 16 below). Even thoughthe corona is generally cooler, it radiates more energy—this is not contradictory, because the coolingrate as determined with the target-temperature method is tied to the orbital time scale, not therelevant thermodynamic time scale. The rate of advection of thermal energy through the black holeevent horizon reduces by one-third after the IC cooling function is switched on, consistent with agreater fraction of injected heat being radiated away. As shown in Figure 1, the accretion rate increases after the IC cooling function is switched on, due to the rapidly cooling gas losing some ofits pressure support against gravity; eventually, however, the IC simulation’s accretion rate returnsto the same value of the continued target-temperature everywhere run (though its heat advectionrate remains lower).Even after the corona has cooled off, the IC cooling function simply radiates more energy per unittime than does the target-temperature cooling function, even when supplied with the same inputs.The magnitude of corona cooling shown in Figure 3 is 15–50 times greater than the target-temperaturecooling function would be for the same simulation data (the ratio is larger closer to the disk); thecooler corona at t = 1000 M is nearer to the target temperature, enhancing this difference.Figure 6 shows the fractional difference between the cooling rate per unit mass of the IC simulationdata and the continued target-temperature everywhere run, both evaluated at t = 1000 M . While thedisk body—which remains subject to the target-temperature cooling function in both instances—isunsurprisingly nearly the same, the IC-cooled corona is, overall, more luminous.The corona is even more magnetically supported ( B /u (cid:29) t = 1000 M than it was at the start, as the plasma β decreases by a factor of ten near the poles and by a factor ofone hundred near the disk. The magnetic field strength varies but generally increases, especially in nverse Compton Cooling r/M . . . . . . . . ∂ ( L c o r / L E dd ) / ∂ l og r ICtarget-temperature
Figure 5.
The radial distribution of the coronal cooling power, time-averaged, for the 0–2000 M intervalfor both the IC-cooled corona simulation and the continued target-temperature everywhere run. Figure 6.
The fractional difference between the ratio of local cooling rate to density, at t = 1000 M , betweenthe IC simulation and the target-temperature everywhere simualtion; azimuthally-averaged. Kinch et al. r/M . . . . . . . . H ph o t / r − M M M Figure 7.
The height of the photosphere, as a function of radius, averaged over azimuth, top and bottom,and for three time windows. While the location of the inner cutoff of the photosphere (the “reflection edge”)varies somewhat with time, it turns out to be close to the ISCO at r = 6 M . and near the disk, and the gas pressure falls everywhere. The increase in the magnetic field strengthis a trend observed both with and without switching on the IC cooling function; the large decrease inplasma β , however, is due to the falling gas pressure as the corona radiates away its internal energymore rapidly with the IC cooling function compared to the target-temperature version. The overallgeometry of the accretion flow is not significantly affected, however—as indicated in Figure 7, thedemarcation between corona and disk remains relatively fixed.The dynamical evolution of the corona is clearly affected by the choice of cooling function. Bythe end of their 2000 M runs, the IC-cooled simulation is advecting thermal energy through theevent horizon at a rate of 0.001 L Edd ; by contrast, the continued target-temperature everywhererun is advecting at a rate 0.0015 L Edd . At the same time, the IC simulation’s volume-integratedluminosity is 0.003 L Edd higher than the target-temperature run’s. The decreased advective loss,therefore, accounts for only one-sixth of the increase in luminosity of the IC simulation relativeto the continued target-temperature simulation. The magnetic heating rate of the coronal plasmaincreases as a consequence of the more realistic coronal cooling function. nverse Compton Cooling VALIDATION AGAINST A RAY-TRACING SOLUTION pandurata launches thermal seed photon packets from the disk surface and follows their trajec-tories through the corona until they either scatter off electrons, re-impinge on the disk surface, arecaptured by the black hole, or escape to an observer at infinity. When a photon packet scatters, itsenergy is convolved (in the fluid rest frame) with a thermal Compton scattering kernel correspondingto the presumed T e at the point of scattering [this procedure is described in Kinch et al. (2019)]. Byfollowing large numbers of photon packets in this manner, the net difference between each pre- andpost-scatter photon packet energy is used to construct the effective fluid rest frame inverse Comptoncooling rate for a given spatial map of T e . Through an iterative procedure, T e is adjusted everywhereuntil it matches harm3d ’s cooling map. At the same time, pandurata computes the energy- andinclination-dependent flux as would be measured by an observer at infinity. Because pandurata ’sprocedure and harm3d ’s IC cooling function calculation share several key assumptions—“fast light”and a local, thermal disk seed photon spectrum (determined for both via equation 12)—we can iso-late effects due to the ways in which they differ: pandurata accounts for all special and GeneralRelativistic effects, and the occlusion of disk radiation by intervening material, while harm3d doesnot.First, we examine harm3d ’s value for the radiation energy density in the corona, equation 15,compared to pandurata ’s. These are shown in Figure 8. pandurata ’s values are in fact calculatedby solving the equation for L IC we derived above for use in harm3d (equation 8) for u rad , using thesupplied ρ values and pandurata ’s values for L IC consistent with its own solution for T e ; pandurata does not explicitly calculate u rad . Immediately apparent in Figure 8 (right) is the poor Monte Carlosampling in the low density region near the poles. Even though a sufficient number of photonpackets are launched (1620 from each photosphere surface element) to ensure that the distant observerspectrum is well resolved, scattering events are simply so unlikely in the jet region that it is difficultto evaluate L IC there with pandurata ; therefore, pandurata values for u rad and T e are unreliable inthis region as well. In the much better sampled regions of the corona, it is apparent that pandurata ’s u rad values agree fairly well with harm3d ’s.In Figure 9, we show the azimuthally-averaged harm3d values for T e compared to pandurata ’sequilibrium T e values for the snapshot at t = 1000 M . In Figure 10, the comparison is presentedas the ratio of harm3d ’s T e values to pandurata ’s, as a function of radius (left) and polar angle(right), averaged, weighted by density, over the full 2000 M run of the simulation. In regions withpoor statistical sampling, pandurata will not adjust the temperature from its initial guess— T C asfound by harm3d —which is in effect a lower bound. Figure 11 shows the azimuthally-averaged ratioof the electron temperature to the Compton temperature. As expected (and required), T e (cid:29) T C .Thus in Figure 10 (right), the T e ratios tend to be very high within 30 ◦ of the z -axis. For much of the(well sampled) coronal volume, however, the two values for T e are within a factor of 2 of each other. harm3d ’s lower values for u rad (and correspondingly higher T e ) in the corona immediately above andbelow the inner disk are consistent with ignoring relativity: the relativistic inner disk orbital speedswill preferentially beam seed photons at angles more nearly parallel to the disk surface; in addition,the curved photon trajectories amplify u rad near the disk surface in the close vicinity of the blackhole. At large radii where relativistic effects are less important, ignoring the occlusion of interveningcorona material enhances u rad (lowering T e ) as compared to pandurata ’s value.8 Kinch et al.
Figure 8.
Azimuthally-averaged radiation energy density: left, from harm3d ’s calculation according toequation 15; right, the value “backed out” from pandurata output via equation 8. The white polar regionsindicate ( r, θ ) for which there were no scattering events during pandurata ’s Monte Carlo post-processing.The white region near the midplane is the disk body.
Figure 9.
Azimuthally-averaged electron temperature: left, from harm3d ; right, pandurata ’s equilibriumvalue. The white region near the midplane is the disk body.
Figure 12 shows the fraction of total coronal cooling which occurs within a certain polar angle fromthe midplane for four contiguous annuli of the simulation volume, averaged over time. We see thatthe majority of the cooling occurs within 45 ◦ of the midplane. And as shown in Figure 5, there islittle cooling within r (cid:39) M . Comparing the spatial distribution of corona cooling to the comparison nverse Compton Cooling r/M h T ha r m e / T p an d u r a t a e i ◦ ◦ ◦ ◦ ◦ ◦ ◦ θ h T ha r m e / T p an d u r a t a e i < r/M ≤ < r/M ≤ < r/M ≤ < r/M ≤ Figure 10.
The ratio of harm3d ’s T e value to pandurata ’s equilibrium T e value, averaged with densityweighting over the full 2000 M run. Left: also averaged over polar angle and azimuth, showing dependence ofratio on radius. Right: averaged over radius and azimuth for four contiguous annuli, showing dependence onpolar angle. The dashed red lines indicate the boundaries between which 1 / < (cid:104) T harm e /T pandurata e (cid:105) < Figure 11.
Azimuthally-averaged ratio of the harm3d -calculated values for the electron temperature tothe Compton temperature. Where the bulk of the luminosity occurs, T e /T C ∼ of T e values in Figure 10, we find that in the region which accounts for the majority of the totalcooling, harm3d ’s estimation of T e agrees fairly well—at the 10–20% level—with pandurata ’s.From this we conclude that, on net, the procedure detailed in the previous section for calculating the0 Kinch et al. ◦ ◦ ◦ ◦ ◦ ◦ ◦ θ (from midplane)0 . . . . . . L I C ( ϑ < θ ) / L I C , t o t < r/M ≤ < r/M ≤ < r/M ≤ < r/M ≤ Figure 12.
The fraction of total corona cooling, for four contiguous annuli, as a function of polar angle asmeasured from the midplane; time-averaged.
IC cooling rate is generally consistent with a more detailed ray-tracing radiation transport (thoughstill fast light) approach, at a fraction of the cost. nverse Compton Cooling − energy [keV]10 ε L ε (cid:2) e r g s − c m − (cid:3) − M M M Figure 13.
The spectral luminosity as seen by a distant observer, averaged over three time intervals.5.
THE EFFECT OF THE INVERSE COMPTON COOLING FUNCTION ON THEOBSERVED SPECTRUMFollowing the procedure laid out in Schnittman et al. (2013) and Kinch et al. (2016, 2019), wegenerate energy- and inclination-dependent simulated spectra by counting those photon packets in pandurata that escape to infinity. The spectral luminosity (integrated over inclination angle) aswould be seen by distant observers, averaged over time, is shown in Figure 13 for the last 1000 M of thestarter target-temperature simulation and for the first and second 1000 M intervals of the simulationwith IC cooling function in place. Note that pandurata is applied to multiple snapshots (20 evenlyspaced in time) in each interval separately—averaging is performed on the output spectrum, not theunderlying simulation data. The features of the predicted spectrum are qualitatively similar betweenthe three intervals: a broad thermal peak centered at (cid:39) (cid:38)
60 keV. However, with the IC cooling function in place, a substantially harder power-lawcomponent is achieved. Figure 14 shows the variation with time of the photon index Γ ( L ε /ε ∝ ε − Γ ),measured in the range 2–30 keV. Note that while the bolometric luminosity declines in the finalinterval, Γ is relatively unchanged.The spectral differences between the target-temperature and IC cooling function simulations areexplained by the distribution in temperature—as calculated by pandurata —of the cooling coronagas. Figure 15 shows, for the last snapshot of the target-temperature simulation and a representativesnapshot of the IC cooling function data, the spread in the pandurata -determined temperature ofthe cooling gas. In both cases, the coronal gas radiates over a broad range in temperature—however,the relative distribution with temperature is very different. At least for this choice of parameters,the ad hoc target-temperature cooling function posited a cooling rate that was unphysically low for2 Kinch et al. − −
500 0 500 1000 1500 2000 t/M . . . . . . Γ Figure 14.
The photon index Γ, fit to the range 2–30 keV of the spectral luminosity, as a function of time. the hot gas, and unphysically high for the cool gas. The result is a softer power-law component inthe X-ray spectrum.Figure 16 shows the coronal mass distribution by temperature for the − M and 1000 M snap-shots, with T e as calculated by harm3d (dashed lines) and T e as calculated by pandurata (solidlines). Note that the agreement between the harm3d and pandurata distributions is substan-tially improved for the 1000 M snapshot, for which harm3d uses the IC cooling function. Whilethe harm3d -determined mass-by- T e distribution has a broader spread in temperature than its pan-durata -determined counterpart for the t = 1000 M snapshot, their peaks—i.e., where most of thematerial lies—agree. By its nature as a Monte Carlo code, pandurata will best sample the mostdense region just outside the disk (recall Figure 8, right); the high altitude regions where harm3d -and pandurata -determined temperatures most differ (Figure 10, right), and where only a smallfraction of the total cooling takes place (Figure 12), account for the wider breadth the harm3d distribution. nverse Compton Cooling T e [eV]0 . . . . . . . . ∂ ( L c o r / L E dd ) / ∂ l og T e M − M Figure 15.
The distribution of coronal cooling with respect to its pandurata -determined electron tem-perature. Two snapshots are shown: just before the IC cooling function is turned on (red) and 1000 M later(black). T e [eV]0 . . . . . . . . ∂ ρ / ∂ l og T e M pandurata M harm3d − M pandurata − M harm3d Figure 16.
The distribution of coronal mass with respect to electron temperature; T e is determined by pandurata for the dashed lines and by harm3d for solid lines. Two snapshots are shown: just before theIC cooling function is turned on (red) and 1000 M later (black). Kinch et al. TWO TEMPERATURE INVERSE COMPTON COOLING FUNCTIONBelow we describe our procedure for relaxing the assumption that T e = T i in a particular scenario:1. The ion and electron populations are individually in thermal equilibrium locally, described by T i and T e at each point in the corona.2. Turbulent energy is dissipated into the ions only .3. Energy is exchanged between the ion and electron populations through Coulomb collisions.4. T e adjusts instantaneously so that the rate at which energy is added to the electron population—through either Coulomb collisions or Compton heating—is equal to the rate at which energy is lost due to inverse Compton cooling.The problem is to find T e and T i given the above assumptions.We begin with the relativistically correct ion-electron energy exchange rate derived in Stepney(1983); Stepney & Guilbert (1983): ddt u e (Coulomb heating) = 32 m e m i σ T c ln Λ n e n i ( k B T i − k B T e ) (33) × (cid:26) K (1 / Θ e ) K (1 / Θ i ) (cid:20) e + Θ i ) + 1Θ e + Θ i K (cid:18) Θ e + Θ i Θ e Θ i (cid:19) + 2 K (cid:18) Θ e + Θ i Θ e Θ i (cid:19)(cid:21)(cid:27) = 32 m e m i σ T c ln Λ n e n i ( kT i − kT e ) f (Θ e , Θ i ) , where u e is the internal energy per unit volume of the electron population, ln Λ is the Coulomblogarithm (the logarithm of the ratio of the maximum to minimum impact parameters; ln Λ ∼ i is the dimensionless ion temperature, equal to k B T i /m i c , and K n is the n th order modified Besselfunction of the second kind. The term in braces constitutes f (Θ e , Θ i ).Electrons gain or lose energy to ions through Coulomb collisions (equation 33), gain energy fromphotons through Compton heating (equation 18) and lose energy to photons through Compton cooling(equation 8). The equilibrium T e is that for which these processes balance, du e /dt = 0. Setting thesum of these three processes to zero, we have:32 m e m i σ T c ln Λ n e n i ( k B T i − k B T e ) f (Θ e , Θ i ) + σ T m e c n e u rad (cid:104) ε (cid:105) − σ T cn e u rad Θ e (1 + 4Θ e ) = 0 . (34)With the identification that n e = χn i and n i = ρ/m i , and some algebraic manipulation, we rewritethe above as: 32 m e m i ln Λ (cid:18) Θ i − m e m i Θ e (cid:19) f (Θ e , Θ i ) + u rad ρc (cid:104) ε (cid:105) m e c − u rad ρc Θ e (1 + 4Θ e ) = 0 . (35)We rearrange the ideal gas law, equation 9, to solve for Θ i ,Θ i = ( c P /c V − uρc − χ m e m i Θ e . (36) nverse Compton Cooling m e m i ln Λ (cid:20) ( c P /c V − uρc − m e m i (1 + χ )Θ e (cid:21) f (cid:16) Θ e , ( c P /c V − uρc − χ m e m i Θ e (cid:17) (37)+ u rad ρc (cid:104) ε (cid:105) m e c − u rad ρc Θ e (1 + 4Θ e ) = 0 . We define three dimensionless quantities, A ≡ uρc , B ≡ u rad ρc , and C ≡ (cid:104) ε (cid:105) m e c , (38)which we substitute into equation 37:32 m e m i ln Λ (cid:20) ( c P /c V − A − m e m i (1 + χ )Θ e (cid:21) f (cid:16) Θ e , ( c P /c V − A − χ m e m i Θ e (cid:17) (39)+ BC − B Θ e (1 + 4Θ e ) = 0 . We arrive at an equation for Θ e , in terms of only dimensionless quantities, which depends on threeparameters which can be read off directly from already-computed values in harm3d . Equation 39 isnot amenable to real-time solution in each coronal cell each timestep; rather, we tabulate its solutionon a grid covering all possible, reasonable values of A , B , and C , and use trilinear interpolationto calculate the appropriate equilibrium value of Θ e from said lookup-table in the course of thesimulation run. From Θ e we calculate Θ i from equation 36. The cooling function which enters into harm3d is the net Compton cooling translated to code units.Figure 17 is a demonstration of our “snap to equilibrium” approximation for T e . It shows asimple forward Euler integration of the Compton and Coulomb heating/cooling equations, with n e = 3 × cm − and u rad = 3 × erg cm − , values typical for the corona of a 10 M (cid:12) blackhole accreting at 1% Eddington. For t < T i = 60 keV. At t = 0, the ion temperature jumps to120 keV. The adjustment of T e to its new equilibrium value takes (cid:39) . M . We approximate suchan adjustment to be instantaneous. Of course, n e and u rad —which dictate this re-equilibration timescale—vary broadly in time and space in the corona. The range of actually encountered values for n e and u rad (in regions where there is any substantial cooling) correspond to adjustment times in therange 0 . M . COMPARISON OF THE 1T AND 2T COOLING FUNCTIONSThe two-temperature (2T) IC cooling function is applied to the same starter simulation to which weapplied the 1T IC cooling function in the previous section. With the 2T IC cooling function switchedon, the system is evolved for 1000 M . In Figure 18 we compare the total and corona-only volume-integrated cooling rates for both approaches. Note that because the IC cooling function is appliedonly in the corona, the disk component for each approach is not shown as they are nearly identical.Not surprisingly, the 2T method results in a lower corona luminosity; the mean disk fraction is, aswith the old target-temperature cooling function, nearly exactly half. The mass accretion rate atthe event horizon is also somewhat lower (Figure 20). Comparing the time-averaged luminosity (the“ray-traced” power—not shown in Figure 18 for clarity, though as in Figure 2, it is nearly the same as6 Kinch et al. − . − . − .
25 0 .
00 0 .
25 0 .
50 0 . t/M . . . . . . T e [ k e v ] Figure 17.
An integration for T e of the Coulomb and Compton rate equations using values typical for thecorona, in a scenario where the ion temperature doubles at t = 0. The electron temperature adjusts to itsnew equilibrium value after approximately 0 . M . Our model treats this as instantaneous (the red dashedline). the total volume-integrated cooling rate), we find that the 2T approach yields a radiative efficiencyof η = 0 . M window). The standard deviationof the moving time-averaged coronal luminosity, L IC / (cid:104) L IC (cid:105) , is twice as large for the 1T simulationthan for the 2T simulation. This is consistent with the demonstration in Figure 17: in that example, T i jumped from 60 keV to 120 keV, while the equilibrium T e value increased by less than 1 keV.Because IC power is a function of the electron temperature—not the ion temperature—the energyexchange between ions and electrons dampens variability.Indeed, Coulomb collisions act as a bottleneck between MHD heating and IC cooling. Figure 21shows the azimuthally-averaged cooling rate for the 2T simulation at t = 1000 M . Compared to thesame plot for the 1T simulation (Figure 3, right), the cooling is confined to a smaller region, closerto the disk. The ion-electron exchange rate is proportional to the square of the density—furtherfrom the midplane, the density is simply too low to support Coulomb heating of the electrons. Thetemperature plots of Figures 22 and 23 show this clearly. Compare to the plot of T e (and, byassumption, T i ) for the 1T simulation in Figure 9 (left): the ions are hotter and the electrons arecooler. The ratio T e /T i declines sharply from near unity just outside the disk photosphere to < − in the jet cone. The decreased ion cooling rate with the 2T method relative to the 1T accounts forthe decreased mass accretion rate: greater ion temperature translates to a stronger pressure support nverse Compton Cooling t/M . . . . . . . L / L E dd total (2T)corona (2T)total (1T)corona (1T) Figure 18.
The total cooling rate, with the contribution from the corona only, for both 2T and 1Tsimulations, as fractions of the Eddington luminosity, as functions of time. t/M . . . . . L I C / h L I C i Figure 19.
The volume-integrated coronal cooling rate, divided by a 25 M window moving average of thesame data, for both 2T and 1T simulations. Kinch et al. r/M − . − . . . . . . ˙ M ( r ) / ˙ M E dd Figure 20.
For both 2T and 1T simulations: the shell-integrated mass inflow rate, ˙ M , as a function ofradial coordinate, averaged in time over the 0–1000 M window, expressed in ratio to the nominal Eddingtonmass accretion rate. against gravity. Note also that while T e > T C in the 2T case as well as the 1T, the ratio here issmaller, T e /T C (cid:38) pandurata to successive snapshots of the 2T simula-tion as well. Figure 24 shows the spectrum as seen by an observer at infinity, for both the 1T and 2Tsimulations, each averaged over the range 0–1000 M . The 2T simulation produces a notably softerX-ray power-law: Γ = 2 .
53 compared to the 1T simulation’s Γ = 2 .
25 (measured on 2–30 keV). Inaddition, the rollover occurs at a lower energy and falls off more sharply. This is consistent with anoverall less-luminous, cooler corona. Compare the distribution of cooling with temperature in Figure25: not only is the 2T curve shifted toward cooler temperatures, but the tail above 100 keV in the1T data—physically located more than 45 ◦ from the midplane (see Figure 9)—is totally absent inthe 2T data. nverse Compton Cooling Figure 21.
Azimuthally-averaged values of the fluid frame cooling rate, for the snapshot at t = 1000 M ofthe 2T simulation. The white lines indicate the ( φ -averaged) photosphere surfaces. Cells with zero coolingare shown in white. Figure 22.
Azimuthally-averaged values of the ion temperature (left) and electron temperature (right).Note the dramatic difference in scales between the two plots. The disk body (where no distinction is made)is shown in white. Kinch et al.
Figure 23.
Azimuthally-averaged values of the ratio of the electron to ion temperatures (left) and electronto Compton temperatures (right). The disk body is shown in white. − energy [keV]10 ε L ε (cid:2) e r g s − c m − (cid:3) Figure 24.
The spectral luminosity as seen by a distant observer, averaged over 0–1000 M for both 1T and2T data. nverse Compton Cooling T e [eV]0 . . . . . . . . ∂ ( L I C / L E dd ) / ∂ l og T e Figure 25.
The distribution of coronal cooling with respect to the gas electron temperature, for the 2Tand 1T simulations each at t = 1000 M . 8. CONCLUSIONWe have shown that our simplified calculation of the radiation energy density in the corona requiredto compute the inverse Compton cooling function is reasonably accurate compared to the much morecareful pandurata ray-tracing calculation—especially so in the regions that account for the majorityof the coronal cooling. Our confidence in the method is further bolstered by the fact that the time-averaged post-processed continuum spectra are qualitatively similar to real X-ray observations. Thepower-law index for the 1T simulation run (Γ = 2 .
25) and its disk fraction (0.38) place it just barelytoo soft to meet the standard parameters for a classical “hard” spectrum (Remillard & McClintock2006), but it is qualitatively similar. The 2T run results in a power-law that is substantially softerthan most X-ray binary spectral observations; to the extent that these simulations adequately capturethe effects real coronal physics has on the observed X-rays, we conclude that it is likely that somecoupling mechanism stronger than Coulomb collisions is at work in real systems.We must note, however, that real X-ray spectra of stellar-mass black holes are typically integratedover tens of thousands of seconds; by contrast, the full 2000 M run with the new cooling function inplace is only 0.01 seconds long. This run cost about 80,000 core-hours to perform on a modern highperformance computing cluster. Thus, simulating for a substantial fraction of a real observation is notcomputationally feasible. It is for exactly this reason, however, that we were motivated to develop arealistic inverse Compton cooling function that does not require the additional computational burdenof real transport: harm3d using the IC cooling function in the corona runs at essentially the samespeed as with the target-temperature cooling function everywhere, but is considerably more physical.The same cannot be said for any treatment involving actual radiation transport.2 Kinch et al.
The success of this method applied to the a = 0 case motivates us to apply it to spinning black holesimulations as well. Past efforts to treat the coronal equation of state for black hole accretion havebeen frustrated by the very large computational expense of true radiation coupling [e.g., Jiang et al.(2019b,a)], while target-temperature cooling functions without real radiation transfer [e.g., Nobleet al. (2009, 2010) and Shafee et al. (2008); Penna et al. (2010)] are physically unrealistic. Our newcoronal cooling function permits efficient computation while simultaneously providing a very goodapproximation to the actual cooling rate of coronal plasma; we have also shown that a more realisticcoronal cooling function can alter coronal dynamics sufficiently to change the luminosity by tens ofpercent. Moreover, when these simulations provide the input data for full-up disk-atmosphere pluscoronal modeling [via pandurata + ptransx : Kinch et al. (2016, 2019)], it becomes possible forthe first time to make credible predictions of spectral features, including Fe K α emission, for a widerange of accretion rates onto black holes of all masses and spin parameters.ACKNOWLEDGMENTSBEK thanks the members of the Center for Theoretical Astrophysics at Los Alamos NationalLaboratory for helpful conversations, and John G. Baker (NASA GSFC) for useful geometric for-mulae. BEK was supported by the U.S. Department of Energy Advanced Simulation and Com-puting Program’s Metropolis Fellowship, through the Los Alamos National Laboratory, and usedresources provided by the Los Alamos National Laboratory Institutional Computing Program. LosAlamos National Laboratory is operated by Triad National Security, LLC, for the National NuclearSecurity Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001).SCN was supported by NSF awards AST-1515982 and OAC-1515969, NASA TCAN award TCAN-80NSSC18K1488, and by an appointment to the NASA Postdoctoral Program at the Goddard SpaceFlight Center administrated by USRA through a contract with NASA. JDS was supported by NASATCAN award TCAN-80NSSC18K1488. JHK and BEK were supported by NSF awards AST-1516299,CDI-1028111, and PHYS-1707826. REFERENCES Abramowicz, M. A., Lanza, A., & Percival, M. J.1997, ApJ, 479, 179, doi: 10.1086/303869Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376,214, doi: 10.1086/170270Balsara, D. S., & Spicer, D. 1999, Journal ofComputational Physics, 148, 133,doi: 10.1006/jcph.1998.6108Blumenthal, G. R., & Gould, R. J. 1970, Reviewsof Modern Physics, 42, 237,doi: 10.1103/RevModPhys.42.237Chandrasekhar, S. 1960, Radiative transferFragile, P. C., Gillespie, A., Monahan, T.,Rodriguez, M., & Anninos, P. 2012, ApJS, 201,9, doi: 10.1088/0067-0049/201/2/9Fragile, P. C., Olejar, A., & Anninos, P. 2014,ApJ, 796, 22, doi: 10.1088/0004-637X/796/1/22 Gammie, C. F., McKinney, J. C., & T´oth, G.2003, ApJ, 589, 444, doi: 10.1086/374594Haardt, F., & Maraschi, L. 1991, ApJL, 380, L51,doi: 10.1086/186171Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376,223, doi: 10.1086/170271Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis,S. W. 2019a, ApJ, 885, 144,doi: 10.3847/1538-4357/ab4a00Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014a,ApJS, 213, 7, doi: 10.1088/0067-0049/213/1/7—. 2014b, ApJ, 796, 106,doi: 10.1088/0004-637X/796/2/106—. 2019b, ApJ, 880, 67,doi: 10.3847/1538-4357/ab29ff nverse Compton Cooling33