{'fraction_non_alphanumeric': 0.02916814357405369, 'fraction_numerical': 0.01567076469600582, 'mean_word_length': 4.564858551777835, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A simple framework for modelling the dependence of bulk Comptonization by turbulence on accretion disc parameters

Warm Comptonization models for the soft X-ray excess in AGN do not self-consistently explain the relationship between the Comptonizing medium and the underlying accretion disc. Because of this, they cannot directly connect the fitted Comptonization temperatures and optical depths to accretion disc parameters. Since bulk velocities exceed thermal velocities in highly radiation pressure dominated discs, in these systems bulk Comptonization by turbulence may provide a physical basis in the disc itself for warm Comptonization models. We model the dependence of bulk Comptonization on fundamental accretion disc parameters, such as mass, luminosity, radius, spin, inner boundary condition, and $\alpha$. In addition to constraining warm Comptonization models, our model can help distinguish contributions from bulk Comptonization to the soft X-ray excess from those due to other physical mechanisms, such as absorption and reflection. By linking the time variability of bulk Comptonization to fluctuations in the disc vertical structure due to magnetorotational instability (MRI) turbulence, our results show that observations of the soft X-ray excess can be used to study disc turbulence in the radiation pressure dominated regime. Because our model connects bulk Comptonization to one dimensional vertical structure temperature profiles in a physically intuitive way, it will be useful for understanding this effect in future simulations run in new regimes.


INTRODUCTION
The soft X-ray excess in active galactic nuclei (AGN) refers to the part of the spectrum below 1keV that lies above the power law fit to the hard (2-10keV) X-rays (Singh et al. 1985;Arnaud et al. 1985). The origin of this component is not well understood. In narrow-line Seyfert 1s (NLS1), which are high Eddington, comparatively low mass (M ∼ 10 6 ) AGN, the disc color temperature is high enough that it contributes to the soft excess, but it cannot fully explain it (Done et al. 2012, hereafter D12). In other systems, the soft excess must entirely be due to some other physical mechanism.
The soft excess is often fitted with warm Comptonization models (D12; Magdziarz et al. 1998;Janiuk et al. 2001 Pal et al. 2016;Petrucci et al. 2017). In this picture, photons from the disc are upscattered by a warm (∼ 0.2keV), optically thick medium. While these models have had considerable success fitting the data, they do not explain the physical origin of such a medium in some sources and therefore cannot predict its properties from the parameters of the underlying accretion disc. Models more tied to intrinsic disc physics include reflection models (Crummy et al. 2006;Ross & Fabian 2005) and absorption models (Gierlinski & Done 2004;Schurch & Done 2007).
In this work we explore the extent to which Comptonization by turbulent, bulk velocities in the disc itself, rather than thermal velocities, can provide a physical basis for warm Comptonization models. Since the ratio of turbulent velocities to thermal velocities for the alpha disc model (Shakura & Sunyaev 1973 we expect Comptonization by turbulence to be relevant in radiation pressure dominated systems. The ratio of radiation to gas pressure increases with mass and accretion rate, so we expect this effect to be most significant in high Eddington accretion discs around supermassive black holes, such as NLS1s. It is likely insignificant, therefore, in lower Eddington flows, even though the data show that in these flows the soft excess carries a more significant fraction of the power Mehdipour et al. 2011Mehdipour et al. , 2015. Preliminary work on bulk Comptonization in accretion discs was done by Socrates et al. (2004) and Socrates (2010). Kaufman & Blaes (2016) studied the physical principles that govern bulk Comptonization in detail. Kaufman, Blaes, and Hirose (2017) (hereafter, K17) used data from radiation magnetohydrodynamic (MHD) simulations to model bulk Comptonization for accretion discs with parameters typical of NLS1s. They characterized their results with temperatures and optical depths to make contact with warm Comptonization models of the soft excess. For a system with disc parameters similar to those fit by D12 to the source REJ1034+396, the temperatures, optical depths, and Compton y parameters found by K17 broadly agree with those measured by D12.
This work simplifies and generalizes the bulk Comptonization model presented in K17 in order to develop greater physical insight into this process and explore a larger space of accretion disc parameters. In K17, bulk Comptonization is modeled by fitting the Comptonization temperature and optical depth parameters to spectra computed with Monte Carlo post-processing of simulations of the turbulence. Here we develop a procedure to infer these Comptonization parameters from the underlying disc vertical structure radiation MHD simulation data without computing spectra. The immediate benefit of this is that we can efficiently explore a larger space of accretion disc parameters. We can also find the time-averaged Comptonization parameters for a given simulation since without computing spectra we can now efficiently calculate Comptonization parameters for multiple timesteps.
More importantly, our model provides a physically intuitive framework for understanding bulk Comptonization. In particular, we show that the variation of this effect with disc parameters can be understood in terms of the vertical gas temperature and "wave" temperature (Section 2.1) profiles. This allows us to determine the dependence of bulk Comptonization on each accretion disc parameter separately without exhaustively exploring a multiparameter space. We can also probe how various physical effects, such as vertical radiation advection (Blaes et al. 2011;Jiang et al. 2013Jiang et al. , 2014, may impact bulk Comptonization, as well as evaluate how robust our model's predictions are to changes in the disc vertical structure. In particular, although the specific bulk Comptonization parameters that we calculate here result from applying our model to limited, scaled data from the shearing box simulation 110304a (K17), we show that our principal findings regarding how bulk Comptonization scales with fundamental accretion disc parameters are likely to be robust to differences in the disc vertical structure seen in other simulations. Furthermore, understanding this framework should be useful for developing physical intuition in new situations in which some of our particular results may no longer hold, such as shearing box or global disc simulations run in radically different regimes.
The structure of this paper is as follows. In Section 2 we describe our model and show why it is effective. In Section 3 we apply our model to data from radiation MHD simulations. We show how the dependence of bulk Comptonization on shearing box parameters can be understood in terms of one dimensional temperature profiles (Section 3.2), and then proceed to examine its dependence on each accretion disc parameter individually (Section 3.3). We estimate bulk Comptonization for an entire disc as well by fixing the radius to the region of maximum luminosity (Section 3.4). We consider the effect of including radiation advection (Section 3.5) and discuss the time variability of bulk Comptonization within a given simulation (Section 3.6). In Section 4 we discuss our results, and we summarize our findings in Section 5.

Overview
Since this work simplifies and generalizes the bulk Comptonization model that we presented in K17, we begin by summarizing how they calculated the bulk Comptonization temperature and optical depth for a given system. At each radius in an accretion disc, they used radiation MHD stratified shearing box simulation data to calculate spectra both including and excluding velocities. Each spectral computation was performed by running a post-processing Monte Carlo simulation on a simulation data snapshot at a particular epoch in time. Since their data was limited, they developed a scheme to scale the original data to the accretion disc parameters of interest. At each radius in the disc, they used the Kompaneets equation to pass the spectrum computed without velocities through a Comptonizing medium with a given electron temperature and optical depth, and the resulting spectra were superposed to obtain the observed spectrum. Meanwhile, the spectra computed with velocities at each radius were superposed to obtain a different observed spectrum. The temperature and optical depth parameters (which are assumed to be the same at all radii) were adjusted until the two spectra match.
Here we develop a more efficient and physically revealing procedure for calculating the Comptonization temperature and optical depth. We focus on computing these parameters for each radius individually rather than for the whole disc at once, a choice that also allows us to study the dependence of bulk Comptonization on radius. This choice does not limit us to studying individual radii, since bulk Comptonization for a whole disc can be estimated by the Comptonization parameters at the radius where the luminosity is greatest (Section 3.4). One other difference from K17 is that we include only turbulent velocities, not shear velocities. This allows us to study the effects of turbulence alone. The preliminary results in K17 suggest that bulk Comptonization by shear is subdominant to bulk Comptonization by turbulence, but to rigorously calculate the effect of shear near the photosphere will require global simulations. Since the scalings for the shear velocities are nearly identical to the scalings for the turbulent velocities (K17), this omission should not affect our conclusions regarding the general dependence of bulk Comptonization on accretion disc parameters.
At a given radius, we define the Kompaneets parameters, consisting of the electron temperature TK and optical depth τK, analogously to how the temperature and optical depth are defined for the whole disc in K17. In other words, spectra are calculated with and without velocities at a given radius, and the spectrum computed without velocities is passed through a Comptonizing medium using the Kompaneets equation. The temperature and optical depth of the medium are adjusted until the resulting spectrum matches the spectrum calculated with velocities. To scale simulation data to different accretion disc parameters we use the scheme developed in K17.
We show that the Kompaneets parameters are approximated by what we will henceforth refer to as the Comptonization parameters. We define these parameters with an efficient and physically revealing procedure that we outline here and then discuss in greater detail in the following sections. First, we map the bulk velocity grid to a temperature grid by defining at each point a bulk "wave" temperature, where (∆v) 2 r is the average square velocity difference between subsequent photon scatterings at r. We note that if instead of applying equation (2) to the bulk velocities in a region we apply it to the thermal velocity distribution at a particular point, then since (3) we find that 3 2 kBTw = 1 2 me v 2 r , as expected. We call this a "wave" temperature and not a "turbulent" temperature because it depends on the power spectrum of the turbulence and is less than the temperature that one usually associates with a turbulent velocity distribution, which is given by 3 2 kBT = 1 2 me v 2 r , analogous to a thermal temperature. We discuss (∆v) 2 r in more detail in Section 2.2.2. Next, we horizontally average all simulation variables, including the newly defined wave temperature, to obtain 1D profiles of the data. For the systems of interest, the wave temperature is negligible compared to the gas temperature at the effective photosphere, and it increases going outward so that near the scattering photosphere it may exceed it. We define the Comptonization optical depth τC as the optical depth of the region in which the wave temperature is at least half the gas temperature, a region that we will henceforth refer to as the bulk Comptonization region. We define the Comptonization temperature TC as a weighted average of the sum of the gas and wave temperatures in this region, given by (7) In the next section we both describe in detail why this procedure approximates the Kompaneets parameters and demonstrate its effectiveness by comparing what it predicts with the results of actual Monte Carlo spectral calculations.

Physical justification for the bulk Comptonization model
To justify our definition of the Comptonization parameters, τC and TC, we start with the procedure that defines the Kompaneets parameters and then incrementally simplify it.
In the following sections we detail each step of this process. Where appropriate we invoke an accretion disc parameter set for which M = 2 × 10 6 M⊙ and L/L Edd = 5 as a test case. The other parameters are given in Table 1. The parameter a is the black hole dimensionless spin parameter, ∆ǫ is the change in efficiency for a non-zero torque inner boundary condition (Agol & Krolik 2000), α is the ratio of vertically integrated stress to vertically integrated pressure, and α0 is the value of α for the original simulation data. We list α/α0 rather than α since it is the former that we can directly adjust with the scaling scheme from K17. Typically α0 ∼ 0.01. We note that these parameters are nearly identical to those of the systems modeled in K17, and were originally chosen to correspond to those fit by D12 to the NLS1 REJ1034+396. The only parameter whose value differs from the value in K17 is L/L Edd , which we set here to 5 rather than 2.5. Since bulk Comptonization increases with L/L Edd , we choose a higher value here so that it is easier to see the effectiveness of our approximations in plots of actual spectra. After describing all steps, we demonstrate the effectiveness of the resulting procedure at six different radii for not only the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set but for two others as well. One is the same except that L/L Edd = 2.5, the original value in K17. For the other, M = 2 × 10 8 M⊙ and L/L Edd = 4.2. All disc parameter sets are given in Table 1. As in K17, all spectra in this section are computed with Monte Carlo post processing simulations (Davis et al. 2009;Pozdnyakov et al. 1983), using data from the 140 orbits timestep of simulation 110304a. We discuss this simulation in more detail in Section 3.2.2.

2.2.1
Step 1 -Truncate the simulation data inside the effective photosphere and turn off emissivity above this surface To begin, we observe that we can modify the defining procedure for calculating the Kompaneets parameters by using simulation data that is truncated inside the effective photosphere and turning off the emissivity everywhere except at this surface, without changing the resulting parameters. The effective photosphere is defined by using the Planck mean opacity. This phenomenon was demonstrated in the context of fitting Kompaneets parameters for an entire disc at once  (Table 1). The green and cyan curves are computed with velocities, and the blue and red curves are computed without velocities. For the cyan and red curves, the simulation data is truncated inside the effective photosphere and the emissivity is set to zero everywhere except at the effective photosphere.
in K17, but it arises from the fact that it is true for individual radii.
For example, we calculate spectra both with and without velocities at r = 14 for the M = 2×10 6 M⊙, L/L Edd = 5 parameter set (Table 1). We note that all numerical radii in this paper are in units of the gravitational radius GM/c 2 of the black hole. We also calculate spectra using simulation data that is truncated at the effective photosphere and in which the emissivity is set to zero everywhere except the base. All four resulting spectra are plotted in Figure 1. We see that except at very low energies, the two spectra calculated with velocities coincide and the two spectra calculated without velocities coincide.
Two reasons underlie this result. First, the emergent spectrum is dominated by photons originally emitted at or near the effective photosphere since free-free emission depends strongly on density. Second, bulk Comptonization is negligible except near the scattering photosphere (as we show in Section 2.2.2) and so its effect does not depend on the precise effective optical depth at which photons are originally emitted. Because all systems in this work meet these conditions, this result is robust. Therefore, modifying the defining procedure for calculating the Kompaneets parameters in this way has a negligible effect on the outcome.
Calculating the temperature and optical depth with this modified method is the first step to simplifying the calculation and developing physical insight into the problem. Not only does this modified procedure run faster, but it shows that in order to model and understand bulk Comptonization we only need to understand the effect of bulk Comptonization on photons emitted at the effective photosphere.

2.2.2
Step 2 -Map the velocity grid to a "wave" temperature grid Next, we map the bulk velocity grid to a temperature grid by defining at each point a bulk "wave" temperature given by equation (2). This definition is motivated by the results of Section 4.2 in K16. There they showed that for a periodic box with statistically homogeneous turbulence and an escape probability, bulk Comptonization can be treated as thermal Comptonization by solving the Kompaneets equation with a temperature given by 3 2 kBTw = 1 4 me (∆v) 2 , where (∆v) 2 is the volume average of (∆v) 2 r , which is in turn given by K16 equation (56): (∆v) 2 r = (∆v(∆r, r)) 2 P∆r(∆r)d 3 ∆r.
Here, ∆v(∆r, r) is the velocity difference between positions r and r + ∆r, and P∆r(∆r) is the probability density that a photon scattering at r subsequently scatters at r+∆r. Hence (∆v) 2 r is the average square velocity difference between subsequent photon scatterings at r. In this work, therefore, equation (2) defines a wave temperature at each point instead of taking a volume average and defining it for an entire box.
To develop physical intuition into equation (2) it is important to understand the dependence of (∆v) 2 r on density. In the high density limit the velocity difference between subsequent photon scatterings is small. In particular, (∆v) 2 r is proportional to the square of the mean free path λp (K16) so Tw decreases significantly with increasing density. In the low density limit, on the other hand, (∆v) 2 r approaches 2 v 2 so that 3 2 kBTw approaches 1 2 me v 2 r . We also define the bulk temperature, given by Applying equation (8) directly to simulation data is somewhat problematic, so we discuss our implementation in detail in Appendix A. For example, in Figure 2 we plot the profile of the density weighted horizontal average of the wave temperature at r = 14 for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table 1). We also plot the gas temperature, Tg, the bulk temperature, T bulk , and the sum of the gas and wave temperatures. The dashed line on the right denotes the location of the scattering photosphere, which we define as the height at which the Thomson optical depth τs = 1. We see that the wave temperature significantly increases with decreasing density (that is, moving rightward) as the photon mean free path grows, and is comparable to the bulk temperature only near the scattering photosphere.
We find that photon spectra computed with simulation data in which the velocities are turned off and the wave temperatures are added to the gas temperatures approximate photon spectra computed with the velocities turned on. For example, in Figure 3 we plot spectra at r = 8. 5, 9.5, 11, 14, 20, and 30 for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table 1), computed with and without velocities. We also plot spectra computed with data in which the velocities are turned off and the wave temperatures are added to the gas temperatures. We see that these spectra approximate the spectra computed with velocities. In other words, bulk  (Table 1) for the gas temperature Tg (blue), bulk temperature T bulk (green), wave temperature Tw (red), and sum of gas and wave temperatures (cyan). The dashed lines denote where τs = 1 and τs = 10.
Comptonization can be modeled by thermal Comptonization in which the temperature is given by equation (2).

2.2.3
Step 3 -Horizontally average the simulation data Once bulk velocities are replaced by wave temperatures, it is straightforward to further simplify the problem by horizontally averaging the simulation data. In order that the effects of bulk Comptonization remain unchanged, the wave temperature data must be density averaged, not volume averaged. For example, we calculate the spectrum at r = 14 for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table  1) using data in which the velocities are turned off and the wave temperatures are added to the gas temperatures, as described in Section 2.2.2. We repeat this calculation using horizontally, density weighted averaged data and plot both spectra in Figure 4. We see that the two spectra coincide. In Figure 5 we again plot the spectrum calculated with the unaveraged data as well as a spectrum calculated with horizontally averaged data, except that this time the wave temperatures are computed with a simple spatial horizontal average instead. We see that simple spatially averaging the wave temperatures with no density weighting overestimates bulk Comptonization.
Density averaging improves the accuracy because the time photons spend in a region increases with the region's density. The reason that volume weighting overestimates bulk Comptonization is that the wave temperature is strongly correlated with density. As we discussed in Section 2.2.2, the wave temperature decreases with density, and so volume averaging gives too much weight to regions where the wave temperature is larger. Because the gas temperature, on the other hand, is not strongly correlated with density, horizontally volume weighting the gas temperature has a negligible impact on the spectrum. For example, we calculate another spectrum with horizontally averaged data, except that this time the gas temperatures are computed with a simple spatial average instead. We plot the result alongside the spectrum calculated from the unaveraged data in Figure  6. We see that the two spectra coincide, which indicates that gas temperature inhomogeneities at a given height are not sufficiently correlated with density inhomogeneities to affect the spectrum.
The fact that we can use horizontally averaged quantities means we can map a three dimensional problem to a one dimensional problem, an important step to efficiently calculating and understanding bulk Comptonization. Plots of horizontally averaged quantities such as Figure 2, first introduced in Section 2.2.2, will be of great use in the remainder of this work.

