The Origin of the Slow-to-Alfvén Wave Cascade Power Ratio and Its Implications for Particle Heating in Accretion Flows

The partition of turbulent heating between ions and electrons in radiatively inefficient accretion flows plays a crucial role in determining the observational appearance of accreting black holes. Modeling this partition is, however, a challenging problem because of the large scale-separation between the macroscopic scales at which energy is injected by turbulence and the microscopic ones at which it is dissipated into heat. Recent studies of particle heating from collisionless damping of turbulent energy have shown that the partition of energy between ions and electrons is dictated by the ratio of the energy injected into the slow and Alfvén wave cascades as well as the plasma β parameter. In this paper, we study the mechanism of the injection of turbulent energy into slow- and Alfvén-wave cascades in magnetized shear flows. We show that this ratio depends on the particular (r ϕ) components of the Maxwell and Reynolds stress tensors that cause the transport of angular momentum, the shearing rate, and the orientation of the mean magnetic field with respect to the shear. We then use numerical magnetohydrodynamic shearing-box simulations with background conditions relevant to black hole accretion disks to compute the magnitudes of the stress tensors for turbulence driven by the magneto-rotational instability and derive the injection power ratio between slow and Alfvén wave cascades. We use these results to formulate a local subgrid model for the ion-to-electron heating ratio that depends on the macroscopic characteristics of the accretion flow.


INTRODUCTION
Recent Very Long Baseline Interferometric (VLBI) observations by the Event Horizon Telescope of the horizon-scale environments of the black holes at the centers of the Milky Way (Sgr A*) and M87 galaxies have offered a unique opportunity to study plasma astrophysics in low-luminosity accretion flows (Event Horizon Telescope Collaboration et al. 2019;Akiyama et al. 2022).These systems are categorized as Radiatively Inefficient Accretion Flows (see Narayan & Yi 1995) owing to their low luminosities and the characteristics of their spectra.They are made up of low density plasmas, with mass accretion rates typically less than 10 −3 times the Eddington rate ṀEdd .The collisional timescales between the ions and the electrons in these flows are much larger than accretion timescale, allowing for the ions and the electrons to co-exist in a two-temperature state.In addition, low densities also cause the cooling processes in these systems to be inefficient, resulting in most of the energy carried by the plasma to be advected into the black hole.
Due to the two-temperature nature of the plasma, a firstprinciples approach to studying accretion flows involves independently evolving the thermodynamics of the ions and the electrons, incorporating physical models for the partition of heat between the species (Ressler et al. 2015;S ądowski et al. 2017).In these accretion flows, heating mechanisms include collisionless damping due to wave-particle resonances (Quataert 1998;Kawazura et al. 2020), magnetic reconnec-tion (Ball et al. 2018;Rowan et al. 2019), and heating resulting from instabilities due to velocity-space anisotropies (Sharma et al. 2007).These processes are typically dominant at the length-scales comparable to the gyroradii of the ions and the electrons, which are much smaller than the scale of the system.As a result of the large scale separation, the physics of heating mechanisms studied at the microscopic kinetic scales need to be incorporated as sub-grid prescriptions in global simulations.
Among these dissipation processes, magnetic reconnection is primarily dominant in localized regions of the flow that have current sheets (Ball et al. 2016) and causes episodic dissipation of energy.On the other hand, velocity space anisotropies are believed to be driven by mechanisms connected to the time-evolving large scale magnetic fields.In the absence of current sheets or local mechanisms driving velocity space anisotropies, the turbulent energy cascades down to scales comparable to the ion gyroradius and undergoes collisionless damping.This is a ubiquitous mechanism of dissipation, both spatially and temporally and is expected to be the primary driver of electron and ion heating.
Recent numerical studies show that the partition of heat between ions and electrons resulting from collisionless damping of the turbulent energy (Kawazura et al. 2020) depends on the plasma β (the thermal to magnetic pressure ratio) and the ratio of the driving power of the compressive slow magnetosonic wave cascade to that of the Alfvén wave cascade (P S /P A ).These studies have been performed in the gyroki-netic limit, where the gyromotion of particles around the mean magnetic field are averaged out.In such a limit, the compressive slow waves and the Alfvén waves in the turbulence decouple at the leading order at scales smaller than the driving scales of the turbulence (Schekochihin et al. 2009).Owing to this inertial range, the ratio P S /P A does not only impact the microphysical partition of heat between ions and electrons but is also evident in the ion and electron temperature profiles in an accretion flow (Satapathy et al. 2023), underscoring the importance of studying the nature of the slow and Alfvén driving caused by the MRI.
The ratio of the driving power of the slow-and the Alfvénwaves, in turn, is determined by the mechanism that drives turbulence in the flow.In black hole accretion flows, the magneto-rotational instability (MRI, see Balbus & Hawley 1991) is believed to drive the turbulence and cause the outward transport of angular momentum.The physics of the injection of energy into the slow and Alfvén wave turbulent cascades caused by the MRI has not been extensively explored.Only one study, performed by Kawazura et al. (2022), addresses the problem numerically using reduced-MHD shearing box simulations.These authors studied the MRI driven by a near-azimuthal mean magnetic field in a Keplerian shear flow and showed that the slow-wave cascade is driven approximately at twice the power as the Alfvén-wave cascade.
In this paper, we study the slow-to-Alfvén wave driving power ratio P S /P A and establish its relation to the largescale properties of the turbulence.We show that the ratio P S /P A depends on particular components of the Maxwell and Reynolds stress tensors that cause the outward transport of angular momentum in the disk, the orientation of the mean magnetic field, and the shearing rate.We then employ local shearing box simulations in presence of different mean magnetic field configurations and shearing rates relevant to accretion flows to compute the two stress tensors resultant due to the MRI and its parasitic instability in a non-linear turbulent state.Subsequently, we infer the slow-to-Alfvén driving power ratio P S /P A for these configurations.In combination with hybrid kinetic scale simulations performed by Kawazura et al. (2020), we formulate a local sub-grid model for the partition of turbulent heat between the ions and the electrons resultant from collisionless damping of slow-and Alfvén-waves.We also resolve a confusion in the literature regarding the compressive nature of slow waves in the reduced-MHD limit and establish that the Helmholtz decomposition is not useful in calculating the ratio P S /P A (see Appendix A).
In Section 2, we review the effects of the large scale turbulent driving on the dissipation scale partition of heat into the ions and electrons resultant from the collisionless damping of the turbulence.In Section 3, we set up equations of the shearing box approximation of MHD and derive the energy equations for the slow magnetosonic and the Alfvén wave fluctuations valid to the first order in the reduced-MHD limit.In Section 4, we discuss the physics of the injection of energy into the slow and the Alfvén wave cascades as caused by the MRI and establish the relation of the injection ratio to the Maxwell and Reynolds stress tensors.In Section 5, we compute the value of P S /P A for a wide range of conditions relevant to the turbulence driven in black hole accretion flows using a suite of shearing box simulations.We summarize our findings in Section 6.

