{'fraction_non_alphanumeric': 0.033999385182908085, 'fraction_numerical': 0.01922328107388052, 'mean_word_length': 4.531172069825437, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local Simulations of Heating Torques on a Luminous Body in an Accretion Disk

A luminous body embedded in an accretion disk can generate asymmetric density perturbations that lead to a net torque and thus orbital migration of the body. Linear theory has shown that this heating torque gives rise to a migration term linear in the body's mass that can oppose or even reverse that arising from the sum of gravitational Lindblad and co-orbital torques. We use high-resolution local simulations in an unstratified disk to assess the accuracy and domain of applicability of the linear theory. We find agreement between analytic and simulation results to better than 10\% in the appropriate regime (low luminosity, low thermal conductivity), but measure deviations in the non-linear (high luminosity) regime and in the high thermal conductivity regime. In the non-linear regime, linear theory overpredicts the acceleration due to the heating torque, which we find to be due to the neglect of non-linear terms in the heat flux. In the high thermal conductivity regime linear theory underpredicts the acceleration, although here both non-linear and computational constraints play a role. We discuss the impact of the heating torque for the evolution of low-mass planets in protoplanetary disks, and for massive stars or accreting compact objects embedded in AGN disks. For the latter case, we show that the thermal torque is likely to be the dominant physical effect at disk radii where the optical depth drops below a critical value.


INTRODUCTION
Planets, stars or compact objects orbiting within accretion disks perturb surrounding gas due to gravitational forces (Goldreich & Tremaine 1980), accretion (Bondi 1952), the release of heat or radiation (Masset 2017), and winds (Gruzinov et al. 2020). It is commonly the case that the resulting leading and trailing density perturbations are asymmetric, producing a gravitational back-reaction and a non-zero torque on the body. The torque leads to orbital migration-either an increase or a decrease in the semi-major axis-and evolution (usually in the sense of damping) of any eccentricity or inclina-tion. In many circumstances of interest the time scale for migration is short compared to the disk lifetime, making it probable that observable properties of the system are substantially shaped by the effects of migration.
The longest-studied torque is that due to the purely gravitational perturbation of the disk-embedded object. It is made up of two independent components, one from waves excited at Lindblad resonances and one exerted in the co-orbital region (Kley & Nelson 2012), both of which scale as the square of the object's mass. The net Lindblad torque (summing the opposite-signed contributions from interior and exterior resonances) has some dependence on disk properties, but is mostly due to intrinsic asymmetries in the interaction and almost always leads to inward migration (Ward 1997). The co-orbital torque, on the other hand, can lead to either inward or outward migration, and depends in a complex way on numerous properties of the disk (including radial gradients of vortensity and entropy, viscosity, thermal diffusivity, and disk winds; Paardekooper et al. 2011;Mc-Nally et al. 2020).
Numerical simulations show that thermal effects, either in the disk gas in the vicinity of the planet or associated with the release of heat or radiation from a luminous body, result in additional torques (Lega et al. 2014;Benítez-Llambay et al. 2015;Chrenko et al. 2017). Unlike the purely gravitational torques, thermal effects can (in principle) remain significant even for very low mass bodies. In particular, Masset (2017), using linear perturbation theory, identified a "heating torque" that arises when an orbiting body injects thermal energy into the surrounding disk. The thermal energy leads to the formation of low-density lobes near the planet, which are generically asymmetric, producing a torque. The heating torque can be on the same order of magnitude as other torques that cause migration (such as the Lindblad torque), and typically leads to outward migration.
The consequences of thermal torques on the migration rate of disk-embedded objects have been studied in the context of low-mass planet formation, where Lindblad torques alone would cause planets with masses of the order of the Earth's mass to migrate toward the central star on a timescale shorter than the disk lifetime. The luminosity on these mass scales typically results from pebble accretion (Ormel & Klahr 2010;Lambrechts & Johansen 2012). The heating torque modifies the predicted map of where in the disk inward and outward migration occur (Guilera et al. 2019), though the consequences for the final population of planets that form may be modest (Baumann & Bitsch 2020).
Heating torques could also impact the migration rate of luminous bodies such as stars and accreting stellarmass black holes, which can be captured (Syer et al. 1991) or form (Shlosman & Begelman 1987;Goodman 2003;Levin 2007;Dittmann & Miller 2020) in the gas disks around supermassive black holes. Heating torques could interact with other gas torques (e.g. Lindblad torque, corotation torques) to form a migration trap-a radius in Active Galactic Nuclei (AGN) disk where the net torque is zero. Such migration traps would host an increased density of objects and provide a possible formation location for intermediate-mass black holes (Bellovary et al. 2016) or for stellar-mass black hole binaries (Secunda et al. 2019;Tagawa et al. 2019). Stellar-mass binaries merging within an AGN disk could contribute to the observed LIGO population (Stone et al. 2017;Abbott et al. 2019), while stellar-mass black holes merging with the central supermassive black hole are future LISA sources, whose detailed properties may be modified by migration torques (Derdzinski et al. 2019).
Heating torques have been studied analytically (Masset 2017) and using global numerical simulations (Lega et al. 2014;Benítez-Llambay et al. 2015;Chrenko et al. 2017). Here, we complement these prior studies using a local shearing box model for the disk. By simulating a luminous body in the limit where its mass goes to zero, using 32 zones per characteristic wavelength of the heating torque, we are able (a) to isolate the heating torque from the cold thermal torque and (b) to fully resolve the influence of the heating torque on the disk. Our work effectively extends the thorough numerical investigation of a luminous body travelling through a homogeneous medium (Velasco Romero & Masset 2019 to the case of a luminous body embedded within a shear flow. The main questions we seek to answer are: 1. What are the numerical prerequisites needed to reproduce the (Masset 2017) linear theory, and how accurate is that theory when the approximations involved are relaxed? 2. When do non-linear effects set in, and how do they change the linear theory's prediction for the thermal torque? 3. Is the heating torque important for stars and accreting compact objects embedded within AGN disks? The structure of the paper is as follows: we summarize the analytic theory for the heating torque resulting from a luminous body in a shear flow in §2.1 and describe our numerical methods in §2.2. Our numerical results are presented in §3. §4.1 discusses the limits to the analytic theory, and §4.2 discusses applications of the model to luminous objects in AGN and protoplanetary disks. We conclude in §5.

Analytic Results
Analysis of the local hydrodynamic equations with thermal conductivity shows that the effects of a massive body's gravitational potential and its luminosity on a surrounding disk can be separated and studied independently in the linear regime (Masset 2017, Eq. 34). We take advantage of this separation to focus solely on the "heating torque", the torque due to the density perturbation that is sourced by thermal energy diffusing outward from a luminous body. Fig. 1 illustrates how the asymmetry in this perturbation, resulting from the displacement of the orbiting body from co-rotation, leads to a net torque. To aid in the interpretation of our numerical results, we summarize the key assumptions and results from Masset (2017).  Figure 1. Illustration of the physics leading to a heating torque. Heat diffusing away from a luminous, disk-embedded body, is sheared by the Keplerian disk flow, forming hot low density lobes. These lobes are asymmetric interior to / exterior to the body, because the body is displaced from exact co-rotation due to the presence of a pressure gradient in the disk gas. The gravitational back-reaction from the heated lobes exerts a torque on the orbiting body.
Masset (2017) linearizes the hydrodynamic equations (see eqs. (9) to (11)), assuming a steady-state, in a local ("shearing box") frame that co-rotates with the orbiting body. In this local frame, x corresponds to the radial direction, y to the azimuthal, and z to the vertical direction (perpendicular to the disk midplane) as illustrated in Fig. 1. If there is a radial pressure gradient in the disk, there is an offset x p between the orbiting body and disk gas that has the same orbital velocity. This distance from co-rotation is given by, Here p 0 (ρ 0 ) is the equilibrium background pressure (density), q is the shearing parameter (equal to 3/2 for Keplerian disks), and Ω 0 is the angular velocity of the local frame. For typical pressure profiles that decrease as a function of radius x p is positive, implying that the body will sit further away from the central body than the gas rotating at the same angular velocity, experiencing a headwind. Three characteristic scales enter the problem: the distance from co-rotation x p , the characteristic size of the density perturbation caused by the body's luminosity λ c , and the pressure scale height of the disk H. In the linear calculation it is assumed that the following hierarchy holds, We refer to the first requirement for scale separation (Eq. 2) as Assumption II and the second (Eq. 3) as As-sumption III (the first assumption is that of linearity). Assumption III allows the vertical density gradient of the local box to be neglected (justifying our use of unstratified simulations), while the small parameter associated with Assumption II is used extensively to expand the expected gravitational force from the under-density caused by the body's luminosity. The relative importance of these two hierarchies and the validity of the predicted net force on the body (Eq. 5) is explored in § 4.1.1. The characteristic size of the disturbance and the net force experienced by the body as a result of the heating torque are predicted to be (Masset 2017, Eq. 83, and Eq. 109), Here γ is the adiabatic index (γ = 5/3 in all of the following simulations), c 2 s = γp 0 /ρ 0 is the equilibrium sound speed, L is the luminosity emitted by the body, χ is the disk's thermal conductivity, and M is the mass of the body. Crucially, Eq. 5 is linear in the mass of the body. This feature of the heating torque allows us to calculate the force per unit mass (i.e. the body's acceleration) without needing to explicitly include the body's mass at all in the simulations.
The heating torque is of interest because it can be the same order of magnitude as the other torques in the system (such as the Lindblad torque). Defining the heating torque can be written as (Masset 2017, Eq. 144) where is of the order of the Lindblad torque. Here r 0 is the semi-major axis of the body, M * is the mass of the central object, and H is the pressure scale height of the disk. The formula given by Eq. 5 for the net force on a body (due to the asymmetric gravitational forces caused by the body's luminosity distributed by differential rotation) is predicted to hold under three conditions: • Assumption I: perturbations of density and pressure should be much less than equilibrium values (linearity, ρ ρ 0 ).
• Assumption II: the offset from corotation x p should be much less than the size of the disturbance λ c ( x p λ c ).
• Assumption III: the disturbance should be much smaller than the pressure scale height of the disk (λ c H).
In this work, we test the validity of the linear theory when one or more of these assumptions is violated.

Simulations
The linear theory is developed in the local "shearing sheet" approximation (Masset 2017), which translates directly into a well-studied numerical set-up. We solve the inviscid hydrodynamic equations in a local approximation of a Cartesian box rotating around a massive body (a star or black hole, for instance) with orbital frequency Ω 0 , and add a source term to the energy density equation to model the luminosity. With ρ the mass density, e the energy density, P the pressure, and V the velocity, the equations read, ∂ρV ∂t where F H = −χρ∇ (e/ρ) is the heat flux and L is the total luminosity emitted by the body. Φ t = −qΩ 2 0 (x−x p ) 2 is the tidal potential due to the central object. The vertical density gradient is neglected, both for consistency with Masset (2017) and for the same physical reasons discussed there, and the radial density gradient is modelled through a non-zero offset from co-rotation (it is neglected in the shearing-box, as justified by assuming that the background pressure does not change significantly over the short radial scales under consideration). The shearing parameter q is equal to 3/2 for the Keplerian flows studied in this work.
The Athena++ code is used to solve the above equations in the luminous body's rest frame on a uniform Cartesian mesh (White et al. 2016;Felker & Stone 2018, Stone et al., 2020 with the HLLC solver. The simulation's origin sits at the radial location where the gas orbits at the same frequency as the luminous body, such that the position of the body is fixed over the course of the simulation (10 orbits). We discuss the consequences of neglecting the radial motion of the body in response to the generated torque in §4.2.2. Both the origin and the body sit at the mid-plane of the disk (z = 0). In the fiducial run L1K1, the simulation domain spans [4.13, 12.4, 4.13]H in the x (radial), y (azimuthal), and z (vertical) directions respectively, where H is the pressure scale height of the disk (defined through H = c s /Ω 0 ). The fiducial run L1K1 has a resolution of [256, 192, 256] cells (i.e. [62.0, 20.7, 62.0] cells/H) and a value x p = 0.097 H. Convergence with resolution and domain size is studied in §3.5.
The body's luminosity is modelled by directly injecting internal energy density into the gas via the energy density equation (11). In the analytic theory the injection term is L δ(x − x p ) (as in Eq. 11). Since injection at a single point is not possible numerically, we implement this term in the simulations by adding an energy density v × ∆t at each time to each cell whose center lies within an injection radius r rad . Here, v is the (constant) luminosity per volume and ∆t is the time step as determined by the Courant condition. The total luminosity L injected at each time step can be calculated as where n is the number of cells included in the injection region and v is each cell's volume; thus the total luminosity L depends on both the luminosity per volume v and the injection radius r rad . Unless otherwise specified, the injection radius is set such that the luminosity is evenly distributed into the eight cells neighboring the body. The effect of a finite injection radius (rather than a strict Dirac delta function) is explored in §3.5. We note again that according to the linear theory the torque due to the gravitational poten-tial of the body can be separated from the torque due to the body's luminosity (Masset 2017). In this work we model only the heating torque, and do not include the gravitational potential of the body.
As is standard for the shearing-box set-up, all simulations use periodic boundary conditions in the azimuthal and vertical directions and shearing-periodic boundary conditions in the radial direction, the effects of which are discussed in §3.5. Generic units H (disk scale height) for length, Ω −1 (2π/Ω is one orbit) for time, and P 0 (background pressure) for energy density are used. The sound speed thus has units of HΩ, acceleration has units of HΩ 2 , thermal diffusivity has units of H 2 Ω, and total injected luminosity has units of P 0 H 3 Ω. These values can be scaled to various astrophysical systems, which are discussed in §4.

Diagnostics
Linear theory provides a prediction for the total gravitational force F y experienced by the orbiting body as a result of the perturbed gas density (purely from the body's luminosity, not its gravity). We calculate this force per unit body mass in the simulation on spherical shells by calculating the distance between the body and each cell, assuming that all of the cell's mass is located at its center, and using the inverse square law with GM = 1 in code-units. The resulting acceleration can be plotted as a function of radius or summed over radius to directly compare to Eq. 5. The noisiness of the data in Fig. 3 is due to the fact that no interpolation was used in the calculation of the acceleration. To avoid introducing artificial asymmetry to the force (i.e. a systematically larger force on the x < 0 side), the summation stops at the shortest radius that fits inside the simulation domain in all directions. The limiting radius is thus L x − x p , where L x is the half-width of the box in the x (radial) direction and x p is the body's distance from corotation. This restriction is not expected to change the results significantly since the excluded portions are not a large fraction of the box, the gas there is not as perturbed, and the force from the gas there is attenuated by the inverse square of the radius).
The fractional changef in a quantity f is useful to establish the linearity of density perturbations, where f (t) = f (t) − f (t = 0) is the perturbation from equilibrium. We describe a simulation as being in the linear regime if the deviation from equilibrium values is no more than 5%.
To facilitate direct comparison with previous work (specifically, Masset (2017)'s Fig. 1), we calculate the perturbation in surface density σ as the k z = 0 mode of the Fourier transformρ(x, y, k z ) of the density perturbation ρ in the z-direction, i.e.
Separating into the effect due to zero offset from corotation (σ (0) ) and the first order effect due to nonzero offset (σ (1) ), we find (Masset 2017, Eq. 114) σ (15) Note that σ is the perturbation of the surface density; the unperturbed surface density σ is constant. 2) lead to perturbations that are less than 5% of the equilibrium values, which we take to be in the linear regime. The measured pressure perturbations are two orders of magnitude smaller than the density perturbations. This is consistent with Masset (2017)'s estimation (Eq. 36) that P H 2 Ω 2 0 ρ ≈ c 2 s ρ . With the value c 2 s = γP/ρ = 1.00 H 2 Ω 2 , P should be much less than ρ . The non-linear regime is illustrated in the right column of Fig. 2, which injects two orders of magnitude more energy per timestep ( v = 142 P 0 Ω; simulation L10K1). The qualitative appearance of the pressure and density perturbations remain similar for this much higher rate of energy injection. Both simulations, as expected, quickly reach an equilibrium within approximately two orbits.

Net Acceleration as a Function of Radius
We take the v = 1.42 P 0 Ω run as our fiducial simulation L1K1 so as to be firmly in the linear regime. Using the technique described in §2.3, we plot as a function of radius the net gravitational force on the body per unit body mass as a result of the gas perturbed by the object's luminosity. Fig. 3 shows the result of summing up gas in front of and behind the body (y > 0 and y < 0, respectively), as well as the total force (green; right-hand scale). The bottom panel shows how the net acceleration differs ahead/behind the body from the initial acceleration (which is non-zero on either side but which sums to zero) as well as how close the net acceleration summed over all radii is to the linear prediction (horizontal lines). We obtain agreement with the linear prediction to better than 10%.
The heating torque was previously studied by Benítez-Llambay et al. (2015) using global simulations. Masset (2017) noted that the numerical torque found in that work was smaller than the analytic prediction by approximately an order of magnitude. Our results suggest that this discrepancy was due to resolution; using the higher resolution possible in a local setup we obtain a good level of agreement with the analytic theory. The remaining discrepancies are usually small, exceed the analytic prediction, and reduce with higher resolution. Exceptions are simulations with thermal conductivity much smaller than the fiducial value. This trend suggests that a more precise value depends on the exact physics close to the body, which would be better modelled with a self-consistent luminosity prescription to capture accretion onto the body.

Perturbation of the Surface Density
In addition to the value of the net force acting on the body, Masset (2017) calculates a map of the surface density perturbation predicted by the analytic theory. Fig. 4 reproduces this map with simulation data (compare to Masset (2017)'s Fig. 1). The upper panel Fourier transforms the density perturbation of a simulation with no offset (x p = 0), as outlined in §2.3. The lower panel Fourier transforms the density perturbation with the equilibrium value of x p , then subtracts the zero offset case and divides by x p to extract σ (1) . The perturbations are smaller in amplitude than expected, largely because the peak expected amplitude is very close to the luminous body and is not as resolved. However, the general shape of the perturbation agrees well with the linear prediction.

Scaling Relations
Linear theory predicts a linear dependence of the net gravitational force on the total luminosity L emitted by the body and a power-law dependence F y ∝ χ −3/2 on the thermal conductivity (Eq. 5). To test these predictions we ran two suites of simulations: one that fixes the thermal conductivity and varies the total emitted luminosity, and one that fixes the total emitted luminosity and varies the thermal conductivity.
For the first suite we fixed χ = 0.017 H 2 Ω and varied v over three orders of magnitude: from v = 0.142 P 0 Ω to 142 P 0 Ω. Simulations are considered to be in the linear regime if the perturbation never exceeds 5% of the equilibrium value. Fig. 5 reveals a tight agreement with the linear prediction even an order of magnitude into the non-linear regime (indicated by green squares). As the injected luminosity increases even more, the linear theory begins to over-predict the measured force because of non-linear effects; this over-prediction is discussed in §4.1.2. In the linear regime at least, we are able to reproduce both the scaling and the normalization of the net force to within 10%. Assessing the validity of the analytic prediction for the scaling of the net force with thermal conductivity is substantially harder, because changing the conductivity also changes the characteristic wavelength λ c . It is difficult to find a numerically tractable set of parameters that both (a) remains in the linear regime and (b) maintains the hierarchy of scales required by Masset (2017), over a substantial range in χ. Fig. 6 shows the measured dependence of the net force as a function of the thermal conductivity, at fixed luminosity. For sufficiently low values of the thermal conductivity, heat cannot diffuse away fast enough, causing the system to enter the non-linear regime (indicated by a green square). For high values of the thermal conductivity, the required scale separation x p λ c H is lost (shown as orange dots). Only the blue crosses, Fitting the data only at intermediate χ, we find that the dependence of net gravitational acceleration on conductivity is close to χ −1 , rather than the expected χ −3/2 . We caution, however, that this fit is made over only a very limited range of χ. If, instead, we consider simulation data at higher values of χ, we find a dependence that appears to be closer to the analytically predicted power-law. The ideal regime for matching the linear theory appears to be around χ = 0.017 H 2 Ω, which in the simulations presented has ratios x p /λ c = 0.187 and λ c /H = 0.519. Simulations using χ = 0.0061 H 2 Ω to obtain ratios x p /λ c = 0.311 = λ c /H measured an acceleration lower than the linear prediction by a fac- Table 1. Summary of simulations used to calculate the scaling of net acceleration with total injected luminosity L (Fig. 5). All simulations have χ = 0.017 H 2 Ω, corresponding to λc = 0.52 H, while all other parameters (e.g. resolution, offset of the body from corotation, described in the text) are that of the fiducial simulation L1K1. This set of simulations keeps λc/H and xp/λc constant at 0.52 and 0.187, respectively. * indicates a simulation that has density fluctuations greater than 5% of the equilibrium value and has thus entered the nonlinear regime. Simulation L1K1 is often referred to as the fiducial simulation, and L10K1 as the high luminosity simulation. Measured values are presented as the average between 1 and 10 orbits plus/minus one standard deviation of the value in time.
Figure 5. Net acceleration due to the gas's gravity as a function of total injected luminosity. Thermal conductivity is fixed at χ = 0.017 H 2 Ω. The linear theory's prediction (slope: 1.0; red dotted line) matches the data well even into the non-linear regime, although the fit (slope: 1.0 ± 6.5 × 10 −7 ; black dashed line) was determined using only the linear data points. Simulations are summarized in Table 1.
tor of two. This suggests that the requirement that x p λ c is more important for matching linear theory than λ c H. This hierarchy of the assumption is reasonable since the former is used in expanding the net force, whereas the latter is used to drop vertical density stratification (Masset 2017); see §4.1.1.

Numerical Considerations
In order to assess the robustness of the numerical results, we explored the dependence of the simulation results on domain size, resolution, boundary conditions, and injection radius. Of these factors, we find that the Figure 6. Net acceleration due to the gas's gravity as a function of thermal conductivity. Total luminosity is fixed with L = 1.96 × 10 −4 P0ΩH 3 ( v = 1.42 P0Ω). The linear theory (red dotted line) predicts a power-law index of −1.5, whereas the fit (black dashed line) determines a power-law index of −1.0±1.6×10 −4 and total luminosity L = 6.7×10 −3 P0ΩH 3 . The fit was determined solely from the simulations satisfying the hierarchies xp/λc < 0.3 and λc/H < 0.6 (blue crosses). Simulations are summarized in Table 2. most important numerical effects are related to the size of the injection region. The analytic assumption that all the body's luminosity is deposited at a single point is both an approximation to the physical situation, and an idealization that cannot be achieved in grid-based numerical simulations. We find that for the fiducial parameters, and a resolution that allows for an injection radius of r rad = 0.04 H, we can obtain agreement with the linear theory at the level of 7.6%. Doubling the injection radius to r rad = 0.07 H, at half the resolution, Table 2. Summary of simulations used to calculate the scaling of net acceleration with conductivity χ (Fig. 6). All simulations have L = 1.9 × 10 −4 P0ΩH 3 ( v = 1.42 P0Ω), while all other parameters (e.g. resolution, offset of the body from corotation) are that of the fiducial simulation L1K1. The value of λ/H is easily read off; the value of xp/λc is obtained by noting that xp = 0.097 H. * indicates a simulation that has density fluctuations greater than 5% of the equilibrium value and has thus entered the non-linear regime. Simulation L1K1 is often referred to as the fiducial simulation, and L1K10 as the high conductivity simulation. Measured values are presented as the average between 1 and 10 orbits plus/minus one standard deviation of the value in time.
To isolate the effect of changing spatial resolution from the effect of differing injection radii, we test for convergence with spatial resolution by keeping the same injection radius and changing the resolution. Due to the discretization of the region around the body, increasing the resolution will result in an injection region that closer approximates a sphere rather than a rectangular prism (as is the case for the low resolution simulation, which injects energy evenly into eight neighboring cells). Because of the slight change in injection volume, the total injected luminosity will also be modified; since we have an excellent prediction of what a simulation with a slightly different total luminosity would be (see Fig. 5), we can control for the difference in total luminosity and isolate the influence of the injection region's shape. We compare two simulations, both with an injection radius of 0.07 H, conductivity χ = 6.1 × 10 −3 H 2 Ω, and injected luminosity per volume v = 1.42 P 0 Ω but one with fiducial resolution and the other with half the fiducial resolution, resulting in total injected luminosity 2.5 × 10 −4 P 0 ΩH 3 and 1.5 × 10 −4 P 0 ΩH 3 , respectively. We find that the net force per unit mass agrees between these runs at approximately the 10% level. (Note that for this value of the conductivity neither the high nor the low resolution simulation recover the analytic prediction to high accuracy.) From Fig. 3 it is apparent that the heating torque arises from within approximately 0.5H of the body, well within the size of our fiducial simulation domain. Nonetheless, the use of periodic boundary conditions does introduce artefacts that are visible in the plots of the pressure perturbations as structures close to the edges of the box that re-appear on the opposite side of the box (sheared, in the case of the y-edges). To test for possible errors introduced by the use of periodic boundaries, we compared simulations in which the box size was increased to twice that of the fiducial simulation L1K1's in each direction, while maintaining all other variables constant. The measured accelerations agreed to better than 1%. Somewhat larger changes, at the 5% level, occur if we compare against a box with half the resolution, but three times the box size, of the fiducial simulation L1K1. We conclude that the low-level periodic pressure perturbations visible in the simulation do not impact the force calculation at the level of accuracy we are interested in here.
Finally, we note that the simulations assume that the luminous body's location remains fixed over a small multiple of the local dynamical timescale. In principle, for sufficiently high luminosities and local disk surface densities, the resulting torque might be able to migrate the body fast enough to invalidate this assumption. Analogous physics has been studied in the context of gravitational torques, where motion of the gravitating body can lead to a dynamical co-rotation torque and "Type III" migration (Masset & Papaloizou 2003;Paardekooper 2014). We do not explore this possibility here, but note that caution and additional study would be needed in any circumstance where the implied migration speed due to the heating torque exceeded a fraction of HΩ.

Limits of the Linear Theory
The analytic theory for the heating torque relies both on linearity, and on satisfying hierarchical separation between the scales of the displacement from co-rotation, the induced density perturbation, and the disk scale height. By numerically solving the full set of hydrodynamic equations, we can test the limits of these various assumptions.

Testing the Hierarchy Requirements
The first set of assumptions are the hierarchies given by Eq. 2 and Eq. 3, i.e. that x p λ c H. Fig. 6 shows how the derived acceleration scales as these assumptions are broken. The simulations closest to the analytic prediction do not have equal ratios of x p /λ c and λ c /H; rather, they prefer a smaller x p /λ c . As thermal conductivity increases, λ c becomes larger whereas the offset from corotation x p and disk scale height H stay constant. This results in a decrease in the ratio x p /λ c and an increase in λ c /H, i.e. Eq. 2 becoming better satisfied and Eq. 3 becoming less satisfied. The result is that the characteristic wavelength of the perturbations are less well-contained by the simulation domain, as illustrated by the high conductivity simulation L1K10's temperature perturbation (green circles) in Fig. 7, which does not decay to zero even at the edge of the simulation domain in either direction (vertical or radial). Extending the simulation domain to a scale height where all three scales are separated by a factor of 10 (i.e. the offset from co-rotation is a factor of 100 smaller than the scale height) would require upwards of 1200 cells in each direction to resolve a single scale height, which is not even large enough to capture the full decay of higher conductivity simulations (Fig. 7). The computational cost of such simulations is beyond the scope of this work. Figure 7. Time-averaged azimuthal profiles at the body's position (not the origin) of the density (dashed line) and temperature (solid line) perturbations δρ and δT for the fiducial simulation L1K1 (blue crosses), and a high conductivity simulation L1K10 (orange dots), and a high luminosity simulation L10K1 (green squares). Each line has been normalized to its maximum value and time-averaged over the last seven orbits, t = 3 to 7 orbits. Filled-in portions denote one standard deviation over time.

Non-linear Effects
Another assumption is that the perturbations are small compared to their equilibrium values: ρ ρ 0 . By increasing the luminosity, we can study how the acceleration departs from the linear prediction as we enter the non-linear regime. The scaling relation of acceleration with fixed conductivity (Fig. 5) shows that higher luminosity simulations that are in the non-linear regime (indicated by green squares) measure a smaller acceleration than the linear theory would predict. This difference could be due to either the computational issues with resolving hierarchies listed in the previous section or to non-linear effects that are not adequately captured by the linear theory. In this section we show that the measurement is physical, i.e. that the linear theory's neglect of higher-order terms leads to an overprediction of the actual non-linear net acceleration.
There are two properties of the steady-state perturbations that could contribute to the final density distribution: their shape and their amplitude. Plots of the density perturbation (Fig. 7) hint at the similarity between the high luminosity (L10K1) and the fiducial simulation (L1K1) for the linear ρ distribution. Fig. 8 shows the non-linear term ∇ · [χρ ∇T ] (the divergence of the heat flux's contributions by non-linear terms; this term is subtracted from the internal energy in Eq. 11). From this figure it is clear that the non-linear term has essentially the same shape for the fiducial (L1K1) and the high luminosity run (L10K1), and a different shape for the high conductivity run (L1K10). The magnitude of the heat flux is of course drastically different; the in-flow of heat to the origin is 1.7×10 4 times higher for the high luminosity run than for the fiducial simulation because heat does not diffuse away very fast (compare this to the high conductivity run L1K10, which has a lower maximum because the heat is more spread out). This similarity shows that the relevant length scales are indeed the same for these two simulations; the differences in the resulting density distribution must therefore be due to the magnitude of the heat flux rather than its shape. In contrast, the high conductivity run L1K10 (which has a larger characteristic wavelength λ c ∝ χ 1/2 ) is fundamentally different in its shape; the differences in the resulting density distribution are due to at least partially to the non-linear term, but is also influenced by the computational difficulties outlined in §4.1.1 of resolving the hierarchies.
The difference in density distribution is evident in Fig. 7, most notably that the asymmetry is not the same. To quantify what this means for the net acceleration of the body, Fig. 9 plots the net acceleration from the gas as a function of distance from the body (very similar to Fig. 3's plot of the sum of the acceleration in either direction), normalized to the values of the fiducial run at every point. Fig. 9 shows that the high luminosity (L10K1) run's acceleration profile is modified only close to the body, where it is lower than might be expected from linear theory. On the other hand, the high conductivity run L1K10 has a long tail of density perturbations that contribute significantly to the net acceleration even at large distances, hence why the linear theory underpredicts the net acceleration.
We briefly note that the separation of the effect of the body's luminosity from its gravitational potential is only valid in the linear regime. In the non-linear regime, interaction between these two effects (which linearly act in the same way to provide a net outward migration) could result in deviation from the linear prediction. Exploring aspects of this interaction is left to future studies.

Stars or Accreting Compact Objects in a Thin Disk
Migration processes may be important in geometrically thin accretion disks in Active Galactic Nuclei (AGN). Stars may form within such disks as a consequence of gravitational instability (Goodman 2003;Levin 2007), and they may also be captured from a cluster whose orbits intersect the disk gas (Syer et al. 1991). Either circumstance could lead to a population of luminous stars, or accreting stellar-mass compact objects, orbiting within an AGN disk.
The full ramifications of having a population of stellarmass objects within AGN disks are complex, and we do Each line has been normalized to its maximum value (0.027 and 1.7 × 10 4 times the fiducial simulation's maximum value, respectively) and time-averaged over t = 3 to 7 orbits. Error bars denote one standard deviation over time. Figure 9. Snapshot of the radial profile of the net acceleration on the body at a time of 9 orbits for a high conductivity simulation L1K10 ( v = 1.42 P0Ω, χ = 0.085 H 2 Ω; orange dash-dot line), and a high luminosity simulation L10K1 ( v = 142 P0Ω, χ = 0.017 H 2 Ω; green dotted line), normalized at every point to the fiducial simulation L1K1. As in Fig. 3, the injection radius r rad = 0.04 H for all three simulations and half the characteristic wavelength λc/2 = 0.26 H for the fiducial and high luminosity simulation are plotted. The characteristic wavelength of the high conductivity run L1K10 extends further to the right, at a distance of λc/2 = 1.8 H from the body.
not discuss them here. Rather, we assume that we have a single luminous object orbiting on a circular, noninclined orbit, with the same sense of rotation as the disk gas. The question we seek to answer is whether the heating torque is large enough, compared to previously studied torques arising at the Lindblad and co-orbital resonances, that it should be included in models of migration within AGN disks.
We assume, consistent with our numerical results, that the analytic result given as Eq. 7 provides a good estimate of the ratio of the heating torque, Γ heat , to the fiducial torque scaling, Γ 0 , In this expression x p is the radial displacement of the object from co-rotation, L is the luminosity of the object, and λ c and L c are characteristic scales for the thermal perturbation and for the luminosity, respectively. We assume a Keplerian disk (q = 3/2) in which the pressure p 0 ∝ r −n . At r = r 0 the relevant quantities can then be written as, Here χ is the thermal diffusivity in the disk gas surrounding the luminous object. If the diffusivity is physically the result of radiative diffusion, we can write (Paardekooper et al. 2011), where κ is the opacity, T is the temperature, and σ is the Stefan-Boltzmann constant. For both massive stars and accreting compact objects, the Eddington luminosity provides a very rough but reasonable estimate of how the likely luminosity scales with the mass. We write, simulations and those presented in this paper. In particular, in Lega et al. (2014) the planet has mass but no luminosity (and hence is affected by only the cold thermal torque), while in Benítez-Llambay et al. (2015) the planet is both luminous and massive (and hence is affected by both the cold thermal and the heating torque). Relevant parameters for all three simulations are summarized in Table 3.
Although the results of both Lega et al. (2014) and Benítez-Llambay et al. (2015) highlighted the importance of thermal effects for an accurate assessment of the migration rate, there was a significant mismatch between the quantitative values obtained and the subsequent analytic theory of Masset (2017). Inspection of Table 3 suggests that this discrepancy may well be due to the difficulties inherent in resolving the relevant scales in a global simulation. Compared to previous simulations (Lega et al. 2014;Benítez-Llambay et al. 2015), our fiducial simulation L1K1 better resolves the characteristic wavelength of the density perturbation. Lega et al. (2014) has approximately two cells spanning λ c (Masset 2017), whereas Benítez-Llambay et al. (2015) has approximately 4.7 cells to resolve λ c . By using a local domain, our fiducial simulation L1K1, which resolves λ c with 41 cells, is better poised to capture the full effect of the thermal torques. Limited resolution could be one of the reasons that Benítez-Llambay et al. (2015) sees a net force about an order of magnitude below the predicted value (Masset 2017). As for physical parameters, our fiducial simulation L1K1's luminosity is three orders of magnitude larger than those presented in Benítez-Llambay et al. (2015), leading to a large value for the heating torque. However, the results should continue to scale down to lower values of luminosity, where the assumption that the planet does not change its distance from corotation over the course of the simulations should be more accurate.
Our simulations also explore a different regime in terms of the scale hierarchy. The simulation of Lega et al. (2014) has H/λ c = 14, whereas our simulations have H/λ c = 1.5; similarly, in Lega et al. (2014), λ c /x p = 2.3 and in our fiducial simulation L1K1 λ c /x p = 6.88 (Table 3). We have better scale separation between λ c and x p at the expense of less scale separation between λ c and H. Our closer-to-linear results support our argument that the first criterion is more essential to the linear theory than the second. As discussed in §4.1.1, the restriction on scale height is not a physical requirement but rather an ease-of-computational one; therefore this work simply explores a slightly different physical parameter regime.
We conclude that the mismatch between previous simulations' measurement of the heating force and the linear theory's prediction is likely due to poor resolution. In local shearing box simulations, the relevant small scales are better resolved and the net acceleration is within 10% of the linear theory in the linear regime.

CONCLUSIONS
In this work, we have used three-dimensional local simulations to test and extend the analytic theory for the heating torque developed by Masset (2017). The heating torque arises from the interaction between a luminous disk-embedded body and Keplerian shear, which distorts low-density regions that were heated by the body into asymmetric lobes that exert a net gravitational force on the planet. In the regime where the resulting density perturbations are linear, we find good agreement between the results of our direct numerical simulations and the analytic theory. We surmise that prior global simulations probably lacked enough resolution of the energy injection region, leading to an underestimate of the magnitude of the heating torque. Going beyond the linear theory, we explored regimes of high thermal conductivity and high luminosity. We find that at high luminosity the derived torque is smaller than the linear prediction, and attribute this as being due to non-linear terms in the heat flux. In the high conductivity regime we infer a higher acceleration than predicted by the linear theory. We find that both the non-linear terms in the heat flux and computational limitations contribute to this larger value.
At the linear level the heating torque can be considered separately from other contributions to the migration of disk-embedded bodies. Although numerically convenient, there are few if any physical circumstances where gravitational (Kley & Nelson 2012) and other thermal torques (Lega et al. 2014) would not also need to be considered. In most cases, study of these torques requires a combination of analytic theory, local, and global numerical simulations, whose results can partially be encapsulated in relatively simple torque formulae (Paardekooper et al. 2011;Jiménez & Masset 2017). In the context of low-mass planet migration, using such formulae, the heating torque is estimated to be most important (relative to the sum of all other torques) for masses ∼ 0.1M ⊕ (Baumann & Bitsch 2020).
A second environment where the heating torque might be important is for the migration of luminous objects (massive stars or accreting compact objects with luminosities of the order of the Eddington limit) in AGN disks. Using simple scaling arguments, and the analytic theory of Masset (2017), we find that the heating Table 3. Summary of important values for different simulations. Lega et al. (2014) studies the cold thermal torque (no luminosity), while Benítez-Llambay et al. (2015) studies the cold and heating torques in a semi-global simulation. The values in the table were obtained either from Sec. 5.3.2 in Masset (2017) ( †), by private communication ( * ; in particular H/R = 0.036), or by calculation from parameters found in the respective work ( §). Unmarked values are straightforward calculations from other values. Physical values for the fiducial simulation L1K1 were calculated assuming a body orbiting at 5.2 AU (the same as Benítez-Llambay et al. (2015)) around a solar-mass central body.
Here λc ≡ k −1 c rather than as in Eq. 4 in accordance with Masset (2017)'s definition. Lega et al. (2014) Benítez-Llambay et al. (2015)  torque is expected to provide a dominant contribution to the total migration torque at disk radii where the optical depth drops below a critical value. The heating torque may therefore impact models for the migration, trapping, and growth of objects embedded within AGN disks, and should be considered in future analyses of such systems. In the case where the disk-embedded body is itself an accreting compact object, the mechanical luminosity of outflows may also modify the local density distribution and generate a migration torque (Li et al. 2020).
Using local simulations on a uniform grid, we have been able to verify the analytic predictions for the strength of the heating torque at approximately (in the most favorable cases) the 10% level. More precise tests would be possible using static mesh refinement methods, which would also allow a fuller mapping of how the thermal torque scales with the control parameters in regimes where the assumptions of the analytic theory fail. It would also be valuable to relax the assumptions of a constant energy injection rate and conductive energy transport. Simulations that consistently resolve accretion onto disk-embedded objects, and the radiative feedback that accretion produces, are challenging but are becoming increasingly feasible.

ACKNOWLEDGEMENTS
We thank Pablo Benítez-Llambay for useful discussions and feedback. AMH acknowledges support from the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1650115. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. PJA acknowledges support from NASA TCAN award 80NSSC19K0639. This research was carried out in part during the 2019 Summer School at the Center for Computational Astrophysics, Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. The numerical simulations and analysis were performed on the Rusty cluster of the Flatiron Institute.