2.2.4
Step 4 -Solve the 1D inhomogeneous thermal Comptonization problem By this point, we have modified the original procedure for calculating the Kompaneets parameters by instead calculating spectra with and without adding the wave temperature profile to the gas temperature profile (Section 2.2.2), using horizontally averaged data (Section 2.2.3) truncated inside the effective photosphere in which velocities are turned off and emission is zeroed everywhere except at the effective photosphere (Section 2.2.1). The temperature and optical depth are adjusted until the spectra match.
To further simplify the problem, we first need to understand the effect of thermal Comptonization on photons emitted at the base of an inhomogenous one dimensional medium. We expect that if the optical depth is not too high, so that the average photon energy is always significantly below the local temperature (i.e. the photon spectrum does not saturate), then this process can always be well described by a homogeneous thermal Comptonization model. Since the number of scatterings is proportional to the square of the optical depth, the appropriate average scattering temperature should be given by We test this description of 1D thermal Comptonization at r = 8. 5, 9.5, 11, 14, 20, and 30 for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table 1). At each radius, we add the wave temperature to the gas temperature and truncate the data inside τs = 10, an optical depth that is large enough to result in significant Comptonization but small enough to prevent the saturation of photon spectra for our purposes. We place a 50eV Planck source at this location and calculate the emergent spectra. At each radius we also use the Kompaneets equation to pass the source through a homogeneous medium with temperature T1D and optical depth τ = 10. We plot the resulting spectra for r = 14 in Figure 7. Spectra at the other radii illustrate the same effect and are plotted in Appendix C ( Figure C1). The value of T1D at each radius is given in Table 2. We see that the spectrum calculated using the Kompaneets equation coincides with the spectrum computed directly from the data, confirming that unsaturated, 1D inhomogeneous thermal Comptonization is well modeled by homogeneous thermal Comptonization, with the temperature given by T1D. . Normalized accretion disc spectra at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1) computed with (green) and without (blue) the velocities. For the red curve, the velocities were not included but the wave temperatures were added to the gas temperatures.