EFFECTS OF TURBULENT DRIVING ON COLLISIONLESS DAMPING
The turbulence in plasmas found in radiatively inefficient accretion is believed to be strongly anisotropic in nature, caused by the presence of strong magnetic fields.It is useful to analyze such plasmas in a reduced-MHD limit, where fluctuation amplitudes, mode anisotropy, and fluctuation timescales are treated at linear order.In this approximation, the slow-wave and Alfvén-wave cascades are energetically decoupled in the inertial range of turbulence (Schekochihin et al. 2009).The slow-and the Alfvén-wave cascades only undergo mode mixing at length-scales comparable to the ion Larmor radius and undergo collisionless damping, channeling a fraction of the turbulent energy into heating the ions.The remaining energy cascades down to smaller length scales and dissipates at the electron Larmor radius, heating the electrons.In Figure 1, we represent the typical scales in which the above processes operate.
In a recent paper, Kawazura et al. (2020) studied the dependence of the heat partition between the ions and the electrons on the composition of the wave cascades driving the turbulence.They established that the heating ratio between the ions and the electrons, Q i /Q e , primarily depends on the ratio of the driving powers of compressible slow waves and the incompressible Alfvén waves P S /P A , and on the plasma β.This dependence can be approximated using the expression In the limit of plasma β ≲ 1, the value of Q i /Q e is almost entirely dependent on the composition of the driving-scale turbulent cascade, with Q i /Q e ∼ P S /P A .Additionally, in plasmas which are dominantly driven by slow waves, nearly all of the energy carried by the slow waves is channeled into heating the ions irrespective of the value of plasma β, causing The second term in equation ( 1) captures the dependence of the heating ratio on the microphysics of collisionless damping of turbulence that occurs at the small lengthscales comparable to the ion gyroradius.On the other hand, the first term involving the ratio of the slow-to-Alfvén driving power, P S /P A , is determined by the large scale processes that drive turbulence at large scales.Our goal is to explore this ratio in the general setting applicable to RIAFs.
In such flows, MRI has been identified as the most likely mechanism that drives the turbulence responsible for the transport of angular momentum (Balbus & Hawley 1998).The net outward transport of angular momentum is caused by particular off-diagonal components of the Maxwell and Reynolds stress tensors.The same components of these stresses are also responsible for the injection of turbulent energy into the accretion flow (Hawley et al. 1995).To that end, we compute the rates of the injection of turbulent energy into the individual slow-and Alfvén-wave cascades, and show that we can relate them to the Maxwell and Reynolds stress tensors in a local region of the accretion flow.
In order to infer properties of the turbulence in accretion flows, we will study the turbulence in the shearing-box approximation (Goldreich & Lynden-Bell 1965).This approximation allows us to capture the mechanism of the injection of turbulent energy, i.e., the MRI, in a local patch of an accretion flow.In this paper, we set up a local shearing box with mean background magnetic fields ( ⃗ B 0 ) and shearing rates (q) relevant to black hole accretion disks, as we describe in the following section.

THE SHEARING BOX ENERGY EQUATIONS IN A REDUCED MHD LIMIT
In this section, we derive the energy equations for the slow magnetosonic and the Alfvén wave fluctuations in a shearing box MHD system.We first set up the equations of MHD in a shearing box configuration.We then expand the equations to first order in the reduced-MHD approximation.In this limit, we discuss the orientations of the eigenvectors of the normal modes of MHD and derive the reduced-MHD energy equations for the Alfvén-and the slow-wave turbulent cascades.

The Shearing Box Approximation in the Reduced MHD Limit
The equations of MHD in the shearing box approximation can be written as (Goldreich & Lynden-Bell 1965) and In the above equations, ρ denotes the plasma density, P the pressure, ⃗ u the velocity perturbation on the background shear flow ⃗ u 0 , ⃗ B the magnetic field strength, and ⃗ Ω the angular velocity of the system.We set up the equations in a Cartesian comoving coordinate system (x, y, z) such that the angular velocity ⃗ Ω ≡ Ω ẑ and x represents the radial outward orientation in a global disk.The background shear flow ⃗ u 0 is given by where q ≡ d(ln Ω)/d(ln r) is the shearing rate.We assume an adiabatic equation of state, P = Kρ γ , where γ is the adiabatic index and K is a constant.In a magnetized system described by the equations of MHD, the linear normal modes consist of the Alfvén, the slow magnetosonic, and the fast magnetosonic waves.The direction of the mean magnetic field (denoted by ⃗ B 0 ) sets the orientation of the fluctuations in each of these normal modes.The orientations of the velocity eigenvectors in Fourier space for each of the normal modes can be written as (see Cho & Lazarian 2003) for the Alfvén wave eigenvector ξA , for the slow wave eigenvector ξS , and for the fast wave eigenvector ξF .In equations ( 6)-( 8), the direction of the mean magnetic field sets the orientation of the parallel wavevector, with k∥ = B0 .The perpendicular The velocity eigenvectors for the Alfvén ( ξA, green), the slow magnetosonic ( ξS, blue), and the fast magnetosonic ( ξF , red) MHD waves.The green dot for the Alfvén eigenvector represents an out of plane orientation.The slow and the fast mode eigenvectors depend on the value of α.The α → 0 and α → ∞ limits of these eigenvectors are indicated by the labels.The reduced-MHD limit of the Alfvén eigenvector ( ξR A ) is the same as the MHD eigenvector, while the same limit of the slow wave eigenvector ( ξR A ) is along the direction of the mean magnetic field (shown in the bottom-right of this figure).wavenumber ( k⊥ ) is defined such that ⃗ k = ⃗ k ∥ + ⃗ k ⊥ .The parameter α is defined as the square of the ratio of the sound speed (c s ) to the Alfvén speed (v A = B 0 / √ 4πρ 0 ), i.e., The quantity D is given by D Figure 2 shows the set of these three orthogonal eigenvectors in Fourier space.As denoted by the green dot, Alfvén waves are always azimuthally polarized with respect to the mean magnetic field (eq.[6]).The slow (shown in blue) and fast waves (shown in red) lie on the k ∥ − k ⊥ plane with the parameter α determining their orientations.In a magnetically dominated plasma (α → 0), the slow waves are polarized nearly parallel to the mean magnetic field and the fast waves are polarized nearly perpendicular to the mean magnetic field.On the other hand, in a pressure-dominated plasma where α → ∞, the fast waves are polarized along k.Importantly, the slow waves in such a limit are nearly incompressible and are polarized in the − θ ≡ ( k⊥ × k∥ ) × k direction.They are also referred to as the pseudo-Alfvén waves (Maron & Goldreich 2001).
The reduced MHD limit of the MHD equations is obtained by expanding the equations to first order in mode anisotropy k ∥ /k ⊥ and in the perturbations of density (δρ/ρ 0 ), velocity (u/v A ), and magnetic fields (δB/B 0 ) (see Schekochihin et al. 2009).The frequencies ω of the fluctuations in all the fluid quantities are considered to be of the order of the frequencies of the Alfvén waves (ω ∼ k ∥ v A ).As a result, the fast wave displacements exhibit timescales much shorter than the characteristic timescale 1/ω and, hence, are ordered out.
We now examine the polarizations of the Alfvén-and the slow-waves in this limit by expanding the eigenvectors ξA and ξS , respectively, to the leading order in k ∥ /k ⊥ .The polarization of the Alfvén modes in the Fourier space are identical to that in a full MHD system, ξR A = − k∥ × k⊥ .The compressible slow wave components, on the other hand, undergo displacements parallel to the mean magnetic field ( ξR S = k∥ , see eq. 7 and also Schekochihin et al. 2009).Henceforth, we denote the slow wave component of the velocity and magnetic fields as ⃗ u ∥ and δ ⃗ B ∥ , respectively, and their Alfvén wave components as ⃗ u ⊥ and δ ⃗ B ⊥ .We show the orthonormal reduced MHD eigenvectors for the slow and the Alfvén waves in Figure 2 with the vectors labeled ξ R S and ξ R A , respectively.1 Expanding to the first order in the reduced-MHD ordering transforms equations ( 2)-( 4) and the equation of state to (see Kawazura et al. 2022) and In the above equations, ∇ ⊥ denotes the component of the gradient operator ∇ perpendicular to the mean magnetic field B0 .Equations ( 10) and ( 11) carry the following key physical manifestations resulting from the reduced MHD ordering: (i) the Goldreich-Sridhar critical balance assumption x to become dynamically unimportant, where ⃗ x represents velocity or magnetic field fluctuations; (ii) the fast magnetosonic wave fluctuations are ordered at timescales too fast to excite, resulting in a perpendicular balance between the pressure fluctuations and the magnetic tension in the direction perpendicular to the magnetic field (δP = − ⃗ B 0 • δ ⃗ B); and (iii) the induction equation is combined with the continuity equation to capture the compressive density fluctuations (which, in turn, is related to the parallel component of the magnetic field via the perpendicular pressure balance).

The Energy Equations for Slow-and Alfvén-Waves
In order to obtain the energy equation for the Alfvén polarized fluctuations, we take the the dot product of equation ( 10) with ρ 0 ⃗ u ⊥ and of ( 11) with δ ⃗ B ⊥ and add the resulting equations.The energy equation of slow wave fluctuations can be similarly obtained by taking the the dot product of equation ( 10) with ρ 0 ⃗ u ∥ and of ( 11) with δ ⃗ B ∥ and adding the resulting equations.We then write the energy equations at the first order in the reduced-MHD approximation as for the Alfvén waves and for the slow waves.In these equations, we denote the total Alfvén-and slow-wave energies by E A and E S , respectively, where and The right hand sides of equations ( 12) and ( 13) contain the energy injection rates into the Alfvén waves (R A ) and the slow waves (R S ), respectively, and are given by ) for the Alfvén waves and for the slow waves.The first two terms in the energy injection rates R A and R S represent injection caused by the background shear and the third term represents the advection of energy caused by the Coriolis force.The left hand side of Equations ( 12) and ( 13) contains the time derivative of the respective energies and additionally, the divergence of fluxes carried by the shear velocity, the wave fluctuations, and the mean magnetic field.
Equations ( 16) and ( 17) give the expressions for the energy injection rates into the slow-and the Alfvén-wave cascades, respectively, correct to the first order in the reduced-MHD approximation.Because the slow-and Alfvén-wave cascades energetically decouple at lengthscales smaller than the injection scales of the MRI, we can infer the ratio of their relative driving power using numerical shearing-box simulations that capture the physics in the injection scales of the MRI.
In the upcoming section, we study the physics of the energy injection described in equations ( 16) and ( 17) and show that energy injection rates in the slow-and Alfvén-wave cascades can be related to quantities measured at large scales, namely the Maxwell and the Reynolds stress tensors in the flow.