2.2.5
Step 5 -Use the solution to the 1D thermal Comptonization problem to model bulk Comptonization Armed with the results of Section 2.2.4, we return to the original problem. We begin by observing that typically in the region between the effective photosphere and the scat-tering photosphere (that is, where τ eff < 1 and τs > 1), the gas temperature does not vary significantly. The wave temperature, on the other hand, changes rapidly with density. It is negligible compared to the gas temperature at the bottom of the effective photosphere and increases moving outward. Near the scattering photosphere it may exceed the gas temperature, depending on the parameters of the problem.  (Table 1) computed by omitting the velocities and instead adding the wave temperatures to the gas temperatures. The blue curve is computed with unaveraged data, and the green curve is computed with horizontally density weighted averaged data.  (Table 1) computed by omitting the velocities and instead adding the wave temperatures to the gas temperatures. The blue curve is computed with unaveraged data, and the green curve is computed with horizontally averaged data. For the latter, the wave temperature averages are volume weighted and the other variable averages are density weighted.
The effect of adding the wave temperature profile to the gas temperature profile, therefore, is to take the spectrum that results from when there is no wave temperature and pass it through a Comptonizing medium of optical depth given by that of the region where the wave temperature is comparable to the gas temperature. We define this to be the region in which the wave temperature is at least half the gas temperature. We refer to it as the bulk Comptonization region and define the Comptonization optical depth parameter τC to be its optical depth. We then define the associated Comptonization temperature parameter TC by equation (10), where T is the sum of the gas and wave temperatures.  Figure 6. Normalized accretion disc spectra at r = 14 for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1) computed by omitting the velocities and instead adding the wave temperatures to the gas temperatures. The blue curve is computed with unaveraged data, and the green curve is computed with horizontally averaged data. For the latter, the gas temperature averages are volume weighted and the other variable averages are density weighted. Table 2. Values of T 1D for vertical structure data truncated at τs = 10 at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1). . Normalized spectrum (red) computed by passing a 50eV Planck source (blue) through vertical structure data truncated at τs = 10 at r = 14 for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1). The velocities are zeroed and the wave temperatures are added to the gas temperatures. The green curve is calculated by using the Kompaneets equation to pass the 50eV Planck source through a homogeneous medium with temperature T 1D , given in Table 2. Spectra at other radii are plotted in Appendix C ( Figure C1).  Figure 8. Normalized spectra at r = 14 for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1) computed with (red) and without (blue) velocities. The green curve is calculated by using the Kompaneets equation to pass the blue curve through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 3. Spectra at other radii are plotted in Appendix C ( Figure C2).
For example, we calculate spectra with and without velocities at r = 8. 5, 9.5, 11, 14, 20, and 30 for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table 1). At each radius we also use the Kompaneets equation to pass the spectrum computed without velocities through a homogeneous medium with temperature TC and optical depth τC. We plot the resulting spectra for r = 14 in Figure 8. Spectra at the other radii illustrate the same effect and are plotted in Appendix C ( Figure C2). The temperature and optical depth parameters at all radii are given in Table 3. We see that the spectrum computed with the Kompaneets equation approximates the spectrum computed with velocities, indicating that the Comptonization parameters TC and τC approximate the Kompaneets parameters, TK and τK.
To demonstrate that the effectiveness of the parameters TC and τC at describing bulk Comptonization is not limited to a narrow mass range, we modify the parameter set by choosing a significantly higher mass, M/M⊙ = 2×10 8 (Table  1). We again calculate spectra at multiple radii, with and without velocities, and plot the results for r = 14 in Figure  9. In the same figure we plot the spectrum predicted by the parameters TC and τC for r = 14, given in Table 4, and see that the resulting spectrum well approximates the spectrum computed with velocities. Spectra at all radii are plotted in Appendix C ( Figure C3).
Finally, we show that the Comptonization parameters TC and τC well describe bulk Comptonization for the M = 2 × 10 6 M⊙, L/L Edd = 2.5 parameter set (Table 1), whose value of L/L Edd is the same as in K17. The corre- Table 4. Comptonization temperatures and optical depths at multiple radii for the M = 2 × 10 8 M ⊙ , L/L Edd = 5 parameter set (Table 1) (Table 1) computed with (red) and without (blue) velocities. The green curve is calculated by using the Kompaneets equation to pass the blue curve through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 4. Spectra at other radii are plotted in Appendix C ( Figure C3).
sponding spectra for r = 14 are plotted in Figure 10, and the Comptonization parameters are given in Table 5. We see that the spectrum predicted by the parameters TC and τC well approximates the spectrum computed with velocities. Spectra at all radii are plotted in Appendix C ( Figure C4).

2.2.6
Step 6 -Qualify the bulk Comptonization model scope In Sections 2.2.1-2.2.5, we justified each step of a process that results in a simplified, physically revealing model for bulk Comptonization by turbulence, and demonstrated the success of this model for multiple radii, masses, and accretion rates. We acknowledge that underlying each step are various assumptions, some of which may not hold over the entire range of accretion disc parameters of interest. This is a limitation only if the sole goal is to reproduce the Kompaneets temperature and optical depth as they are originally defined. But in this work our primary goal is rather to characterize bulk Comptonization in a physically revealing way so that we can easily map out its dependence on Table 5. Comptonization temperatures and optical depths at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 2.5 parameter set (Table 1).  (Table 1) computed with (red) and without (blue) velocities. The green curve is calculated by using the Kompaneets equation to pass the blue curve through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 5. Spectra at other radii are plotted in Appendix C ( Figure C4). a wide range of disc parameters and understand how this dependence itself may change depending on the robustness of certain features in the disc vertical structure. Therefore, each step of this process should be viewed more as a search for parameters that are physically revealing and easily calculated rather than as an attempt to merely speed up the calculation of the Kompaneets parameters. For example, we may find that in some regimes the calculated Comptonization region optical depth is sufficiently large that photon spectra saturate, which violates an assumption we made in Section 2.2.4. In this case, the Comptonization temperature and optical depth will probably differ somewhat from the Kompaneets parameters. But since they would still, by definition, tell us the optical depth of the region in which bulk Comptonization is significant as well as the weighted sum of the gas and wave temperatures in this region, they would still provide a useful characterization of bulk Comptonization.

Overview
The independent variables in radiation MHD shearing box simulations are the surface density Σ, the vertical epicyclic frequency Ωz, and the strain rate ∂xvy. Since our simulation data is limited, we use the scheme developed in K17 to scale data from one set of independent variables to another. This scheme also allows for the variation of α (Shakura & Sunyaev 1973), defined as the ratio of the vertically integrated total pressure to the vertically integrated total stress.
In Section 3.2 we calculate the dependence of bulk Comptonization on the four shearing box parameters Σ, Ωz, ∂xvy, and α. In particular, we show that these four parameters can in practice be reduced to two parameters, Σ and α 3 Ωz, so that the dependence of bulk Comptonization on shearing box parameters can be illustrated in a single figure with multiple curves. In Section 3.3 we show the dependence of bulk Comptonization on accretion disc mass, luminosity, radius, spin, and inner boundary condition. To do this, we examine the dependence of the shearing box parameters Σ and Ωz on these parameters. In Section 3.4 we estimate bulk Comptonization for an entire disc by setting the radius equal to the value that contributes maximally to the luminosity. The scaling scheme from K17 assumes that the radiation energy flux is carried by radiation diffusion, but it does allow for variation in the opacity κ. In Section 3.5 we show that vertical radiation advection can be included indirectly by varying κ, and we examine the effect of this on bulk Comptonization. Finally, in Section 3.6 we examine how bulk Comptonization is effected by time variability in the simulation data.

Reduction of four shearing box parameters to two
To simplify the problem, we observe that for Newtonian disc scalings For Kerr disc scalings, the strain rate scale factor is nearly equal to the vertical epicyclic frequency scale factor (K17).
To show this, in Figure 11 we plot the ratio of these quantities for multiple values of black hole spin. We see that at worst they agree to within ∼ 6%. For the purpose of understanding bulk Comptonization, then, we can set these factors equal to each other. Next, we eliminate the dependence on α by proving the following statements for the bulk Comptonization parameters TC, τC, and yp,C (the Compton y-parameter): For any To prove equation (12), we first show that for fixed Σ the wave temperature scales as the bulk temperature. The wave temperature depends not only on the bulk temperature but on the density via the ratio of each length scale in the turbulence λw to the photon mean free path λp (K16). The scalings for these parameters are λw ∼ h and λp ∼ ρ −1 , where h is the disc scale height. Since for fixed Σ the scaling for the density is given by ρ ∼ h −1 , varying only other parameters leaves the ratio λw/λp unchanged and hence the wave temperature scales as the bulk temperature.
Next, we observe that the ratios of the turbulent velocities to the thermal velocities in the midplane and the photosphere are given by and respectively. We see that if we vary the vertical epicyclic frequency inversely to α 3 , then the ratio of the turbulent kinetic energy to the thermal kinetic energy remains unchanged everywhere. Since for fixed Σ the wave temperature is proportional to the gas temperature, it follows that the ratio of the wave temperature to the gas temperature is also everywhere unchanged. Then, since τC is defined as the optical depth of the region in which the wave temperature is comparable to the gas temperature (Section 2), under these circumstances it can change only if the overall density or the scale height changes. That is, But ρch ∼ Σ, and Σ is held constant since it is an independent variable, so ρc varies inversely to h. In other words, for fixed Σ the optical depth for any length scale is invariant. (Note that this is the exact same reason that the wave temperature scales as the bulk temperature.) It follows that Equation (12) follows directly from equation (18).
To prove equation (13), we start with the turbulent and thermal velocity scalings individually (K17), rather than the ratio of scalings: We observe that if we vary the vertical epicyclic frequency inversely to α 3 , as we just showed we must do in order to leave τC unchanged, then the scalings for the individual variables are We see that all velocities scale inversely to α. Since the Comptonization temperature is defined as a density weighted average of the sum of the wave and gas temperatures (Section 2), it follows that Equation (13) follows directly from equation (25). Finally, the definition of the Compton y-parameter is where N is the average number of scatterings. For a plane parallel geometry with τ > 1, 1 N = 1.6τ 2 , so and equation (14) follows directly from equations (12) and (13). Therefore, for the purpose of understanding bulk Comptonization, we can regard Σ and α 3 Ωz as the fundamental shearing box parameters and αTC, τC, and αyp,C as the Comptonization parameters.