PHYSICS OF ENERGY INJECTION INTO SLOW AND ALFVÉN WAVES
The injection of energy and transport of angular momentum in a local region in an accretion flow is caused by the Maxwell and Reynolds stress tensors in that region.The i jth component of the Maxwell stress tensor M i j in such a local domain of volume V in the accretion flow is defined as and the i j-th component of the Reynolds stress tensor R i j in that domain is defined as note the minus sign in the definition of the Maxwell stress tensor.In a local shearing-box configuration, the spatial volume is chosen such that it includes the energy injection scales.
In a shearing-box system, the background shear velocity field ⃗ u 0 (eq.[5]) is oriented along the azimuthal direction (ŷ) and exhibits a gradient along x.Consequently, the x − y components of the Maxwell and Reynolds stress tensors are responsible for the net injection of energy and the outward transport of angular momentum in the local shearing domain (Hawley et al. 1995;Pessah et al. 2006a).Our aim, however, is to compute the rate of energy injection into the individual slow-and Alfvén-wave cascades, and relate them to the Maxwell and Reynolds stress tensors.
The total energy injection rates into the Alfvén-and slowwave cascades can be obtained by integrating equations ( 12) and ( 13), respectively, over the spatial volume of the shearing box.Upon carrying out the integral, the terms involving the divergences vanish as they represent transport of fluxes by the background incompressible shear flow and, hence, do not contribute to energy injection.As a result, we obtain the injection powers into the individual slow-and Alfvén-wave cascades.
The slow wave injection power P S is given by and the Alfvén wave injection power P A is given by The background shear flow causes the injection of turbulent energy into the system, with the injection being captured by the terms in the square brackets in equations ( 20) and ( 21).The second term in each of equations ( 20) and ( 21) describes the effects of the Coriolis force in facilitating advection of energy in the plane perpendicular to the angular velocity of the shearing box (the x − y plane).
In order to express P S and P A in terms of components of the Maxwell and Reynolds stress tensors, we now invoke the polarizations of the slow-and Alfvén-wave fluctuations in shearing-box configurations relevant to black hole accretion disks into equations ( 20) and ( 21).The orientation of the mean magnetic field determines the polarizations of the slow-wave and the Alfvén-wave fluctuations (see discussion surrounding Figure 2).
In the disk regions of accretion flows, the radial components of the magnetic field undergo strong shearing in the azimuthal direction.As a result, the mean magnetic field in the disk is expected to be primarily along the azimuthal direction.In addition to large scale azimuthal fields, the turbulent structure in the disk is expected to develop vertical components of the magnetic fields that are determined by the local dynamics in the flow.The MRI is driven primarily by these weak vertical components of the magnetic field.Motivated by these considerations, we setup the mean magnetic field in the shearing box to lie on the y − z plane such that ⃗ B 0 = B y0 ŷ + B z0 ẑ.
Using expressions for the polarizations of the slow and the Alfvén wave fluctuations (equations [7] and [6]) under the reduced-MHD approximation, the injection rates P S and P A can be written as and respectively, where the pitch angle ϕ is defined such that tan ϕ ≡ B z0 /B y0 ≡ Γ.
These equations lend themselves to a straightforward interpretation.The azimuthal background shear flow causes the energy injection into the shearing box to be experienced by the components of the Alfvén and slow velocities and the magnetic field fluctuations along ŷ.The injection rates into these individual wave cascades are proportional to the Maxwell and Reynolds stress tensors with factors of sin ϕ and cos ϕ that capture the orientation of the wave polarizations with respect to the shear.At the same time, the effect of the Coriolis interactions is apparent on the projections of the Alfvén and slow velocity fluctuations on the x and z axes, captured by the Reynolds stress tensors.It is noteworthy to observe from equations ( 23) and ( 22) that the total energy injection rate into the combined slow and Alfvén wave cascade, as expected, depends only on the xy components of the Maxwell and Reynolds stress tensors, i.e.P S + P A = qΩ(M xy + R xy ).The x − z components of the stresses and the effects of the Coriolis force determine individual injection powers into the slow and Alfvén wave cascade but do not contribute to the total power injected into the system.
Combining expressions for the injection powers into slow and Alfvén wave cascades, we now write the slow-to-Alfvén driving power ratio P S /P A as Equation ( 24) represents the ratio P S /P A in a local domain of an accretion flow, that captures the lengthscales of energy injection into the domain.This ratio depends on the off-diagonal x − y and x − z components of the Maxwell and Reynolds stress tensors, the local shearing rate q, and the orientation of the mean magnetic field in the domain with respect to the shear (which appears as the factor Γ).
In the presence of a predominantly azimuthal mean magnetic field, the slow wave displacements are nearly along the direction of the shear (u ∥ ∝ ξ R S ∝ ŷ) and the Alfvén wave displacements lie in the plane perpendicular to the direction of the shear velocity (u ⊥ • ŷ = 0).The orientation of the slow wave displacements along the shear results in them receiving all the energy injection brought into effect by the shear.The Alfvén waves do not experience direct injection of energy via the shear velocities but only experience the Coriolis interaction via the xy component of the Reynolds stress tensor.Hence, in the limit of the magnetic field of the MRI-driven system to be purely azimuthal, relevant to the midplane regions of accretion disks, the slow-to Alfvén-wave energy injection ratio depends only on the ratio of the xy components of the Maxwell and Reynolds stresses and can be written as In principle, the ratio P S /P A can be inferred locally in global simulations of black hole accretion flows using equation (24).However, in the interest of examining the dependence of the P S /P A on the background plasma parameters, namely the plasma β and the shearing rate, we approach the problem by inferring P S /P A using the Maxwell and Reynolds stress tensors using numerical simulations of the shearingbox configuration.

NUMERICAL COMPUTATION OF SLOW-TO-ALFVÉN WAVE INJECTION RATIO
We perform a suite of magnetohydrodynamic shearingbox simulations to study the injection of slow and Alfvén waves during the linear and non-linear turbulent stages of the growth of the MRI.In the previous sections, we derived equations for the injection rates into these turbulent cascades that are correct to first order in the reduced-MHD approximation.We now infer the value of P S /P A at this order, using numerical MHD simulations that resolve the scales of the injection of the MRI.
Our approach is different from the one undertaken by Kawazura et al. (2022) to compute P S /P A .These authors performed numerical shearing-box simulations by solving the dynamical equations at first order of the reduced-MHD approximation.In contrast, we solve the MHD equations without employing the reduced-MHD approximation, but only infer the value of P S /P A using the same approximation.
Because Kawazura et al. (2020) solve the reduced-MHD equations, the energy cascades of the slow and the Alfvén waves are identically decoupled in their setup past the injection scales.In our setup, there can potentially be some coupling between the slow and the Alfvén wave cascades that are of higher order in the reduced-MHD approximation.Such higher order terms might also cause non-local exchanges in energy in the Fourier space between the velocity and magnetic field cascades (Lesur & Longaretti 2011).These effects on the slow-and Alfvén-wave driving powers in MRI-driven MHD turbulence have not been explored and will be the subject of a later study.