Dependence of bulk Comptonization on Σ and Ωz
The original data we use is from ZEUS simulation 110304a (K17). The shearing box parameters for this simulation are given in Table 6. The time-averaged α parameter is α0 = 0.01 (which we do not list in Table 6 since it is not an independent variable). These correspond to an accretion disc annulus with parameters given in Table 7. The opacities included are electron scattering and free-free. These should be good approximations for the opacities in AGN in the near and super-Eddington regimes of interest in this paper. At a given timestep, we scale the data to a range of values of Σ and Ωz with the scheme in K17, and then calculate the resulting Comptonization parameters TC, τC, and yp,C with the procedure detailed in Section 2. We repeat this for 21 timesteps spaced 10 orbital periods apart and plot the time-averaged results in Figure 12. The error at each point is estimated by dividing the sample standard deviation by the square root of the number of timesteps. For clarity we add the subscript "fid" to the shearing box parameters of the fiducial system given in Table 8 to distinguish them from the original simulation parameters which we denote by the subscript "0". The fiducial shearing box parameters correspond to r = 20 for the M = 2 × 10 6 M⊙, L/L Edd = 2.5 parameter set (Table 1). Following K17, we chose these parameters to be similar to those fit to the NLS1 REJ1034+396 by D12.  We can intuitively understand these results by looking at equations (19)-(21) and the horizontally averaged temperature profiles shown in Figure 13. In particular, we show below why τC and TC strongly increase with increasing Σ −1 , while τC increases weakly and TC decreases weakly with increasing Ω −1 z . It then follows that since yp,C depends more strongly on τC than on TC, yp,C increases strongly with increasing Σ −1 and weakly with increasing Ω −1 z . Since yp is generally used as a proxy for the overall magnitude of Comptonization, we conclude that bulk Comptonization increases strongly with increasing Σ −1 and weakly with increasing Ω −1 z . Because of this as well as the fact that Ω −1 z ∼ M and Σ −1 ∼ L/L Edd (which we discuss in Section 3.3), we treat Ω −1 z and Σ −1 as the fundamental parameters rather than Ωz and Σ.
Dependence of τC on Ω −1 z We start by considering the dependence of the bulk Comptonization optical depth τC on Ω −1 z . Equations (19)- (21) show that the result of varying Ω −1 z is to multiply the entire gas temperature profile by a constant and leave the bulk temperature profile unchanged. Since the wave temperature scales as the bulk temperature for fixed Σ −1 (Section 3.2.1), the wave temperature profile is also unchanged. Increasing Ω −1 z , therefore, corresponds to moving the gas temperature profile downward in Figure  13, increasing the optical depth of the region in which the wave temperature is comparable to the gas temperature, consistent with the results shown in Figure 12.
To understand the dependence of the Comptonization temperature TC on Ω −1 z , we first need to look at the bulk Comptonization region weighted average gas temperature TC,g and wave temperature TC,w, individually. We plot the time-averaged dependence of these two parameters on Ω −1 z in Figure 14. As we already showed, increasing Ω −1 z moves the gas temperature profile downward in Figure 13. Since the gas temperature profile does not spatially vary significantly in this region, equations (20) -(21) imply that TC,g will be approximately proportional to Ω 1/4 z . Typically the gas temperature profile is slightly decreasing at the lower boundary of this region so that the dependence is slightly shallower than Ω 1/4 z , which is what we find in Figure 14. Since increasing Ω −1 z has no effect on the wave temperature profile, the only effect of increasing Ω −1 z on TC,w Table 8. Fiducial shearing box parameters, corresponding to r = 20 for the M = 2 × 10 6 M ⊙ , L/L Edd = 2.5 parameter set (Table  1). . Horizontally averaged temperature profiles for the fiducial shearing box parameters, given in Table 8, at the 20 orbits timestep. The dashed line denotes where τs = 1. α/α fid kT C,w (eV) Figure 14. Dependence of Comptonization gas and wave temperatures on shearing box parameters. The blue, green, red, and cyan curves correspond to (Σ/Σ fid ) −1 = 1, 2, 3.3, and 5, respectively.
is to decrease the height of the lower boundary of the bulk Comptonization region. As it decreases, TC,w also decreases since not only does the wave temperature profile decrease with increasing optical depth, but equation (7) gives greatest weight to Tw in the region where the optical depth is the largest. Therefore, TC,w also decreases with increasing Ω −1 z but less so than TC,g, as we see in Figure 14. Since TC = TC,g + TC,w, TC also decreases with increasing Ω −1 z at a rate slightly faster than TC,w but slower than TC,g, as we see in Figure 12.
Dependence of τC on Σ −1 From equations (19)-(21) we see that v 2 turb depends much more strongly on Σ −1 than does v 2 th , so decreasing Σ −1 moves the bulk temperature profile downward in Figure 13 relative to the gas temperature profile. The wave temperature profile moves downward even more than the bulk temperature profile, since the wave temperature depends on the velocity difference between scatterings, which decreases with decreasing Σ −1 . Therefore, the Comptonization optical depth decreases with decreasing Σ −1 . We expect the dependence of τC on Σ −1 to be much stronger than its dependence on Ω z , which does not even take into account the fact that the wave temperature depends more strongly on Σ −1 than does the bulk temperature. These conclusions are consistent with the results shown in Figure 12.
Dependence of TC on Σ −1 To understand the dependence of the Comptonization temperature TC on Σ −1 , we first look at the dependence of TC,g and TC,w, individually. If the gas and wave temperature profiles both decreased proportionally to the same power of Σ −1 , then both TC,g and TC,w would also decrease in proportion to this power of Σ −1 because the size of the bulk Comptonization region would remain unchanged. But since the wave temperature profile decreases faster than the gas temperature profile, the effect on TC,g and TC,w also depends on other factors, such as the slopes of the gas and wave temperature profiles in the region. Since v 2 turb ∼ Σ −2 gives a fairly strong dependence on density, we expect TC,w to uniformly decrease with decreasing Σ −1 . But since v 2 th,ph ∼ Σ −1/4 gives a very weak dependence on density, it is hard to see whether TC,g will increase or decrease with decreasing Σ −1 . Either way, we expect the dependence of TC,g on Σ −1 to be weaker. Figure 14 confirms these expectations. Finally, since TC = TC,g +TC,w, and since by definition the main contribution to TC is from TC,w, we expect TC to strongly increase with Σ −1 . Figure  12 confirms this expectation.

Dependence of bulk Comptonization on the Reynolds stress fraction
So far we have assumed that the β parameter, defined as the ratio of the vertically integrated Reynolds stress to the vertically integrated total stress, is held constant. We note that this is not to be confused with the plasma β, which is the ratio of the plasma pressure to the magnetic pressure. In radiation MHD simulations it is typically found that β ∼ 0.2. We now show how varying β affects bulk Comptonization. The turbulent velocity scaling, equation (19), becomes (K17) while the thermal velocity scalings, equations (20)-(21), remain unchanged. Since the dependence of v th,ph on Σ −1 is weak, we expect the dependence of bulk Comptonization on β to be similar to its dependence on Σ −2 . In Figure 15 we plot the dependence of the bulk Comptonization parameters on Ω −1 z for β/β0 = 1 and 4. As expected, we see the resulting curves are similar to those in Figure 12 corresponding to (Σ/Σ fid ) −1 = 1 and 2, respectively. We note that unlike the scaling for Σ −2 , the scaling for β is restricted to a much narrower range since β 1. In the rest of this work we suppress the dependence on β for clarity.

Dependence of bulk Comptonization on accretion disc parameters
Now that we have analyzed in detail the dependence of bulk Comptonization on the shearing box parameters Ω −1 z and Σ −1 , we proceed by relating these to the underlying accretion disc parameters. First we write Σ −1 in terms of the local flux F (K17): This says that for fixed Ω −1 z , Σ −1 is simply proportional to the local flux. Next, we need the scalings for Ω −1 z and F as functions of the accretion disc parameters. For Newtonian discs, they are 2 and where η is the efficiency assuming a no torque inner boundary condition, ∆ǫ is the change in efficiency due to a non-zero torque inner boundary condition (Agol & Krolik 2000), rin is the inner radius of the disc, anḋ 2 K17 is missing the η+∆ǫ η 0 +∆ǫ 0 factor in both the Newtonian and Kerr scalings for Σ.
For Kerr discs they are (Riffert & Herold 1995, K17) Ωz and where B, C, and D are functions of r and the spin parameter a, and go to unity for r ≫ rin. In order that the scalings for both Newtonian and Kerr discs be functions of the same underlying parameters, for Newtonian discs we set rin equal to the innermost stable circular orbit, which is in turn a function of the black hole spin parameter a. We note that sinceṁ = L/L Edd , one should not think ofṁ as the mass accretion rate. For example, for fixed mass M and fixedṁ, if we vary η + ∆ǫ (by varying the spin or the inner boundary condition) then the luminosity is unchanged since L =ṁL Edd and L Edd is proportional only to M (L Edd = 4πGM mpc/σT, so all the other parameters are constants of physics). But the mass accretion rate is NOT unchanged sinceṀ = L/ (η + ∆ǫ) c 2 . In other words, when varying other parameters (except for the mass) at fixedṁ, one should think of this as varying the mass accretion rate at fixed luminosity.
For the Kerr scalings, and 3 We now examine the dependence of the Comptonization parameters TC, τC, and yp,C on the accretion disc parameters.
Dependence on mass The dependence of the bulk Comptonization parameters on mass is straightforward. Since Ω −1 z is directly proportional to mass and Σ −1 is independent of mass, the dependence of the bulk Comptonization parameters on mass is identical to their dependence on Ω −1 z . That is, TC decreases weakly, τC increases weakly, and yp,C increases weakly with increasing mass. Furthermore, we can immediately regard the Ω −1 z axis in Figures 12 and 14 as the mass axis.
Dependence on luminosity The dependence of the bulk Comptonization parameters on L/L Edd is straightforward. Since Σ −1 is directly proportional to L/L Edd and Ω −1 z is independent of L/L Edd for both Newtonian and Kerr discs, the dependence of the bulk Comptonization parameters on L/L Edd is identical to their dependence on Σ −1 . That is, TC, τC, and yp,C all increase strongly with increasing L/L Edd . Furthermore, we can immediately regard the curves corresponding to different values of Σ −1 in Figures 12 and 14 as corresponding to different values of L/L Edd . In Section 5 we reproduce the plots from Figure 12 with the independent variables relabeled in order to summarize the dependence of bulk Comptonization on mass and luminosity.
Dependence on radius Both Ω −1 z and Σ −1 depend on r. But since all bulk Comptonization parameters depend strongly on Σ −1 and weakly on Ω −1 z , their dependence on r is almost entirely explained by the dependence of Σ −1 on r. For r ≫ rin, Σ −1 and hence the bulk Comptonization parameters increase with decreasing r for both Newtonian and Kerr discs. For ∆ǫ ≪ 1, Σ −1 eventually begins to decrease as r approaches rin, after which the bulk Comptonization parameters begin to decrease with decreasing r. The precise value of r below which the bulk Comptonization parameters begin to decrease differs slightly from the value of r at which Σ −1 begins to decrease because the bulk Comptonization parameters also depend on r through Ω −1 z , albeit weakly. For example, in Figure 16 we plot the dependence of Σ −1 and bulk Comptonization on r for ∆ǫ = 0. We see that the dependence of bulk Comptonization on r is well predicted by the variation in Σ −1 .
If, on the other hand, ∆ǫ is large enough, then both Σ −1 and the bulk Comptonization parameters monotonically increase with decreasing r, just as they do for r ≫ rin. This holds true for both Newtonian and Kerr discs. For example, in Figure 16 we also plot the dependence of Σ −1 and the Comptonization parameters on r for ∆ǫ = 0.05. We see that both Σ −1 and the bulk Comptonization parameters uniformly increase with decreasing r.
Dependence on spin For Newtonian discs Ω −1 z is independent of the spin parameter a, and for Kerr discs Ω −1 z depends on a only for r very close to rin through the functions C and B. But since all bulk Comptonization parameters depend strongly on Σ −1 and weakly on Ω −1 z , their dependence on a is almost entirely explained by the dependence of Σ −1 on a. For r ≫ rin, the dependence of Σ −1 on a is given by where for both Newtonian and Kerr discs η is a monotonically increasing (albeit different) function of a. We see that Σ −1 decreases with increasing a. The reason for this is straightforward. We recall that for fixed Ω −1 z , Σ −1 is proportional to the flux. As the spin and efficiency increase, the flux increases in the inner radii so for fixed luminosity L/L Edd the flux must decrease at large radii. The bulk Comptonization parameters therefore decrease at large radii with increasing spin. For example, in Figure 17 we plot the dependence of flux and bulk Comptonization on spin for r = 20, 12, and 7. We see that for r = 20, both flux and bulk Comptonization decrease with increasing spin.
For r sufficiently close to rin, on the other hand, since the flux increases with spin, so do the bulk Comptonization parameters. For example, in Figure 17 we see that for r = 7 flux and bulk Comptonization increase with spin until a ≈ 1. This is expected because as a approaches 1, rin approaches 1 and so r = 7 is no longer close to rin.
For an intermediate value of r (at which flux does not monotonically increase or decrease with spin), the dependence of the bulk Comptonization parameters on spin can still be understood by simply plotting the flux as a function of spin. For example, in Figure 17 we see that for r = 12 the dependence of bulk Comptonization on spin tracks the variation in flux.
Dependence on inner boundary condition The dependence of bulk Comptonization on the inner boundary condition is very similar to the dependence on spin. Since Σ −1 is proportional to the flux for fixed Ω −1 z , the dependence of bulk Comptonization on the inner boundary condition follows the variation in the flux. The inner boundary is parameterized in terms of ∆ǫ, the change in efficiency due to a non-zero inner torque. Since increasing ∆ǫ increases the flux in the inner radii, bulk Comptonization increases with increasing ∆ǫ in this region. At large radii, increasing ∆ǫ at fixed luminosity decreases the flux so that bulk Comptonization also decreases. For example, in Figure 18 we plot the dependence of flux and bulk Comptonization on ∆ǫ at large (r = 20) and small (r = 7) radii. In both cases we see that bulk Comptonization follows the variation in the flux.
Dependence on α For fixed Σ −1 and Ω −1 z , the variation of bulk Comptonization with α, given by equations (12)-(14), is reflected in Figure 12. For α = α fid , these plots are uncomplicated. Multiplying α by a constant k > 1 translates each curve for τC to the right on a log scale. For TC and yp,C, multiplying α by k not only translates each curve to the right but also multiplies each curve by 1/k. Since we plot yp,C on a log scale, for this variable multiplying α by k is equivalent to moving each curve to the right and downward.
Alternatively, one can think of multiplying α by a constant k > 1 as moving leftward along each curve for τC. For a given value of TC and yp,C, multiplying α by k is equivalent to not only moving leftward but also dividing the resultant value by k. To develop physical intuition, we plot the depen- dence of bulk Comptonization on Ω −1 z for multiple values of α in Figure 19. We see that bulk Comptonization overall decreases moderately with increasing α.
But in an accretion disc, we see that Σ −1 itself is directly proportional to α via the flux. The effect of varying α as an accretion disc parameter, then, affects bulk Comptonization primarily by varying Σ −1 . Since bulk Comptonization increases strongly with Σ −1 , increasing α generally increases bulk Comptonization. For example, in Figure  20 we plot the dependence of the flux and bulk Comptonization on α. We see that bulk Comptonization increases with increasing α for fixed accretion disc parameters, unlike for fixed shearing box parameters.