Numerical Setup
We carry out local MHD simulations of unstratified shearing boxes using the Godunov-type finite volume code Athena++ (Stone et al. 2020).We set up the box to be periodic in the y (azimuthal) and z (vertical) directions, and shear-periodic in the x (radial) direction (see Hawley et al. 1995, for a discussion on the shear-periodic boundary conditions).The box dimensions are set to a size (L x , L y , L z ) = (4L, 4L, L), where L = c s /Ω 0 .The resolutions of the box in the x, y, and z directions are equal, with the number of grid points being (N x , N y , N z ) = (4N, 4N, N).We perform our simulations with N = 64.In our simulation set, timeintegrations are carried out using the second order van Leer predictor-corrector methods and spatial reconstructions are performed using piecewise linear interpolation.
The MRI is primarily an instability driven by weak vertical magnetic fields in an accretion disk.The instability amplifies velocity and magnetic field perturbations in the radial direction as vertical modes.However, the direction of the mean magnetic field in the accretion disk midplane is expected to be predominantly azimuthal, as a result of the global scale shearing caused by the background azimuthal velocity field on radial magnetic fields.The strong azimuthal fields have the potential to become dynamically important, as they cause We define the vertical and azimuthal plasma β's using the respective components of the mean magnetic fields over the volume of the shearing box.We set the radial component of the mean magnetic field to zero, resulting in the mean magnetic vertical and azimuthal fields to be constant in time.
the growth of large-scale azimuthal modes2 .These considerations motivate us to set up the shearing box in our simulations with a strong magnetic field in the azimuthal direction combined with a weak magnetic field in the vertical direction.
The midplane regions of accretion disks typically exhibit shearing rates equal to the Keplerian value (q = 1.5).However, in regions away from the midplane, non-Keplerian shearing rates are expected due to non-trivial effects of pressure support arising in thick accretion flows.Additionally, in the inner regions of the flow, general relativistic corrections to the azimuthal flow profile can result in non-Keplerian shearing rates.In the interest of allowing our study to be broadly applicable to global-scale studies of accretion flows, we study the local shearing box simulations for a wide range of shearing rates.We also set the equation of state in the system to be isothermal.In Table 1, we list the parameter space of the background conditions of the shearing box explored in this paper.

Linear MRI Growth and Turbulence
The linear growth of the MRI in the presence of a weak, vertical mean magnetic field is well understood both analytically and numerically (Balbus & Hawley 1991, 1998).In a numerical simulation, the linear stage of growth is dominated by the dynamics of the fastest growing mode of the MRI.The simulations at this stage take the form of channel solutions where the velocity and the magnetic fluctuations in the xy plane undergo exponential growth with a rate equal to the growth rate of the fastest growing mode.These channel solutions are vertical, axisymmetric modes, and also satisfy the complete non-linear equations of motion.At later times of their exponential growth, the structure of the channel solutions makes them unstable to Kelvin-Helmholtz-like 0.5 1.0 1.5 2.0 Shear (q) 10 0  27) on the shearing rate q, for the fastest growing vertical mode kM = kmax = √ q(4 − q)/2.The red dotted line represents the same analytically obtained ratio, but instead assuming that the ratio in the non-linear turbulent stage is captured by an effective wavenumber keff = 0.42kmax.parasitic instabilities, which breaks the solutions into a nonaxisymmetric turbulent state (Goodman & Xu 1994;Pessah 2010).In the presence of strong azimuthal magnetic fields, large wavelength non-axisymmetric modes are also unstable to the MRI.
The outward transport of angular momentum is brought into effect by the xy components of the Maxwell and Reynolds stress tensors.Pessah et al. (2006b) showed that the ratio of the volume-integrated Maxwell to Reynolds stress tensor at late times of linear MRI growth in the presence of a weak vertical background magnetic field depends on the wavenumber k M of an unstable mode and the shearing rate (q), and is given by In equation ( 27), the quantity γ(k M ) represents the growth rate of the unstable mode with wavenumber k M , and is given by Assuming that the growth of the linear MRI is dominated by the fastest growing mode (k M = k max = √ q(4 − q)/2), the ratio of the Maxwell to the Reynolds stress tensor in this stage can be approximated as In our set of numerical simulations, we observe a period of channel-solution-like growth of the linear MRI during the early stages of each simulation.During these periods, we numerically calculate the ratio M xy /R xy .In Figure 3, we plot this ratio as a function of the shearing rate using triangular dots.We observe that the ratio is nearly independent of the plasma β but primarily depends on the shearing rate q.Additionally, we see that the ratio closely follows the analytical result for M xy /R xy for the fastest growing mode given by equation ( 28).These observations establish that the signature of the fastest growing vertical mode on the volume integrated stress ratio M xy /R xy is persistent even in the presence of significant azimuthal background magnetic fields.
In order to study the non-linear turbulent stage, we compute the ratio M xy /R xy at later times of our simulations (t = 230 − 300Ω −1 )3 .At this stage, the ratio attains values higher than the corresponding ratio at the the linear stage of growth.This behavior is observed for all values of shearing rate, independent of the background plasma β.The numerically computed ratios are shown as the circular dots in Figure 3, with the overplotted error bars showing the measured standard deviation in the ratios.While this behavior is known in shearing box systems in the presence of a weak vertical mean magnetic fields (see Figure 3 in Pessah et al. 2006b), the addition of strong azimuthal magnetic fields does not seem to alter their conclusion.
We now aim to provide a simplistic model for the sheardependence of the ratio M xy /R xy during the non-linear turbulent stage.We posit that at the late times, the dynamics of the turbulence can be captured by an effective MRIunstable mode.We model this effective mode such that its wavenumber k eff = ak max , where a is a scaling parameter and k max = √ q(4 − q)/2 is the wavenumber of the fastest growing linear mode.In Figure 3, we plot the theoretical shear dependence of the ratio M xy /R xy for a wavenumber k eff such that a = 0.42, demonstrating that the dependence of M xy /R xy on the shearing rate is well captured by this choice.
While the contribution of the azimuthal component of the mean magnetic field is subdominant in the dynamics of the MRI, it does have important implications on the injection of energy into the slow and the Alfvén wave cascades.First, as discussed in §3, the large-scale mean field in the flow determines the polarizations for the slow-and the Alfvénwave displacements.The relative orientations of the polarizations of modes and the direction of the shear is key to determining the energy injection rates into the individual cascades.Additionally, the persistent large-scale mean field in 0.5 1.0 1.5 2.0 Shear (q) 10 0 the disk causes the cascade of turbulence to be anisotropic.The strong anisotropy in the turbulence results in the decoupling of the wave cascades past the injection scales.In the following subsection, we infer the slow-to Alfvén-wave driving power from the Maxwell and Reynolds stress tensors computed here.

Slow-to Alfvén-Wave Driving Power Ratio
Having numerically computed the volume integrated Maxwell and Reynolds stress tensors in MRI driven turbulence with predominantly azimuthal mean fields, we now translate them into the slow-to-Alfvén injection power ratio using equation (25).In Figure 4, we plot the inferred P S /P A as a function of shearing rate for all the background magnetic field strengths explored in our set of simulations.
The ratio P S /P A assumes a value of near unity during the linear stage of MRI growth, for almost all the our simulations spanning various shearing rates and values of plasma β.This result can be analytically explained by combining the theoretical value of M xy /R xy for the fastest growing vertical mode (eq.[28]) with equation (25) as Remarkably, even though the ratio M xy /R xy during the linear growth of the MRI carries a strong dependence on the shear-ing rate q (eq.[28]), the value of P S /P A during this stage is independent of the shearing rate.
In the non-linear stage of the MRI, the ratio P S /P A takes values higher than in the linear stage, with a weak dependence on shear.The values of P S /P A range between ∼ 2 − 8 depending on the shearing rate.The effects of non-linear interactions in the turbulence cause the ratio M xy /R xy to be significantly higher than the corresponding linear limit at low values of q.Hence, the ratios P S /P A are significantly higher at lower values of shearing rates.
The background plasma β in the shearing box does not have a strong effect on the ratio P S /P A at a given value of a shearing rate.This is a result of the fact that the signature of MRI-driven turbulence on the ratio of Maxwell and Reynolds stress tensors, even in the non-linear stages of growth, follows the trend determined by the fastest growing linear modes (Figure 3).As discussed in the previous subsection, the ratio of the stress tensor depends only on shear and is independent of the background magnetic field strength.
For shearing rates equal to a Keplerian value, q = 3/2, which is expected in the midplanes of black-hole accretion flows, the value of P S /P A is approximately 2. This computation is in agreement with the calculation of the same quantity performed using reduced-MHD simulations by Kawazura et al. (2022), which we show as a black star in Figure 4.The ratio computed by Kawazura et al. (2022) is also independent of plasma β, in agreement with our findings and physical explanation as described in §4.
While our calculation of P S /P A is a first order estimate of the value under the reduced-MHD approximation computed from MHD simulations, the agreement with Kawazura et al. (2022) indicates that the approximation is largely valid in the parameter regime that we consider.From our analysis, it is evident that irrespective of the background plasma β and the shearing rate in an MRI driven turbulent flow, the fraction of turbulent energy injected into the slow waves is significantly higher than the energy injected into the Alfvén waves.

A Sub-grid Model for the Slow-to-Alfvén Driving
Power and Ion-to-Electron Heating Ratio In the discussion surrounding Figure 3, we suggested that the value of M xy /R xy during the non-linear stage of the turbulence can be captured by the theoretical estimate of the ratio for an effective wavenumber k eff .We demonstrated that an approximate value for k eff = 0.42k max , where k max is the wavenumber of the fastest growing vertical MRI mode.Using equation (25), we then write In equation ( 30), the value of k eff is given by q = 0.50 q = 0.75 q = 1.00 q = 1.25 q = 1.50 q = 1.75

RIAF Conditions
Figure 5.The electron heating fraction, Qe/ (Qi + Qe), inferred from eq. ( 33) as a function of plasma β, plotted for different values of the shearing rate (q).The electron heating fraction dependence on plasma β corresponding to a Keplerian shearing rate is highlighted in blue.The dashed line depicts the electron heating fraction at the limit in which the turbulence is purely driven by Alfvén waves.The typical values of plasma β encountered in GRMHD simulations of accretion flows are highlighted in the shaded region.and the corresponding growth rate γ(k eff ) for the mode k eff is given by In the dotted red lines shown in Figure 4, we depict the ratio P S /P A given by equation (30).Clearly, this expression captures the dependence of the numerically inferred values of P S /P A on the shearing rate q.We now combine this with microphysical studies by Kawazura et al. (2020) for the ionto-electron heating ratio Q i /Q e resultant from collisionless damping of slow-and Alfvén-wave turbulence to write (see eq. [1]) as The plasma β-dependent second term captures the effects of mode mixing on the heating ratio at microphysical scales comparable to the ion gyroradius.The first term represents the impact of large scale turbulent driving on the heating ratio.Notably, this impact of the large-scale driving is independent of plasma β but primarily depends on the shearing rate q.
In Figure 5, we plot the electron heating fraction, Q e / (Q i + Q e ), as a function of plasma β for different values of the shearing rate q, using equation (33).For comparison, we plot the same fraction for a purely Alfvénically driven turbulence.In the shaded region, we depict the values of plasma β that black hole accretion disks typically exhibit in GRMHD simulations.Clearly, the shearing rate-dependent 10 −1 10 0 10 1 10 2 plasma β 10 0 10 1 10 2 T i /T e q = 0 .5 0 q = 0 .7 5 q = 1 .0 0 q = 1 .2 5 q = 1 .5 0 q = 1 .7 5 RIAF Conditions Figure 6.The temperature ratio, Ti/Te, as a function of plasma β, in the inner region of a radiatively inefficient accretion flow (at ∼ 10 gravitational radii), for different values of the shearing rate, calculated using the two-temperature, covariant, semi-analytic model of Satapathy et al. (2023) and the ratio PS/PA given by equation ( 30).For typical values of the plasma β parameter and for a Keplerian shear (q = 1.5), we find Ti/Te ∼ 5 − 40, with the ratio increasing to ∼ 80 as the shear becomes shallower.nature of turbulent driving strongly determines the electron heating ratio in the disk regions of accretion flows.
The sub-grid model for the ion-to-electron heating ratio given by equation ( 33) can be integrated with twotemperature GRMHD simulations in order to model the electron thermodynamics in radiatively inefficient accretion flows.It can also be used to infer electron temperatures from one-temperature GRMHD simulations in a way which is selfconsistent with the large-scale turbulent driving.Satapathy et al. (2023) used a semi-analytic, height-and azimuthallyaveraged two-temperature model for an accretion flow to calculate the ion-to-electron temperature ratio (T i /T e ) taking into account the heat partition due to the slow-and Alfvenwave turbulent driving as well as compressional heating effects (see, eq. [37] in that paper).We use this calculation combined with the shear-dependent model Q i /Q e presented above to estimate the ion-to-electron temperature ratio in the inner accretion flow around a black hole and plot it in Figure 6 as a function of plasma β, for different values of the shearing rate.For this plot, we have set the radial distance to 10 gravitational radii and the radial profile of the radially infalling velocity to follow a power law with an index -1.5, which was obtained via calibration with GRMHD simulations by Özel et al. (2022).For typical values of the plasma β parameter and for a Keplerian shear (q = 1.5), we find T i /T e ∼ 5 − 40, with the ratio increasing to ∼ 80 as the shear becomes shallower.
The primary underlying assumption in the model for Q i /Q e given by equation ( 33) is that the mean magnetic field is nearly in the azimuthal direction in the region of interest in the accretion flow.This is a reasonable assumption for the disk regions of accretion flows.In certain configurations of global simulations, however, the disk evolves into a magnetically arrested state (MAD, see Narayan et al. 2012).In such a state, significant vertical components of magnetic field configuration might be observed in the inner regions of the disk.Additionally, the x − z components of stress tensors, which are dynamically unimportant in a local shearing box, can potentially impact the partition of energy into slow and Alfvén wave cascades due to the effects of vertical stratification.The impacts of these effects on driving the slow and Alfvén waves need to be further explored.