Dependence of bulk Comptonization on accretion disc parameters at the radius of maximum luminosity
To estimate the magnitude of bulk Comptonization for an entire accretion disc, we calculate the bulk Comptonization parameters at the radius rmax where the local luminosity is maximized. The luminosity at r is given by and rmax is determined by maximizing this function with respect to r. Since rmax is a function only of the spin a and the inner boundary condition parameter ∆ǫ, the dependence of bulk Comptonization on mass, luminosity, and α at this radius is the same as for fixed r, described in the previous section. The dependence on a and ∆ǫ, however, is different.
Dependence on spin We attempt to determine the dependence of bulk Comptonization on spin by analyzing how the flux depends on spin, as we did earlier. But in this case we have to be careful. Since rmax depends on a, r is not held constant. Given that we see that what really matters is the dependence of r 3/2 maxF on a. Fortunately, since the dependence of r 3/2 F on r is qualitatively similar to that of F , this does not change our physical intuition. As the spin parameter a increases, rmax decreases, and as the region of maximum luminosity becomes smaller, we expect that the flux at each point in this region must increase in order for the total luminosity L/L Edd to remain the same. In Figure 21 we plot the dependence of r 3/2 F and bulk Comptonization on a at rmax. We see that r 3/2 F and the bulk Comptonization parameters increase with spin, in agreement with our expectations. We note that for ∆ǫ > 0, it may be the case that rmax = rin rather than a value of r at which dL/dr = 0. Since rin decreases with a, however, this does not effect our conclusions. For example, in Figure 22 we plot the dependence of rmax on a for multiple values of ∆ǫ. We see that for sufficiently large ∆ǫ, rmax tracks rin until a is large enough that the value of r at which dL/dr equals zero is greater than rin. For very large ∆ǫ, rmax tracks rin for almost all values of a.
Dependence on inner boundary condition The dependence on the inner boundary condition parameter ∆ǫ is similar to the dependence on spin. As ∆ǫ increases, rmax decreases, and as the region of maximum luminosity becomes smaller, we expect that the flux at each point in this region must increase in order for the total luminosity L/L Edd to remain the same. Even once ∆ǫ is sufficiently large that rmax = rin, we expect the flux at rin to continue to increase since the increase in efficiency parameterized by ∆ǫ should result in an increase in flux at all radii near rin. In Figure 23 we plot the dependence of r 3/2 F and bulk Comptonization on ∆ǫ at rmax. We see that r 3/2 F and the bulk Comptonization parameters increase with ∆ǫ, in agreement with our expectations. For ∆ǫ slightly greater than zero, flux increases both because the flux at all inner radii increases with ∆ǫ and because rmax itself decreases towards smaller radii where the flux is larger. For larger values of ∆ǫ, rmax is fixed to rin, so the flux increases only due to the first effect and therefore increases at a slower rate. In Figure 24 we plot the dependence of rmax on ∆ǫ for multiple values of a. We see that for each curve rmax decreases until rmax = rin.

Effect of vertical radiation advection on bulk Comptonization
The scheme from K17 that we use to scale simulation data assumes that the flux is carried by radiation diffusion. Since we also see substantial vertical radiation advection in some radiation MHD simulations (Blaes et al. 2011), we attempt here to incorporate this process into our analysis of bulk Comptonization. Vertical radiation advection has both a di- rect and indirect impact on bulk Comptonization. The direct effect is to transport photons through the bulk Comptonization region faster so that they scatter fewer times, reducing bulk Comptonization. But the direct effect is typically negligible since vertical advection is significant only deep inside the photosphere (Blaes et al. 2011), outside the bulk Comptonization region. This is physically intuitive since radiation advection assists radiation diffusion in transporting photons out of the disc in order to maintain thermal equilibrium. Vertical advection is therefore most significant deep inside the photosphere where the photon diffusion time is comparatively large. The indirect effect of vertical radiation advection on bulk Comptonization, on the other hand, is to modify the underlying vertical structure gas and wave temperature profiles, which in turn either increases or decreases bulk Comptonization depending on whether shearing box parameters or the accretion disc parameters are held constant. In order to study this effect we need to incorporate vertical advection into the scaling scheme. One way to do this is to rederive the shearing box scalings without assuming that the flux is carried by radiation diffusion. We take this approach in Appendix B. Although this gives physical insight and is necessary to implement a specific model of advection, it is unnecessarily complex for our purpose here. Instead, we begin by simply adding an advection term Fa to the radiation diffusion equation (K17), which gives We observe that the only effect of adding Fa at fixed surface density Σ is to increase the total flux. Since this is also the equation that introduces the opacity parameter κ into the scaling scheme, it follows that the effect of adding Fa is the same as decreasing κ, as far as our scaling scheme is concerned. Conveniently, we see that the scheme in K17 already allows for scaling with respect to κ, even though for clarity we have suppressed the dependence on this parameter until now. We therefore proceed to determine the effect of advection on bulk Comptonization by varying κ.

Effect of advection for fixed shearing box parameters
To determine the effect of including vertical radiation advection on bulk Comptonization for fixed shearing box parameters, we need the turbulent and thermal velocity scalings with κ included (K17): We see that decreasing κ primarily affects the turbulent velocity magnitude. In particular, it moves the bulk and wave temperature profiles upward in Figure 13, increasing the size of the bulk Comptonization region. As advection increases at fixed Σ, therefore, we expect τC, TC, and yp,C to increase. For example, in Figure 25 we plot the dependence of bulk Comptonization on Ω −1 z for multiple values of κ. We see that this result is consistent with our expectations.

Effect of advection for fixed accretion disc parameters
To determine the effect of including vertical radiation advection on bulk Comptonization for fixed accretion disc parameters, we write Σ −1 in terms of the local flux F , this time allowing for the variation in κ (K17): Since neither Ω −1 z nor F depends on κ, combining this with equations (43)-(45) results in the following dependence on κ: We see that decreasing κ primarily affects the turbulent velocity magnitude, but in the opposite direction to the one in the previous section where the shearing box parameters are fixed. It moves the bulk and wave temperature profiles downward in Figure 13, decreasing the size of the bulk Comptonization region. As advection increases, therefore, we expect τC, TC, and yp,C to decrease. For example, in Figure  26 we plot the dependence of bulk Comptonization on Ω −1 z for multiple values of κ. We see that this result is consistent with our expectations.

Time variability of bulk Comptonization
We now explore the effect of the time variability of the vertical structure on bulk Comptonization. We stress that the specific numerical results of this section should not be directly compared to observations of real discs for two primary reasons. First, variability is an inherently global phenomenon which shearing box simulations therefore cannot effectively capture. Second, shearing box simulations with narrow box widths have been found to overestimate variability at a particular radius in the disc, so even not taking into account global phenomena we would expect this analysis to overestimate the variability of bulk Comptonization. The purpose of this section, therefore, is only to demonstrate how the time variability of bulk Comptonization depends on the time variability of the vertical structure profiles. To model the latter will require global disc simulations. We plot the standard deviation of the bulk Comptonization parameters over the 21 equally spaced timesteps in Figure 27. We also plot the standard deviation of the Comptonization gas and wave temperatures individually in Figure 28. In order to understand these results, we plot the standard deviations of the time-averaged temperature and density profiles in Figures 29 and 30, respectively.
The time variation of the density profile alone may cause significant variations in the bulk Comptonization optical depth. For example, in Figure 31 at each timestep we plot the Comptonization optical depth τC and the optical depth for a region of fixed size. The bottom of this region is taken to be the point at which the time-averaged gas and wave temperature profiles intersect in Figure 29. Variations in the optical depth of this region result in variations in τC that are caused by changes in the density profile alone, rather than changes in the size of the region itself. We see in Figure 31 that the overall variance of the optical depth of the region of fixed size is similar to the variance in τC. But we also see that the two quantities are only weakly correlated, which means that the variation in the bulk Comptonization optical depth must be due to other factors as well. For example, increasing the density may indirectly decrease the bulk Comptonization optical depth by reducing the wave temperature (see Section 2.2.2) and therefore the size of the bulk Comptonization region. Because the spatial variation of the gas temperature profile is so small in the bulk Comptonization region, we expect that its time variation should correlate with the time variation of the bulk Comptonization parameters in a predictable way. We expect that increasing the gas temperature decreases the size of the bulk Comptonization region, thereby increasing the Comptonization wave temperature and decreasing the Comptonization optical depth. In Figure  Figure 27. Dependence of bulk Comptonization on shearing box parameters for (Σ/Σ fid ) −1 = 1 (blue) and 5 (green). The shaded region corresponds to points within 0.675σ (i.e. 50% of the data for a Gaussian distribution) in the time variability. 32 at each timestep we plot the Comptonization gas temperature and the Comptonization wave temperature. We see that the two temperatures are strongly correlated in the direction we expect. In Figure 33 we plot the Comptonization gas temperature and optical depth. In this case, the correlation is also in the direction we expect, but it is weaker since density variations (among other factors) also play a significant role in determining the Comptonization optical depth.
We can also estimate the variability of the luminosity powered by bulk Comptonization according to the shearing box simulation data. The fraction of the luminosity powered by bulk Comptonization is just the total fractional photon energy change, which for unsaturated spectra is approximately ∆ǫ ǫ ≈ e y p,C − 1.
For yp,C > 1 we must check that spectra is unsaturated. For yp,C ≪ 1 spectra is always unsaturated, and in addition the fractional energy change simplifies to The variability of the luminosity powered by bulk Comptonization is therefore characterized by the fractional rms (root mean square) yp,C, which is the standard deviation Figure 28. Dependence of the Comptonization gas and wave temperatures on shearing box parameters, separately, for (Σ/Σ fid ) −1 = 1 (blue) and 5 (green). The shaded region corresponds to points within 0.675σ (i.e. 50% of the data for a Gaussian distribution) in the time variability. Figure 29. Horizontally and time averaged gas (blue), bulk (green), and wave (red) temperature profiles for the fiducial shearing box parameters (Table 8). The shaded region corresponds to points within 0.675σ (i.e. 50% of the data for a Gaussian distribution) in the time variability.  (Table 8). The shaded region corresponds to points within 0.675σ (i.e. 50% of the data for a Gaussian distribution) in the time variability. Figure 31. Time variability of the Comptonization optical depth τ C (blue) and the (normalized) optical depth of a nearly identical region whose physical size is defined to be constant (green). The surface density here is (Σ/Σ fid ) −1 = 2 and the other parameters are the fiducial shearing box parameters (Table 8). The timesteps are spaced 10 orbital periods apart, and each orbital period is 5535 seconds.
divided by the mean of yp,C, fractional rms = σy p,C yp,C .
We plot the fractional rms for (Σ/Σ fid ) −1 = 1 in Figure 34, and see that it is consistent with Figure 27. In particular, as the values on the x axis increase from 10 0 to 4 × 10 2 , we see that σy p,C decreases while yp,C increases so that the fractional rms increases. Since Ω −1 z ∝ M (Section 3.3), the fractional rms seems to vary insubstantially with mass. We note that the fiducial shearing box parameters (Table 8) correspond to r = 20 for the M = 2 × 10 6 M⊙, L/L Edd = 2.5 parameter set (Table 1). Since r ≈ 20 in the region of the disc that contributes most to the luminosity for this parameter set (i.e. for a = 0, ∆ǫ = 0), Figure 34 also characterizes Figure 32. Time variability of the Comptonization gas (blue) and wave (green) temperatures, normalized to the average Comptonization gas temperature. The surface density here is (Σ/Σ fid ) −1 = 2 and the other parameters are the fiducial shearing box parameters (Table 8). The timesteps are spaced 10 orbital periods apart, and each orbital period is 5535 seconds. Figure 33. Time variability of the Comptonization gas temperature T g,C (blue) and inverse optical depth τ −1 C (green), normalized to the average Comptonization gas temperature. The surface density here is (Σ/Σ fid ) −1 = 2 and the other parameters are the fiducial shearing box parameters (Table 8). The timesteps are spaced 10 orbital periods apart, and each orbital period is 5535 seconds.
the variability of the luminosity powered by bulk Comptonization for the entire accretion disc.
We note that the time-averaged bulk Comptonizaton parameters are not equal to the bulk Comptonization parameters computed from the time-averaged temperature profiles. For example, in Figure 35 we plot the bulk Comptonization parameters computed with the time-averaged data. To do this, we first time average the gas and wave temperature profiles and then compute the bulk Comptonization parameters. In the same figure we also plot the timeaveraged parameters, originally plotted in Figure 12. We see that the parameters computed from the time-averaged profiles significantly overestimate the time-averaged param- eters. This result is important because it means that the time-averaged profiles, while often useful, should not directly be used to model bulk Comptonization.
This work is based on only 21 simulation snapshots spaced 10 orbital periods apart, but we note that by applying our model to a complete set of simulation data one could also calculate how the power spectra of the vertical structure profiles affect that of the fraction of the luminosity powered by bulk Comptonization and other bulk Comptonization parameters.

Comparison of results with previous work and observations
Aside from the fact that our model implements a simplified version of the procedure used in K17 to calculate the bulk Comptonization parameters, our work differs from K17 in two important ways: In K17 bulk Comptonization is modelled for an entire accretion disc at once rather than at each radius individually, and shear velocities are included in addition to turbulent velocities. Because of these differences, the bulk Comptonization parameters found in that work depart slightly from those found here. But the dependence of bulk Comptonization on accretion disc parameters detailed in Sections 3.3 and 3.4 is consistent with the results of K17.
In particular, the bulk Comptonization y parameter for the overall disc increases with α and mass, while the bulk Comptonization temperature decreases with increasing mass. And within a given disc, bulk Comptonization is greatest at intermediate radii where the flux is also near maximal.
To make contact with observations, K17 modelled bulk Comptonization for a few systems with accretion disc parameters similar to those fit by D12 to REJ1034+396, a narrow-line Seyfert 1 (NLS1) with L/L Edd = 2.4. The bulk Comptonization parameters found by K17 broadly agree with those fit by D12. This agreement suggests that the large soft X-ray excess seen in REJ1034+396 may at least in part be due to bulk Comptonization. By generalizing the results of K17, our work provides a physical basis for more widely connecting warm Comptonization models of the soft excess to underlying accretion disc parameters.

The importance of the disc inner boundary condition and implications for black hole X-ray binaries
An important consequence of our results is that bulk Comptonization is strongly dependent on the disc inner boundary condition parameter, ∆ǫ. Before proceeding, however, we provide context for the range of ∆ǫ since it is not a widely used parameter and we need to have a sense of what it means for it to be large. To start, we observe that since the efficiency for a no torque inner boundary condition, zero spin system is η = 0.057, any value of ∆ǫ > 0.01 is relatively large. Even for spin a = 0.9, the efficiency with no inner torque is η = 0.16, so ∆ǫ = 0.1 corresponds to a substantial physical change. Another way to understand the effect of ∆ǫ is to examine how it affects the disc scalings presented in K17. We see that this parameter arises in the equations for the flux scalings, reproduced in Section 3.3, equations (31) and (34). To understand why ∆ǫ appears here, we observe that for ∆ǫ = 0 in both equations the purpose of the final factor is to ensure that the flux goes to zero as r approaches rin rather than continue to increase as r −3 . It follows that we can regard ∆ǫ large to the extent that it reverses the effects of this factor. For example, for the Newtonian scalings we see that setting ∆ǫ = 1/rin removes the dependence on r of this term altogether so that F ∼ r −3 . For zero spin, rin = 6 so we should regard ∆ǫ = 0.17 as very large. For a = 0.9, rin = 2.32, so the critical value of ∆ǫ is 0.43. Therefore, values of ∆ǫ anywhere from 0.1 to 0.4 should be viewed as very large, depending on the spin parameter a. The Kerr scalings lead to similar conclusions.
As ∆ǫ approaches infinity the flux scaling asymptotes to a fixed value rather than continuing to increase. At r = rin, we see from the Newtonian scalings that this limit is reached when ∆ǫ ∼ 1. Beyond this point, therefore, bulk Comptonization hardly varies at all with ∆ǫ. The reason for this is that ∆ǫ changes the distribution of flux throughout the disc at a fixed overall luminosity L/L Edd . For ∆ǫ ≫ 1, the distribution at the inner radii is fixed and the flux distribution continues to change only for r ≫ rin.
In Figure 36 we plot the dependence of bulk Comptonization on mass for several values of ∆ǫ for a system with moderate (a = 0.5) spin. We see that all bulk Comptonization parameters strongly increase with increasing ∆ǫ. Note that this dependence holds only for the region where the disc is brightest, not radii for which r ≫ rin (see Sections 3.3 and 3.4).
An important implication of this result is that bulk Comptonization is likely insignificant in black hole X-ray binaries unless the luminosity greatly exceeds Eddington or ∆ǫ is large (i.e. ∆ǫ > 0.1). In Figure 36, for which L/L Edd = 2.5, for example, we see that for ∆ǫ = 0 bulk Comptonization is non-existent for M < 10 5 M⊙. For M = 10M⊙, we see a significant effect only for ∆ǫ > 0.1.

Robustness of bulk Comptonization results to variations in the disc vertical structure
In Section 3.2.2 we showed that the dependence of bulk Comptonization on Σ −1 and Ω −1 z , plotted in Figure 12, can be understood in terms of the shearing box temperature profiles, plotted in Figure 13. Since these profiles correspond to scaled data from a single radiation MHD simulation, we need to examine the extent to which our results are robust to changes in the disc vertical structure that may occur in shearing box simulations run in different regimes or global simulations. Certainly the exact values of the bulk Comptonization parameters are sensitive to such changes (Section 3.6), but we now show that the overall dependence on Σ −1 and Ω −1 z (and therefore on the accretion disc parameters) is more robust.
We first consider the dependence of bulk Comptonization on Ω −1 z . In Section 3.2.2, using the profiles shown in Figure 13, we showed that since only the gas temperature profile varies with Ω −1 z , the Comptonization temperature decreases and the Comptonization optical depth increases with increasing Ω −1 z . There is considerable uncertainty in the shape of the gas temperature profile outside the scattering photosphere, but fortunately the contribution of this region to bulk Comptonization is negligible since the bulk Comptonization temperature is weighted by the optical depth fac- tor τ dτ . The greatest uncertainty in this analysis, therefore, is the bulk velocity field, which determines the shape of the wave temperature profile. But since the wave temperature is defined to strongly decrease with increasing density (Section 2.2.2), we expect that even for significantly different velocity fields the wave temperature profile will increase near the scattering photosphere and that the resulting dependence of the Comptonization parameters on Ω −1 z will be unchanged. Since our conclusions in Section 3.2.2 regarding the dependence of bulk Comptonization on Σ −1 also rely primarily on the fact that the wave temperature profile strongly increases near the photosphere, we also expect them to be robust to changes in the vertical structure.
In addition to the above concerns, we must also check that as the size of the bulk Comptonization region increases it remains outside the effective photosphere. Otherwise, only part of the bulk Comptonization region will contribute to bulk Comptonization (since photons are emitted at the effective photosphere). However, this condition is likely always satisfied since the size of the bulk Comptonization region increases most significantly as Σ −1 increases, which simultaneously moves the effective photosphere inward. In particular, for our data we find that the vertical structure becomes ef-fectively thin well before the bulk Comptonization region optical depth is more than a small fraction of the optical depth of the half-thickness of the disc.