DISCUSSION
In this paper, we derived an analytical expression for the ratio of the slow-wave to Alfvén-wave driving powers (P S /P A ) in an MRI-driven turbulence.Working in the linear order limit of the reduced-MHD approximation, we established the connection of the ratio P S /P A to the macroscopic properties of the turbulence.We show that the ratio depends on the components of the volume integrated Maxwell and Reynolds stress tensors in the turbulence that are responsible for the outward transport of angular momentum, and on the orientation of the mean magnetic field with respect to the direction of shear.We argue that the underlying physics that determines this ratio is the impact of shear injection of energy by the MRI and the impact of Coriolis force on the projection of the slow-and the Alfvén-wave displacements with respect to the orientation of the shear.
Using numerical shearing box simulations with magnetic field configurations relevant to black hole accretion disks, we compute the ratio P S /P A during the non-linear saturated turbulent state of MRI-driven turbulence.During this stage, we find the ratio between the volume-integrated Maxwell and Reynolds stress tensors to be nearly independent of the background magnetic field strength.As a result, the ratio of the driving power of compressive slow-wave cascade to that of the Alfvén-wave cascade is nearly independent of plasma β.The only dependence seen in P S /P A is a relatively weak dependence on the shearing rate.We combine our inference of a shear-dependent model for P S /P A with hybrid kineticscale studies of ion and electron heating (Kawazura et al. 2020) to derive a local sub-grid model for the ion-to-electron heating ratio in black hole accretion flows.Two-temperature GRMHD simulations that evolve the global thermodynamics of the ions and the electrons allow for the direct study of the impact of plasma heating (Ressler et al. 2015;S ądowski et al. 2017).Recent studies, such as the ones performed by Chael et al. (2018) and Dexter et al. (2020) explore the impact of the electron heating physics on the global thermodynamics in GRMHD simulations and their implications on observational signatures.These studies show that the subgrid models for particle heating have a significant impact on not only the dynamics of the accretion flows but also on their implied images and spectra.However, in all of the above studies, subgrid heating of the plasma species due to collisionless damping has been performed assuming pure Alfvén-wave cascades in the turbulence.Our study offers a way to accurately model the ion-to-electron heating ratio considering the injection of energy into slow and Alfvén waves locally in an accretion flow based on components of the Maxwell and Reynolds stress tensors.
A number of additional dissipation mechanisms can effect the partition of heat between the ions and the electrons.Such mechanisms include episodes of magnetic reconnection (Ball et al. 2018(Ball et al. , 2019) ) and dissipation of the stress resulting from anisotropies in velocity space (Sharma et al. 2007;Kunz et al. 2018).The effects of these mechanisms and their connections to global quantities in accretion flows need to be further explored.
Using a height-and azimuthally-averaged covariant twotemperature model for an accretion disk, Satapathy et al. (2023) showed that the electron temperature in the inner regions of the disk is determined primarily by the slow waves in the turbulence.Incorporating local stress-tensor dependent models of sub-grid heating in accretion flows will make the small-scale dissipation processes self consistent with the large scale physics of energy injection.In this Appendix, we show that, even though both sets of studies discuss "compressive" and Alfvénic cascades, they actually decompose the turbulent cascades into different sets of modes.As a result, the only difference between the two sets of studies appears to be their definitions of the nature of the "compressive" cascade.Zhdankin (2021) and Bacchini et al. (2024) defined compressive fluctuations in the sense of the Helmholtz decomposition, i.e., following the literal meaning of the word compressive and distinguishing modes based on their velocity divergence.Consequently, their compressive-wave cascade consists primarily of fast wave fluctuations, which, as noted by the authors, are expected to result from spiral shock waves in the large-scale system.However, fast magnetosonic waves are ordered out in the reduced-MHD calculation of Kawazura et al. (2022) so that their label of "compressive" cascade consists purely of slow magnetosonic fluctuations.It is this ratio of this slow magnetosonic cascade to the Alfvénic one that determines the partition of heat between ions and electrons.
In this paper, we have followed the decomposition of the field into slow, fast, and Alfvén modes derived by Cho & Lazarian (2003).We use this to demonstrate the reason why the Helmholtz decomposition is not useful to separate slow and Alfvén modes in an MHD system, and to re-iterate that the nature of turbulent driving should be measured by the power in the individual slow, fast, and Alfvén wave cascades.
The Helmholtz decomposition of a velocity field in Fourier space ⃗ ṽ( ⃗ k) can be written as In the right hand side of equation (A1), the first and the second terms represent the curl-free and the divergence-free components of the vector field respectively.If the vector field represents a plasma velocity the curl-free component represents the compressible part of the flow, and the divergence-free component represents the incompressible part of the flow.
In Figure 7, we use the black arrows to represent the basis of the Helmholtz decomposition.At a given point in the Fourier space defined by the parallel (k ∥ ) and perpendicular (k ⊥ ) wavenumbers, the first term in equation (A1) picks out components of the velocity field along k.The second term on the other hand, represents components of the velocity field that lie on the plane perpendicular to k, i.e. the θ − φ plane.
We now re-examine the eigenvectors of the linear modes of MHD, namely the slow, fast, and Alfvén waves, that we introduced in §3.Following the notation of Figure 7, the orientations of the slow and the fast wave eigenvectors are represented by the blue and red arrows respectively.These eigenvectors and the basis defined by the Helmholtz decomposition are represented at points in the Fourier space such that they have the same k ∥ /k ⊥ = 1/3.Depending upon the value of α (equation 9), the slow waves are oriented between k∥ and − θ, and the fast waves are oriented between k⊥ and k.The Alfvén wave eigenvector orientation, shown in the green dot in Figure 7, is independent of α and is always in the φ direction.
We now turn to compare the eigenvectors of the linear MHD modes with the basis defined by the Helmholtz decomposition.It is to be noted that the Helmholtz decomposition only picks out the components of a velocity field along and perpendicular to k (shown by the black arrow in Figure 7), and define a solenoidal basis in the plane perpendicular to k (shown in orange).The Alfvén waves, being polarized φ, however, always lie in the plane perpendicular to k.This is consistent with the physical picture of Alfvén waves being incompressible modes.
On the other hand, both the slow and fast wave eigenvectors in general have components of velocity along k and in the plane perpendicular to k.In the limit of α → ∞, the slow waves lie entirely in the plane perpendicular to k and the fast waves are polarized along k.These properties of the orientations of the linear MHD modes clarifies that the Helmholtz decomposition in general cannot distinguish the Alfvén modes from the components of the slow and the fast waves perpendicular to k. Particularly, in anisotropic plasmas where most of the power lies in the region k ∥ ≪ k ⊥ , the fast mode carries most of the compressible power, in the k direction.The slow modes are nearly incompressible, with their polarisation nearly along − θ and for this reason, they are also referred to as pseudo-Alfvén modes (Maron & Goldreich 2001).
In the reduced-MHD approximation, most of the power is assumed to lie in the k ∥ ≪ k ⊥ region.Additionally, the fast waves are ordered out, leaving the only normal modes to be the slow and Alfvén waves.Consequently, the Helmholtz decomposition is not useful in separating out these modes and necessitates decomposing the velocity field along the slow and the Alfvén eigenvectors to analyze their injection rates.Indeed, because both the linear MRI and its linear order parasitic (Kelvin-Helmholtz) instability are fundamentally incompressible in nature, the turbulence is not expected to produce significant power in fully compressive modes (i.e., polarized along k), as found in Bacchini et al. (2024).This, however, is not in contradiction to the result of Kawazura et al. (2022), i.e., that the slow-to-Alfvénic power in the same turbulence is larger than unity.