Limitations to the scaling scheme parameter range
We now make note of a subtlety that limits the applicability of the scaling scheme from K17: The scheme can scale data to lower surface densities, but not higher ones. To understand why, we examine how the gas temperature profile scales with decreasing Σ. First we note that since the gas temperature profile below the photosphere is significantly different from the profile above it, the regions must be scaled separately and then joined together. Next, we observe that the photosphere is not defined to be at a set number of scale heights h away from the midplane but rather at the point at which the scattering optical depth is unity. As a result, as the surface density Σ decreases the photosphere moves inward in z/h. Since h is the fundamental length scale for variations in the vertical structure and since fewer scale heights of data are needed outside the photosphere, fewer grid cells of data are needed to fill the region outside the photosphere of the scaled disc. This truncated data is scaled appropriately and then joined to the scaled data from above the photosphere. We see, therefore, that scaling to smaller surface densities requires deleting grid cells from the original simulation data. By the same reasoning, this scheme cannot scale to larger surface densities since it would require data from more grid cells than already exist. It immediately follows that this scheme cannot scale data to any set of accretion disc parameters for which Σ/Σ0 > 1. In particular, since Σ −1 is always directly proportional to the luminosity L/L Edd , we can never scale to smaller values of L/L Edd unless they are offset by simultaneously scaling to, for example, smaller radii or greater ∆ǫ. We note that this scaling scheme is, therefore, useful for scaling lower Eddington ratio simulations to higher ones, as we do in this work, but not the other way around.
We showed in Section 3.2.2 that bulk Comptonization increases strongly with increasing Σ −1 . For the curve in Figure 12 with the smallest value of Σ −1 , Σ −1 = Σ −1 fid = 4Σ −1 0 , we see that 0.1 < yp,C < 0.3. Therefore, the fact that we cannot scale to values of Σ −1 smaller than Σ −1 0 is not a significant limitation since it appears that bulk Comptonization is negligible for such values anyway. But this analysis assumes that the scaling scheme in K17 is valid over an arbitrarily large parameter range. If we want to scale to a regime with significantly different opacities, for example, then we really should use data from simulations with the relevant opacities included. For example, if the vertical structure is significantly different for sub-Eddington AGN because of changes in the opacities that occur in such regimes, then bulk Comptonization could be larger than we would infer from our analysis of the 110304a simulation data. On the other hand, this seems unlikely given that absorption opacities will substantially increase in this regime.

Effect of bulk Comptonization on disc spectra
The effect of bulk Comptonization on disc spectra cannot only be to upscatter photons to higher energies because we also must take into account the back reaction on the disc vertical structure. Since energy conservation fixes the flux as a function of radius and the other accretion disc parameters, we expect that bulk Comptonization will be accompanied by a decrease in the gas temperature at the effective photosphere so that the total emitted flux will remain unchanged. For significant bulk Comptonization, the effect of this is to move the Wien tail to higher energy while moving the spectral peak to lower energy, broadening the spectrum. For moderate bulk Comptonization, the effect of lowering the gas temperature may not translate into a leftward shift of the spectral peak, but the spectrum will still be broadened in such a way that the total flux remains unchanged.
A decrease in the effective photosphere gas temperature is the simplest conceivable back reaction. This would occur if the only effect of bulk Comptonization on the gas is to remove kinetic energy from the turbulent cascade through radiation viscous dissipation (Kaufman & Blaes 2016) so that less kinetic energy is dissipated and converted to gas internal energy. But to self-consistently model this phenomenon, bulk Comptonization must be implemented in the underlying radiation MHD simulations. The shearing box simulations used in this paper ), for example, do not include bulk Comptonization since it is primarily a second order effect in velocity (see K17 Section 4.2), and the flux limited diffusion approximation does not capture second order effects (Kaufman & Blaes 2016).

Effect of the horizontally averaged z component of the velocity field on bulk Comptonization
Bulk Comptonization includes effects that are both first and second order in the velocity field (K16) and the wave temperature defined in Section 2.2.2 captures only the second order effects. The first order effect is non-zero only for compressible modes and is negligible when the photon mean free path is large relative to the mode wavelength. As long as the mean photon energy is less than 4kB (Tg + Tw), the second order effect always results in upscattering, analogous to thermal Comptonization. But the first order effect can result in either upscattering or downscattering depending on whether the velocity field is converging or diverging, respectively (K16). It follows that only long wavelength compressible modes should result in a non-negligible first order effect, since for shorter wavelength modes either the first order effect is negligible or upscattering in one region is offset by downscattering in another. Therefore, the variations with respect to z of the density weighted, horizontal average of the z component of the velocity field may result in a non-negligible first order effect. As in the case of the wave temperature profile (Section 2.2), density weighting is appropriate because photons scatter more times in higher density regions. We expect such long wavelength variations to exist since the vertical structure is stratified. For example, in Figure 37 we plot this profile at the 140 orbits timestep for the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table  1) with r = 14. In the remainder of this section we show that this effect is discernable but subdominant to the second order effect. We also show that once this effect is taken into account the slight discrepancy in Figure 3 between the spectra calculated directly with the turbulence and the spec- tra calculated by modeling the turbulence with the wave temperature vanishes. We therefore conclude that the wave temperature models the second order effect more accurately than we originally had reason to believe based on the preliminary analysis in Section 2.2.2. All data in this section are scaled to the M = 2 × 10 6 M⊙, L/L Edd = 5 parameter set (Table 1). As in Section 2.2, all spectra and vertical structure profiles correspond to the 140 orbits timestep.
To begin, we calculate spectra with the original velocity field, both with and without subtracting off the horizontally averaged z component, and plot the results for r = 14 in Figure 38. We see that the spectrum computed with the horizontally averaged z component included is shifted to slightly lower energies. The spectra at the other radii illustrate the same effect. Since any additional second order effect associated with this component can only increase upscattering, this energy shift must either be due to the first order effect or vertical radiation advection. As explained in Section 3.5, vertical radiation advection transports photons through the bulk Comptonization region faster, which decreases the number of photon scatterings in the region and may therefore reduce the overall second order effect. In order to show that the energy shift is predominantly due to the first order effect, not radiation advection, we calculate spectra with uniform temperature profiles both for the case of no velocities and for the case where only the horizontally averaged z component is included, and plot the results for r = 14 in Figure 39. The spectra at other radii illustrate the same effect. Since a uniform temperature profile with no velocity field has no effect on the base spectrum there is no second order effect, and so adding in the horizontally averaged z component of the velocity field can shift the resulting spectra to lower energies only through the first order effect, not radiation advection. Since the spectrum in Figure 39 is shifted by the same amount as in Figure 38, we conclude that the original shift is predominantly due to the first order effect, not vertical radiation advection. As a check on this analysis, we repeated the uniform temperature profile spectral calculations but included instead the negative of the  (Table 1) computed with (green) and without (blue) the velocities. For the red curve, the spectrum was computed with velocities but the horizontally averaged z component of the velocity field was subtracted from the total z component.  (Table 1), computed with data truncated at τs = 20. All gas temperatures were set to the horizontally averaged value at the base. For the blue curve, the velocities were not included, and for the green curve only the horizontally averaged z component of the velocity field was included.
horizontally averaged z component of the velocity field, and found that the energy shifts were opposite in direction and equal in magnitude. We also estimate the energy shift due to the first order effect heuristically and check that the result is consistent with our spectral calculations. The fractional energy change per scattering due to this effect is (K16) In this case, therefore, the region where this effect is greatest is near the photosphere (Figure 37), where we conveniently just confirmed that vertical advection is dominated by diffusion. It follows that the average number of scatterings dN in a region of optical depth dτ is approximately (Section 3.2.1) equal to 1.6(2τ dτ ). The total approximate fractional energy change f in this region is then Since the fractional energy change per scattering is much smaller than unity, In this case we find that at all radii f ≈ −0.1, consistent with the results in Figures 38 and 39. We note that if the total fractional energy change is also much less than unity, such as in this case, then In order for Monte Carlo calculations to selfconsistently capture the first order effect, one must take into account the time dependent nature of the problem (K16), either by performing time dependent simulations or by careful analysis of the results. This is because this effect can result in either upscattering or downscattering depending on whether a region is converging or diverging, and a diverging region will typically evolve into a converging one on the flow timescale, given by If the region is near the photosphere, such as in the case examined here, then the photons escape the region on the diffusion timescale, which is shorter, and the first order effect will on average broaden the spectrum. If the region is sufficiently deep inside the photosphere that the diffusion timescale, exceeds the flow timescale, then the flow will change significantly before photons can diffuse very far. This is the case for standing acoustic modes, for example (Blaes et al. 2011). In this case, results from time independent Monte Carlo simulations can be trusted only if the spectrum is negligibly affected by the upscattering or downscattering in such regions. In order to capture only second order effects in a Monte Carlo simulation, we can first subtract off the horizontally averaged z component of the velocity field, but we can do this only if its second order effect is negligible. To investigate this, in Figure 40 for r = 14 we plot the original wave temperature profile along with the wave temperature profile computed by first subtracting off the horizontally averaged z component of the velocity field. We see that the resulting two curves are essentially identical except in a small region where bulk Comptonization is negligible since τs ≪ 1.  The horizontally averaged vz profile contributes negligibly to the wave temperature both because it contributes negligibly to the underlying bulk temperature profile, plotted in Figure 41, and because the wave temperature downweights long wavelength variations (Section 2.2.2). To confirm that the contribution of the horizontally averaged vz profile to the second order effect is negligible, we calculate spectra in which we model the turbulence with wave temperatures calculated both with and without including the horizontally averaged vz and plot the results for r = 14 in Figure 42. We see that there is no discrepancy between the respective curves, consistent with the wave temperature profiles in Figure 40. The spectra at the other radii illustrate the same effect.
In Section 2.2.2 we showed in Figure 3 that spectra computed with data in which the velocities were turned off and the wave temperatures were added to the gas temperatures approximated spectra computed with the velocities. Since the wave temperature captures only second order effects and since the horizontally averaged z component of the velocity field results in a non-negligible first order effect, we expect that when this component is subtracted off the approximation will improve. We perform this comparison in Figure 43 for r = 14. We see that in this case the approximation is so good that the respective spectra are indistinguishable from each other. The spectra at the other radii illustrate the same effect. This is not only consistent with our prediction but shows that the wave temperature captures second order effects even more accurately than we originally had reason to believe based on the preliminary analysis in Section 2.2.2. Given these results, we also expect that spectra computed with the velocities turned off except for the horizontally averaged z component and with the wave temperatures added to the gas temperatures will coincide with spectra computed with the velocities turned on, since both sets of spectra should capture both first and second order effects. We plot these spectra in Figure 44 for  (Table 1) for the gas temperature Tg (blue), bulk temperature T bulk (green), and bulk temperature computed by first subtracting off the horizontally averaged z component of the velocity field (red). The dashed lines denote where τs = 1 and τs = 10.  (Table 1) computed without the velocities. For the green and red curves, the wave temperatures were added to the gas temperatures. For the red curve, the wave temperatures were computed by first subtracting off the horizontally averaged z component of the velocity field. r = 14 and see that they agree with our prediction. The spectra at the other radii illustrate the same effect.

SUMMARY
We have simplified the bulk Comptonization model of K17 in order to explore a larger space of accretion disc parameters and develop greater physical insight into this phenomenon. Rather than fit the temperature and optical depth to spectra computed with Monte Carlo post processing simulations, we developed a procedure to calculate the Comptonization temperature and optical depth directly from the underlying vertical structure data (Section 2). Using this, we plotted  (Table 1) computed with (green) and without (blue) the velocities. For the spectrum computed with velocities, the horizontally averaged z component of the velocity field was subtracted from the total z component. For the red curve, the velocities were not included but the wave temperatures were added to the gas temperatures.  (Table 1) computed with (green) and without (blue) the velocities. For the red curve, the z component of the velocity field was horizontally averaged, the x and y components were set to zero, and the wave temperatures were added to the gas temperatures.
the dependence of the Comptonization parameters on the shearing box parameters and showed how these results can be understood in terms of the one dimensional temperature profiles (Sections 3.1 and 3.2). We then showed how we can analytically determine the dependence of bulk Comptonization on each accretion disc parameter individually (Sections 3.3 and 3.4). Our principal results are as follows. The primary independent variables in a shearing box are the surface density Σ, the vertical epicyclic frequency Ωz, and the strain rate, ∂xvy. We also allow α, the ratio of the vertically integrated stress to the vertically integrated total pressure, to vary. For Kerr discs the scalings for the strain rate and vertical epicyclic frequency are always nearly equal (equation 11), which leaves three independent parameters. Using the velocity scalings (equations 19, 20, and 21), we showed that the dependence of the Comptonization parameters on α can be subsumed into the other parameters (equations 12, 13, and 14), which reduces the parameter space to two variables, Σ and Ωz.
We plotted the dependence of the bulk Comptonization temperature, optical depth, and y parameter on Σ and Ωz ( Figure 12). We showed that these results can be understood from analyzing the one dimensional temperature profiles ( Figure 13) and the velocity scalings (equations 19, 20, and 21). In particular, the Comptonization optical depth and y parameter increase strongly with increasing Σ −1 and weakly with increasing Ω −1 z . The Comptonization temperature also increases strongly with increasing Σ −1 , but decreases weakly with increasing Ω −1 z . To determine the dependence of bulk Comptonization on accretion disc parameters, we write Σ in terms of F and then write the scalings for F and Ωz in terms of mass, luminosity, radius, spin, and inner boundary condition (Section 3.3). Since Ω −1 z is directly proportional to mass, and Σ is independent of mass, the dependence of bulk Comptonization on mass is identical to its dependence on Ω −1 z . Similarly, since Σ −1 is directly proportional to luminosity, and Ω −1 z is independent of luminosity, the dependence of bulk Comptonization on luminosity is identical to its dependence on Σ −1 . Therefore, Figure 12 also summarizes the dependence of bulk Comptonization on mass and luminosity. Here, for clarity, we reproduce the plots from Figure 12 in Figure 45 with the independent variables labeled as mass and luminosity.
The dependence of bulk Comptonization on the other accretion disc parameters is inferred by analyzing how they affect Σ −1 since bulk Comptonization depends much more strongly on Σ −1 than it does on Ω −1 z . Since Σ −1 is proportional to the flux F (equation 29), we showed that the dependence of bulk Comptonization on the other disc parameters can be understood intuitively in terms of how they effect F . In particular, at large radius (i.e. r ≫ rin) bulk Comptonization always decreases with increasing radius. At small radius, whether bulk Comptonization increases or decreases with radius depends on the inner boundary condition parameter ∆ǫ. This parameter is defined to be the increase in the efficiency η for a non-zero stress inner boundary condition, i.e. L = (η + ∆ǫ)Ṁ c 2 . Using the same line of reasoning, we showed that bulk Comptonization increases with both spin and the inner boundary condition parameter ∆ǫ at small radius (r ≈ rin), and decreases with those parameters at large radius. Finally, we showed that bulk Comptonization increases with α, since once the accretion disc parameters are substituted in for Σ and Ωz, Σ −1 itself becomes proportional to α (equation 29) and this outweighs the dependence on α discussed earlier.
Next we studied bulk Comptonization for an entire accretion disc by examining how it varies when the radius is fixed to the region of maximum luminosity (Section 3.4). The dependence of bulk Comptonization on mass, luminosity, and α is unchanged from above since the radius of maximum luminosity does not vary with these parameters. But since this radius does depend on the spin and inner boundary condition parameter ∆ǫ, the dependence of bulk Comp- . Dependence of bulk Comptonization on mass. The blue, green, red, and cyan curves correspond to L/L fid = 1, 2, 3.3, and 5, respectively, where L fid /L Edd = 2.5. The parameters held constant are r = 20, a = 0, ∆ǫ = 0, and α/α 0 = 2 (i.e. the other parameters in the M = 2 × 10 6 M ⊙ , L/L Edd = 2.5 parameter set (Table 1), with r = 20). The only difference between this figure and Figure 12 is the labeling of the axes.
tonization on these parameters required a new treatment. We showed that in this case bulk Comptonization always increases with spin and ∆ǫ. In Section 3.5 we showed that the effect of including vertical radiation advection at a fixed radius in an accretion disc is to decrease bulk Comptonizaton. We discussed how to include advection in our model more formally in Appendix B.
In Section 4.1 we pointed out that our results broadly agree with the results of K17, which in turn agree with the analysis by D12 of the narrow-line Seyfert 1 REJ1034+396.
An important result of this work is that bulk Comptonization is strongly dependent on the disc inner boundary condition (Section 4.2). In particular, the larger that ∆ǫ is, the lower the luminosity can be without bulk Comptonization being negligible. Figure 46 summarizes the dependence of bulk Comptonization on mass and luminosity for ∆ǫ = 0.2. In both Figures 45 and 46 the radius is fixed to the region where the luminosity is greatest, but for ∆ǫ = 0.2 this radius corresponds to the innermost stable circular orbit. By comparing these two figures we see that for a given luminosity bulk Comptonization is significantly greater for ∆ǫ = 0.2. We also showed ( Figure 36) that bulk Comptonization is negligible in black hole X-ray binaries unless the disc inner boundary condition parameter is very large (∆ǫ ∼ 0.1) or the luminosity greatly exceeds Eddington. We expect that in a real disc bulk Comptonization at a given radius will be accompanied by a decrease in the gas temperature at the effective photosphere in order to leave the flux unchanged, which is required by energy conservation (Section 4.5).
Because our model connects the bulk Comptonization parameters to the disc vertical structure one dimensional temperature profiles in a way that is physically intuitive, it provides a useful framework for understanding bulk Comptonization even in situations in which some of our specific results may not hold, such as shearing box or global radiation MHD simulations run in entirely different regimes.
Since our results outline how bulk Comptonization depends on fundamental accretion disc parameters, an observer who fits the soft X-ray excess with a warm Comptonization model can use them to distinguish contributions to the soft X-ray excess due to bulk Comptonization from those due to other physical mechanisms. In high Edding-ton sources, this can help provide a physical basis for and therefore constrain warm Comptonization models of the soft excess. Bulk Comptonization is likely insignificant in lower Eddington flows, on the other hand, even though the data show that in these flows the soft excess carries a more significant fraction of the power Mehdipour et al. 2011Mehdipour et al. , 2015. Moreover, since bulk Comptonization depends on the properties of MRI turbulence through α (Sections 3.3 and 3.4) and through the time variability of the temperature and density profiles (Section 3.6), our work indicates that observations of the soft X-ray excess may in turn advance our understanding of disc turbulence in the radiation pressure dominated regime.
Here we derive the shearing box scalings presented in K17 without assuming that the flux is carried by radiation diffusion. We give scalings for ρ, Tg, v turb , and vs in terms of Σ, Ωz, ∂xvy, α, β, f col , and h. The result, therefore, of not assuming radiation diffusion is to leave the scale height h as a free parameter.
To begin, the scaling for the density profile is still given by K17 equation (13), except that h is now a free parameter: The flux scaling is determined by K17 equations (2) The scalings for the shear velocity, pressure and gas temperature profiles are unchanged, except that h is now a free parameter. The scaling for the shear velocity profile is given by K17 equation (22): vs (x) = ∂xvy ∂xvy,0 h h0 vs,0 (h0x/h) .
The pressure profile is given by K17 equation (17), The scaling for the gas temperature profile is given by K17 equations (19), (20), and (21): T 4 g,in (z) =T 4 g,ph + Pc Pc,0 T 4 g,0 (h0z/h) − T 4 g,0 (h0z ph /h) T 4 g,ph = P ph,in P ph,in,0 T 4 g,ph,0 T 4 g,out (z) = P ph,in P ph,in,0 T 4 g,0 (z ph,0 + h0(z − z ph )/h) . (B10) B2 Modelling radiation advection with an effective κ The scaling for h depends on how the radiation is vertically transported. For radiation diffusion, for example, the scaling for h is given by K17 equation (7): In Section 3.5 we pointed out that including radiation advection is equivalent to simply decreasing κ as far as the shearing box scalings are concerned. Now we also see that including radiation advection (at fixed surface density Σ) is therefore equivalent to simply increasing the scale height relative to the value set by radiation diffusion alone. We note that although including advection increases the scale height in a shearing box, it has the opposite effect at a fixed radius in an accretion disc. This is because in a shearing box the surface density Σ is fixed and the total flux F is allowed to vary. In this case, including advection increases both the scale height and therefore the flux (equation B2). But if we substitute in equation (B2) everywhere for Σ, we can instead regard F as a free parameter instead of Σ. Since F is a function of accretion disc parameters such as the mass, mass accretion rate, radius, etc., this procedure gives the shearing box scalings in terms of accretion disc parameters rather than shearing box parameters. The scalings that result from this procedure, assuming radiation diffusion, are given in K17 Appendix C. The scaling for the scale height, for example, is given by We see, therefore, that including advection decreases the scale height at a fixed radius in the disc. In this process, the flux is held constant and the surface density, therefore, increases (equation B2).

APPENDIX C: ADDITIONAL FIGURES
In this section we show the plots of spectra at multiple radii omitted from Section 2.  Figure C1. Normalized spectra (red) computed by passing a 50eV Planck source (blue) through vertical structure data truncated at τs = 10 at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1). In all cases the velocities are zeroed and the wave temperatures are added to the gas temperatures. The green curves are calculated by using the Kompaneets equation to pass the 50eV Planck source through a homogeneous medium with temperature T 1D , given in Table 2. The spectra for only r = 14 were originally plotted in Figure 7.  Figure C2. Normalized spectra at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 5 parameter set (Table 1) computed with (red) and without (blue) velocities. The green curves are calculated by using the Kompaneets equation to pass the blue curves through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 3. The spectra for only r = 14 were originally plotted in Figure 8.  Figure C3. Normalized spectra at multiple radii for the M = 2 × 10 8 M ⊙ , L/L Edd = 5 parameter set (Table 1) computed with (red) and without (blue) velocities. The green curves are calculated by using the Kompaneets equation to pass the blue curves through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 4. The spectra for only r = 14 were originally plotted in Figure 9.  Figure C4. Normalized spectra at multiple radii for the M = 2 × 10 6 M ⊙ , L/L Edd = 2.5 parameter set (Table 1) computed with (red) and without (blue) velocities. The green curves are calculated by using the Kompaneets equation to pass the blue curves through a homogeneous Comptonizing medium with parameters T C and τ C , given in Table 5. The spectra for only r = 14 were originally plotted in Figure 10.