Figure 1 .
Figure1.A representative figure of the turbulent power spectrum in black hole accretion flows (where k represents the wavenumber in units of the gravitational radius Rg).We depict the scales numerically captured by GRMHD, local shearing box, and kinetic simulations.The large scale dynamics and global angular momentum transport is captured by GRMHD simulations.The properties of turbulence such as the power spectra and the dynamics of stresses can be studied using local shearing box approaches.Past the injection scales, the slow and Alfvén wave cascades of turbulence are energetically decoupled from each other until reaching kinetic scales compared to the ion Larmor radius.At these dissipation scales, kinetic studies help in studying the partition of heating between the ions and the electrons.

Figure 3 .
Figure 3.The ratio of the volume-integrated xy component of the Maxwell stress tensor to that of the Reynolds stress tensor, plotted as a function of shearing rate, for different values of βy, computed from our simulation set.The triangles indicate the stress ratios at the stage of linear growth of the MRI (channel solutions), while the filled circles represent the same ratio computed during a saturated turbulent stage at later times in the simulation (t = 230 − 300Ω −1 ).The black dotted line represents the analytical calculation (Pessah et al. 2006b) of the dependence of Mxy/Rxy given by equation (27) on the shearing rate q, for the fastest growing vertical mode kM = kmax = √ q(4 − q)/2.The red dotted line represents the same analytically obtained ratio, but instead assuming that the ratio in the non-linear turbulent stage is captured by an effective wavenumber keff = 0.42kmax.

Figure 4 .
Figure 4.The ratio of slow wave injection rate to the Alfvén wave injection rate, PS/PA, plotted as a function of shearing rate, computed for our simulation set.The triangular dots indicate the inferred ratio at the stage of linear growth of the vertical mode MRI (channel solutions), and the circled dots represent the same ratio computed during a saturated turbulent stage at later times in the simulation (t = 230 − 300Ω −1 ).The dotted line represents the value of PS/PA calculated analytically for the linear stage using Equation (29).The black star indicates the calculation of PS/PA by Kawazura et al. (2022) who solved the dynamical equations of MHD at the first order in the reduced-MHD approximation.

Figure 7 .
Figure7.A representative plot to illustrate the difference between the decomposition of the plasma velocity field into the slow, fast and Alfvén modes, and the Helmholtz decomposition into compressive and solenoidal modes.The slow and fast wave velocity eigenvectors are shown in the blue and red arrows respectively.Their orientations are dependent on the parameter α ≡ c 2 s /v 2 A , where cs is the sound speed and vA is the Alfvén speed.The Alfvén wave eigenvector is shown the green circled dot, and is independent of α.The basis for the Helmholtz decomposition is shown for a equal value of k ∥ /k ⊥ .The solenoidal components of a vector field resulting from the Helmholtz decomposition lie in the θ − φ plane indicated in orange, while the compressive components lie along k, shown in black.

Table 1 .
Background conditions explored in shearing-box simulations