Photon Momentum Coupling

# Numerical Problems in Coupling Photon Momentum (Radiation Pressure) to Gas

Philip F. Hopkins, Michael Y. Grudić
TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
E-mail:phopkins@caltech.edu
Submitted to MNRAS, February 2018
###### Abstract

Radiation pressure (RP; or photon momentum absorbed by gas) is important in a tremendous range of astrophysical systems. But we show the usual method for assigning absorbed photon momentum to gas in numerical radiation-hydrodynamics simulations (simply integrating over cell volumes or evaluating at cell/particle centers) can severely under-estimate the correct RP force around sources, unless photon mean-free-paths are highly-resolved in the fluid grid. The existence of this error is independent of the numerical radiation transfer (RT) method (even in exact ray-tracing/Monte-Carlo methods), because it depends on how the RT solution is interpolated back onto fluid elements. Brute-force convergence (resolving mean-free paths) is impossible in many cases (especially where high-opacity UV/ionizing photons are involved). Instead, we show a “face-integrated” method – integrating and applying the momentum fluxes at interfaces between fluid elements – reproduces the correct solution at all resolution levels. The “fix” is simple and we provide example implementations for ray-tracing, Monte-Carlo, and moments RT methods in both grid and mesh-free fluid schemes. We consider an example problem of star formation in a molecular cloud with single-scattering UV/ionizing RP. At state-of-the-art resolution, cell-integrated methods under-estimate the effects of RP by an order of magnitude, leading (incorrectly) to the conclusion that RP is unimportant, while face-integrated methods predict strong self-regulation of star formation and rapid cloud destruction via RP.

###### keywords:
galaxies: formation — galaxies: evolution — galaxies: active — stars: formation — cosmology: theory

## 1 Introduction

Radiation-hydrodynamics (RHD) is fundamental to the behavior of a huge range of astrophysical systems, including galaxies; star-forming giant molecular clouds (GMCs) and star clusters; the inter-stellar (ISM) circum-galactic (CGM) and inter-galactic (IGM) medium; proto-planetary disks (PPDs) and sites of planet formation; stellar structure and evolution; accretion physics, compact objects, and active galactic nuclei (AGN); supernovae; and more. In most of these systems the “radiation pressure forces” – meaning both optically thin transfer of momentum from the radiation field to the gas via absorbed photon momentum, and optically thick (multiple scattering) local momentum transfer generating a true “pressure” – cannot be neglected. However, because these are complicated, non-linear, turbulent, multi-physics systems, numerical simulations are often necessary, and RHD is a computationally challenging problem for many reasons. As such, a great deal of literature has focused on numerical methods for solving (or approximately solving) the radiative transfer equation (RTE) on-the-fly in simulations (for some recent examples, see e.g. Hopkins et al., 2011; Davis et al., 2012; Kuiper et al., 2012; Bate, 2012; Wise et al., 2012; Kolb et al., 2013; Rosdahl et al., 2013; Davis et al., 2014; González et al., 2015; Roth & Kasen, 2015; Tominaga et al., 2015; Buntemeyer et al., 2016; Zhang & Davis, 2017; Rosen et al., 2017; Kim et al., 2017; Foucart, 2018)

But comparably little attention has been paid to how, given a solution to the RTE, radiation pressure forces are coupled onto the gas or fluid (although see e.g. Lowrie et al., 1999). In this paper, we show that the most common methods used in the literature can artificially suppress the radiation pressure forces by orders-of-magnitude, unless the photon mean free paths are very well-resolved. This is especially problematic in the vicinity of sources, so for physical systems with many discrete sources (e.g. star and galaxy formation) can be particularly important. But in many of these contexts, the mean-free-paths of especially ionizing and UV photons in the dense gas around sources (where radiation pressure is most important) can be many orders-of-magnitude smaller than state-of-the-art numerical resolution.

In § 2 we briefly review the equations solved for the RTE and radiation pressure forces, and common classes of numerical implementations. § 3 considers the coupling of photon momentum in the single-scattering limit, and shows how the typical “cell-integrated” approach fails (§ 3.1), and why the resolution requirements to “brute force” convergence with this method are impractical (§ 3.2). We propose instead, in § 3.3, an alternative, face-integrated approach, which resolves these errors and properly treats the absorbed photon momentum around sources even in the limit where photon mean free paths are unresolved. We show how to implement this in various (ray-tracing, Monte-Carlo, moments-based) RTE methods, and various (fixed grid, mesh-free) hydrodynamic methods (additional details in App. A). § 4 generalizes to multiple-scattering, and § 5 demonstrates how these errors can have a dramatic impact on real astrophysical problems, using a simulation of star formation in a molecular cloud as an example. We conclude in § 6.

## 2 Review of RHD Methods

Recall, the radiative transfer equation (RTE), to the order of interest (neglecting higher-order relativistic & Hubble-flow effects unimportant here), is

 (1c∂∂t+^Ω⋅∇)Iν =jν−ανIν (1) −∫dΩ′[Iνσν(Ω,Ω′)−Iν(Ω′)σν(Ω′,Ω)]

(Mihalas & Mihalas, 1984), where is the speed of light, the unit vector in the direction of some differential solid angle , and the absorption and scattering coefficients, the emissivity, and the intensity. Ignoring scattering (or assuming it is isotropic) and sources, this is

 (1c∂∂t+^Ω⋅∇)Iν →−ρκνIν≡−IνλMFP (2)

where is the gas density, the opacity, and the photon mean free path (MFP). Within an infinitesimally small differential volume , the rate of momentum transfer from radiation to gas is:

 d˙pd3x =∂(ρu)∂t=∫dνdΩρκνIν^ΩdΩc=∫dνρκνFνc (3)

where is the flux. We can also define the energy density , and source function .

A range of methods in the literature attempt to solve Eq. 1 “on the fly” in simulations, subject to various approximations. The problems highlighted here apply to all methods, but because details of the necessary correction differ, we briefly review broad classes of popular methods for solving the RTE.

### 2.1 Collisionless (Ray-Tracing & Monte-Carlo) Methods

Ray-tracing with “long characteristics” (RT-LC) and/or Monte Carlo (MC) methods, which explicitly track along discrete “rays” or photon packets from all sources, can (at least in principle) exactly solve the RTE. This is often prohibitively expensive so it is common to make a variety of approximations, e.g. (1) using a limited (finite) number of rays/photon packets; (2) assuming either an infinite-speed-of-light/equilibrium solution (dropping the term) or “reduced speed of light” () in explicit non-equilbrium solutions, both of which allow for larger timesteps; (3) neglecting scattering (or assuming it is isotropic); (4) neglecting velocity-dependence in (for continuum transfer) and/or higher-order relativistic terms; (5) limiting to a small number of discrete frequency “bins” by integrating over some finite ; (6) replacing direct rays from all sources in the volume to all points with “adaptive ray tracing” (ART; Abel & Wandelt 2002) where rays are “split” or “merged” to sample the domain more/less accurately in some regions, or “short-characteristics” (RT-SC; Olson & Kunasz 1987) where only a fixed set of angles within each cell (which must be interpolated between cells) are used; (7) following direct rays but collapsing non-local shadowing/extinction to a local isotropic region around each source (LEBRON; Hopkins et al. 2014, 2017a); (8) treating absorption and re-emission as elastic scattering (in e.g. “implicit Monte Carlo” methods; Fleck & Cummings 1971).

For our purposes, these approximations are not important – these methods are all similar in the key respects that they explicitly follow along different angles , and (at least in principle) can have structure below the fluid grid scale.

### 2.2 “Moments” or “Fluid” Methods

In moments methods, one integrates the RTE over , with , to take the moments. The hierarchy of moment equations never closes (each depends on the next-higher moment), so one truncates the series by adopting an ad-hoc closure ansatz.

Flux-limited diffusion (FLD; Levermore & Pomraning 1981) closes the hierarchy at , giving , with the ansatz with Eddington tensor ( is the identity matrix) and the “flux-limiter” (with ) chosen to interpolate between free-streaming and isotropic-diffusion-like behavior based on an estimate of the local optical depth. The “optically thin variable Eddington tensor” (OTVET; Gnedin & Abel 2001) approach is identical to FLD but with where is the flux that would be present with no obscuration anywhere. The “moment-one” (M1; Levermore 1984) method closes at , keeping the energy equation above and adding the flux equation (so is explicitly evolved), with the ansatz , , , which again interpolates between free-streaming and isotropic diffusion. Other assumptions for the closures are of course possible but change nothing for our study here.

Independent of the closure, the salient feature of these methods is that each moment is single-valued at a given – one has made the fluid approximation for the radiation. This gives the well-known result that moments methods do not, in fact, converge to the correct RTE solutions (e.g. anti-parallel photon streams will “collide”) except in the limit when the optical depths are everywhere infinite. These are popular methods because of their simplicity, but introduce unique complications discussed below.

## 3 Errors & Corrections: Single-Scattering

First consider “single-scattering”: photons free-stream until they are absorbed (no re-emission or scattering). To demonstrate the numerical issues, consider an idealized problem shown in Fig. 1. A single monochromatic isotropic source with luminosity sits at position , centered on a computational domain which we take to be a regular Cartesian grid with cells of side-length (our conclusions are independent of the mesh geometry but this is convenient). Furthermore, assume that (regardless of RT method used) the RTE is solved perfectly at all points in the domain.

### 3.1 The Problem: Cell-Integrated Coupling

In most RHD implementations, regardless of the method used to solve the RTE, the radiation force on the gas is computed by integrating the momentum of absorbed photons over the volume of a given domain (cell) “”:111Even simpler than the “cell-integrated” Eq. 4, some RHD implementations (in e.g. moments methods or some SPH/finite-point methods) adopt a “cell-centered” approach where the acceleration is evaluated at the cell center/particle location and then assigned to the whole cell/particle. As this is equivalent to evaluating Eq. 4 using only the cell-centered (likewise for , ), it produces the same errors as cell-integrated approaches. In addition it is noisier than cell or face-integrated methods, and can violate linear momentum conservation, so we will not discuss this particular case in more detail.

so the total momentum is simply added to the gas in each timestep . For discrete (e.g. Monte Carlo) methods, this integral is given by summing over all absorptions in a given cell: , where is the total photon energy absorbed in photon packet/ray , within domain , over timestep .

Now consider the case where the photon MFP is un-resolved by the fluid grid, i.e. . Essentially all the photons should therefore be absorbed inside the single cell surrounding the source. Eq. 4 then trivially evaluates to because the source emits isotropically.

In the absence of other physics, the system would remain (incorrectly) perfectly static as . The correct solution (with negligible gravity and pressure forces) is that a shell of material moves away from the source, sweeping up gas (leaving an empty cavity) with radial momentum flux (so the total radial momentum ). If and are constant the shell should expand indefinitely with radius .

### 3.2 Solution 1: Increasing Resolution (Is Not Practical)

One solution to the problem above is to increase the resolution of the hydrodynamic grid until the photon mean-free paths are all well-resolved, i.e. .222Another possibility would be to artificially decrease the opacity in different bands, so that mean free paths are always resolved. This amounts to capping the opacity at . This immediately creates a number of problems: (1) Most important, this can lead to photons being absorbed at the wrong physical locations, far further from their sources than they should be. (2) It is very difficult to implement such a prescription in anything but an idealized simulation, without risking allowing photons to escape entirely from dense regions where they should have been absorbed (e.g. “under-shooting” the opacity). (3) The central cell around the source is still not correctly being “swept up” in a shell in the problem described above, because the photons are free-streaming out to neighboring cells. So a shell will form, but only external to a central dense region from which the photons escape – thus the solution at finite resolution will still not correctly represent the converged solution. This would eventually converge to the correct behavior (assuming a perfect RTE solution).

But this is not possible for many real simulations. For example, in simulations of star formation, galaxies, or AGN (outside the accretion disk), the single-scattering absorption comes primarily from neutral hydrogen absorption of ionizing photons, with (for neutral gas), and/or near-UV with , in relatively dense gas around the sources (e.g. , in HII regions, or in the obscuring “torii” around AGN). In a uniform grid, this would require

 Δx≪λMFP∼5×10−6n−14pc (5)

(where ) for ionizing photons in neutral gas. Even for just near-UV photons this gives , so a grid of elements is needed for a pc-radius GMC (even if the maximum density were capped at ). In Lagrangian or AMR methods where the mass resolution is approximately fixed and spatial resolution is adaptive (), the required mass resolution would be

 Δm≪λ2MFP/κν∼3×10−14n−24M\sun (6)

(for ionizing photons in neutral gas) or (for near-UV). This is wildly beyond the resolution of state-of-the-art simulations.

### 3.3 Solution 2: Face-Integrated Coupling

Consider instead a face-integrated momentum coupling. Instead of considering only the volume element , note that is surrounded by a set of faces333In “mesh-free” methods (e.g. SPH, FPM) one can always define “effective” faces by reference to equation-of-motion and point locations: we provide a generic method for this in Appendix A. (between domain and each neighboring element ) each with an oriented vector area . Now simply integrate the absorbed flux moving “towards” each face :

where if the face is the “intercepted face” (i.e. if is the first face crossed, along the ray originating at in direction ), and otherwise. Equivalently, if “points to” face from within . Thus this is the integral over the domain where is the range of solid angle subtended by face from (i.e. angles where ). In (exact or approximate) ray-tracing methods where absorption is calculated along rays with fixed , it is straightforward to evaluate Eq. 7 exactly.

In the optically-thick limit of particular interest, Eq. 7 is approximately . In discrete (Monte-Carlo) methods, this becomes .

Since the momentum flux is defined now at faces, rather than cell centers, it can be incorporated either as a simple operator-split flux (i.e. cell “” receives a momentum from cell “”), or as a force in the Riemann problem or momentum equation solved between cells. Note that in general , so the net momentum flux between cells, , is what will ultimately matter. But for a single source in our example problem, , so the force is purely “outwards” from the source (as it should be), while for a perfectly symmetric or homogeneous source distribution, and the net forces correctly vanish.

Return to our test problem: it is easy to verify that this produces an outward force from the central cell. The momentum flux into each neighbor follows the exact solution (assuming a spherical shell propagates out from the center), averaged over the face. While initially the solution cannot, of course, be exactly spherical (owing to the grid geometry), it will rapidly converge to such as the shell expands.

#### 3.3.1 Implementation in Moments-Based Methods: Additional Complications and Errors

In moments-based methods, two additional, independent and very important errors arise.

1. Because the “mesh” on which radiative properties are computed is identical to the fluid mesh, the flux is single valued within a given cell. So we must insert an explicit model for around each source. Since the error here, where important (), is dominated by the cell surrounding the source, this can be accomplished during the “photon deposition” step. Recall, in moments methods photons must be “deposited” onto the grid around a source, normally via a scalar kernel/weight function: e.g.  where . So during this step, one can explicitly calculate the integral in Eq. 7 with from each source, integrated over domain within cell towards face . A detailed example of how to do this numerically is given in Appendix A. Outside the cell hosting a source, it is a much smaller error to simply adopt the cell-centered (if , the solution is converged, if , the flux outside the origin cell is negligible).

2. In moments methods, the flux does not directly follow from , but is sourced by the numerical gradient of . This means that even outside the “origin” cell, it requires several resolution elements to resolve any gradient in (consider, at the origin, the numerical gradient must vanish if the source is isotropic). So even when , essentially no photon momentum can be transferred to the gas within the central cells in any direction because the numerical gradient is smoothed by a finite discrete “kernel.” This is true even for higher-order gradient estimators in regular Cartesian-mesh codes (see Fig. B1 in Rosdahl et al. 2015). This means a mass has its momentum under-estimated, which makes the resolution criterion in Eqs. 5-6 much more challenging (the mass resolution must be at least times better, or spatial resolution times better, to converge). One fix to this is to extend the explicit integration method for in (i) above, to a -cell radius around each source – essentially, performing a “mini ray-trace” in each sphere around each source. This is the most accurate method, but is computationally much more complex. A much simpler, albeit significantly less accurate, fix is proposed and adopted by Rosdahl et al. (2015): simply replace whenever calculating the photon momentum/radiation pressure terms. This is exact in the single-source, optically-thin limit; however if there are multiple sources with intersecting rays this can over-estimate the true momentum transfer (although the error is usually small unless the sources should exactly cancel).

### 3.4 Numerical Example

Fig. 2 demonstrates the errors and their fixes here with a variant of the simple test problem in Fig. 1. We place a single, isotropic, monochromatic point source of constant luminosity randomly on an effectively infinite grid of uniform cartesian cells of size with constant density and opacity , and calculate the equilibrium flux at all points (assuming the light-crossing time is much shorter than all other timescales in the problem) according to the labeled RHD method. Again we stress the gas properties are “frozen” - we are just solving the RTE. We then calculate the momentum flux either according to the cell-integrated or face-integrated approach. For the moments methods our “face-integrated” method also implements the additional fixes from § 3.3.1. We calculate the total “outward”/radial momentum flux from the source, where is the unit vector pointing from source to the cell center in the cell-integrated case; for the face-integrated case we sum over where points to the face-center. This has the exact solution . We compare for each method, as a function of – in these units, the absolute units of the problem scale out completely.

As expected, when the momentum is deposited in a “cell-integrated” fashion, even exact ray-tracing or MC methods severely suppress the momentum unless . The “cell-integrated,” “default” moments (M1/FLD) methods fare even worse, requiring . With “face-integrated” methods, however, the errors are radically reduced, and depend only very weakly on resolution (owing to geometric effects discussed above).

At low resolution (), even with an exact RT solution, the face-integrated methods converge to a coupled radial momentum , not (i.e. the total radial momentum is lower than the converged solution). This is a real geometric effect.444In fact the rectilinear Cartesian mesh we adopt gives nearly the worst-case scenario for this geometric effect. If we instead assume a glass configuration of mesh-generating points with a Voronoi tesselation or the mesh-free MFM/MFV methods from Hopkins (2015) used to calculate the faces, the larger number of faces (cells are closer to regular polyhedra with faces) means the mean face-integrated momentum as is . Imagine a perfectly-resolved spherical shell emerges from with total radial momentum – but upon entering cell , we must integrate the total momentum entering , which means averaging the momentum of the parts of the shell crossing (producing some cancellation of the momentum components transverse to the face). So the solution is still “exact” in that it reflects the converged solution, “averaged down” to the grid scale after propagating through each cell. More important, this difference is small and constant (while cell-integrated methods have decreasing exponentially as ).

For this simple test problem, approximate-ray (LEBRON/ART) methods give the same result as exact ray-tracing/MC methods. Moments methods always produce larger errors for the reasons in § 3.3.1 (because we allow the radiation field to reach equilibrium, FLD and M1 are identical here). Even with the face-integrated method and fixes in § 3.3.1, in moments methods are not perfect: at they converge (much faster than with cell-integrated coupling), and at the solution occurs entirely “within one cell” so they are identical to the exact methods, but when , the gradient errors noted in § 3.3.1 somewhat suppress the coupled momentum.

## 4 Extension to Multiple-Scattering

Briefly, consider multiple-scattering: photons are absorbed and re-emitted (or scattered) many times before escaping/being destroyed. Where this is of interest, opacities are usually low : e.g. in IR dust re-processing (), free-electron scattering (), X-ray metal-line absorption ( at the Fe edge). So resolution requirements are less extreme: or . But even these are not always met; moreover multiple-scattering can also be important in resonance lines where is much larger.

It is trivial to see that in this limit once again, a “cell-integrated” coupling of photon momentum to gas (Eq. 4) gives essentially vanishing radiation pressure in the cell around a source when (precisely the opposite of the correct behavior). But it is also easy to verify that the “face-integrated” coupling (Eq. 7) gives the correct solution in this limit.

We can repeat the numerical experiment in Fig. 2, for example, but instead of single-scattering assume multiple-scattering with a gray (constant) opacity , take compared to resolved scales, and assume perfect re-emission (all absorbed energy is re-radiated). In that case the exact solution is now , where is the total radial momentum flux integrated over all cells within a radius , and is the optical depth out to that radius. If we then compare to the correct solution , as a function now of (since by assumption), we obtain essentially identical results to those shown in Fig. 2 for single-scattering.

## 5 Effects in an Example Problem: Star Cluster Formation

To explore this in an astrophysical context, we consider as an example a simulation of star cluster formation in an individual GMC following Grudić et al. (2018). The simulations use the code GIZMO (Hopkins, 2015),, with the meshless-finite-mass (MFM) Godunov MHD solver, and fully-adaptive/Lagrangian force-softening. We initialize a solar-metallicity cloud of mass and size , with intentionally low resolution (fixed ), and an initially turbulent velocity and magnetic field spectrum with virial parameter of unity and mean plasma , and evolve it including self-gravity, ideal MHD, radiative cooling (with the cooling curves from Hopkins et al. 2017a from  K, assuming a universal Milky-Way like radiation field), and star formation into sink particles in locally self-gravitating, Jeans-unstable, converging flows (all details in Grudić et al., 2018). While the simulations in Grudić et al. (2018) included stellar mass-loss, SNe (Ia & II), and multi-band RHD with single and multi-scattering RP & photo-heating, here we disable these terms and consider only monochromatic single-scattering RP (to clearly isolate the effects of interest). To further simplify, each star particle (once formed) is assigned a constant light-to-mass ratio and we adopt a constant (approximately appropriate for near-UV luminosities & dust opacities around young populations).

We solve RHD using one of two approximate methods: (1) the ray-based LEBRON (for implementation details, see Hopkins et al. 2017a), or (2) moments-based M1 (following Rosdahl et al. 2013). For each, we consider both “cell-integrated” and “face-integrated” implementations as Figs. 1-2 (with the “face-integrated” M1 also including the Rosdahl et al. 2015 fix to the gradient error discussed in § 3.3.1). We also run a test with no radiation for comparison. We emphasize that both of these methods will have large errors here: LEBRON is exact in the optically-thin limit, M1 in the optically-thick multiple-scattering limit, so neither solves the RTE exactly in this problem where there are both thin & thick single-scattering regimes. So since their errors and regimes of applicability are quite different, it is useful to consider how both methods are influenced by these errors.

Fig. 3 plots the “star formation efficiency” (SFE; fraction of the mass turned into stars): over a few dynamical times, gas collapses and rapidly turns into stars until it is exhausted or expelled by RP. In both cell-centered runs, the RP is strongly suppressed at this resolution and has no effect on the SFE. In both face-centered runs, the SFE is suppressed by a significant factor – i.e. RP has a dramatic effect on the cluster evolution.666The Fig. 3 runs with “fixed” schemes give roughly the expected result from simple analytic arguments: comparing the strength of the total RP force to gravity shows that RP should over-whelm gravity on the cloud scale when . This may explain several apparently discrepant results in the literature. For example, the default Grudić et al. (2018) simulations adopted the “face-integrated” LEBRON method, while a similar study by Raskutti et al. (2016) used the “cell-integrated” (and otherwise “un-corrected”) M1 method. As expected from our test here, Raskutti et al. (2016) find an order-of-magnitude higher SFE (much weaker effects of radiation pressure), for otherwise similar clouds.

It is worth noting that the errors here can be “hidden” in numerical studies which include e.g. both RP and other feedback (e.g. photo-ionization or stellar winds or SNe). Consider: if some other (non-RP) physics first creates a low-opacity “bubble” around the source (e.g. stellar winds sweeping gas into a shell surrounding a very low-density cavity, or SNe destroying dust, or photo-ionization generating a Stromgren sphere which eliminates the neutral hydrogen opacity inside the sphere and/or pushes gas again into a shell), then one can have a mean free path immediately around the source which becomes resolved (), even though when the photons eventually encounter a distant high-opacity “shell” (or neutral gas or dust) they are absorbed in a thin layer ( “in the shell”). In that regime, as long as the shell radius is well-resolved, the errors here are not large (the momenta do not cancel inside a single cell). However, in these cases, the early RP effects (before the shell/Stromgren sphere reaches large radii ) will still be severely suppressed (if ). Moreover, if one takes one of these simulations and “turns off” the other sources of feedback, then RP alone would be dramatically suppressed, leading to the incorrect conclusion that it is un-important.

## 6 Conclusions

We have shown that “cell-integrated” coupling of absorbed photon momentum (radiation pressure) in radiation-hydrodynamics treatments – the most common approach used in the literature – severely under-estimates the true momentum flux around sources, unless the photon mean-free-paths are well-resolved in the hydrodynamic grid (fluid spatial resolution or mass resolution ). But the required resolution is often impossible – for example, in simulations of star or galaxy formation, proper treatment of the UV & ionizing photons (most of the single-scattering radiation pressure) formally requires mass resolution (or in FLD/M1 “moments” methods)!

Fortunately, we show that adopting a “face-integrated” coupling – in which the momentum flux is integrated towards each hydrodynamic face, and treated as part of the usual fluxes or Riemann problem, instead of being integrated over the entire cell volume and added to the cell-centered momentum – resolves these errors. We show that this produces good convergence (provided the radiative transfer equation is properly solved) even when the gas grid fails to resolve the photon mean free paths ().

We stress that the errors identified here, and their fixes, are independent of the radiation methods, even if the RTE is solved perfectly (with e.g. exact Monte-Carlo or ray-tracing methods). The error arises not from the solution of the RTE, but from how the absorbed photon momentum is assigned to the gas. As such the errors and their fixes are extremely general. We do also identify some additional (related but distinct) errors in common implementations unique to moments-based methods (FLD/OTVET/M1) for RHD and discuss potential fixes. We also demonstrate that the errors and fixes are qualitatively independent of the hydrodynamic method – we outline how to implement these for both grid-based codes and mesh-free (e.g. SPH) methods.

We show that the erroneous formulation can fail entirely to capture the effects of radiation pressure (under-estimating the true photon momentum by orders-of-magnitude) when . We illustrate the effects of this in a fully non-linear example of state-of-the-art simulations of star cluster formation, where we show an improper cell-integrated treatment of the single-scattering photon-momentum leads to incorrect conclusions. Specifically, with a correct treatment of the radiation pressure (in either M1 or ray-based methods), radiation pressure from massive stars rapidly disrupts the star-forming GMC and greatly suppresses the star formation efficiency. With an incorrect (cell-integrated) treatment (again in both M1 and ray-based methods) the radiation pressure (erroneously) does very little to the cloud.

We stress that simple “resolution tests” (running the same ICs with increasing resolution) applied to this problem will not reveal the problem with cell-integrated approaches, because the behavior will change sharply (and converge to the correct solution) only when the mass resolution becomes smaller than the ludicrously-small value above (which would require particles in our example star-cluster simulation). At achievable resolution, one will instead see “false convergence” until this threshold is reached.

## Acknowledgments

Support for PFH and MYG was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.

## References

• Abel & Wandelt (2002) Abel, T., & Wandelt, B. D. 2002, MNRAS, 330, L53
• Bate (2012) Bate, M. R. 2012, MNRAS, 419, 3115
• Buntemeyer et al. (2016) Buntemeyer, L., Banerjee, R., Peters, T., Klassen, M., & Pudritz, R. E. 2016, New Astronomy, 43, 49
• Davis et al. (2014) Davis, S. W., Jiang, Y.-F., Stone, J. M., & Murray, N. 2014, ApJ, in press, arXiv:1403.1874
• Davis et al. (2012) Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9
• Fall et al. (2010) Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJL, 710, L142
• Fleck & Cummings (1971) Fleck, Jr., J. A., & Cummings, J. D. 1971, Journal of Computational Physics, 8, 313
• Foucart (2018) Foucart, F. 2018, MNRAS
• Gammie & Ostriker (1996) Gammie, C. F., & Ostriker, E. C. 1996, ApJ, 466, 814
• Gnedin & Abel (2001) Gnedin, N. Y., & Abel, T. 2001, New Astronomy, 6, 437
• González et al. (2015) González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12
• Grudić et al. (2018) Grudić, M. Y., Hopkins, P. F., Faucher-Giguère, C.-A., Quataert, E., Murray, N., & Kereš, D. 2018, MNRAS, 475, 3511
• Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53
• Hopkins et al. (2014) Hopkins, P. F., Keres, D., Onorbe, J., Faucher-Giguere, C.-A., Quataert, E., Murray, N., & Bullock, J. S. 2014, MNRAS, 445, 581
• Hopkins et al. (2011) Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950
• Hopkins et al. (2017a) Hopkins, P. F., et al. 2017a, MNRAS, submitted, arXiv:1702.06148
• Hopkins et al. (2017b) —. 2017b, MNRAS, in press, arXiv:1707.07010
• Jumper & Matzner (2017) Jumper, P. H., & Matzner, C. D. 2017, ArXiv e-prints
• Kim et al. (2017) Kim, J.-G., Kim, W.-T., Ostriker, E. C., & Skinner, M. A. 2017, ApJ, 851, 93
• Kolb et al. (2013) Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80
• Kroupa (2002) Kroupa, P. 2002, Science, 295, 82
• Kuiper et al. (2012) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2012, A&A, 537, A122
• Levermore (1984) Levermore, C. D. 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149
• Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
• Lowrie et al. (1999) Lowrie, R. B., Morel, J. E., & Hittinger, J. A. 1999, ApJ, 521, 432
• Mihalas & Mihalas (1984) Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York, Oxford University Press, 731 p.)
• Olson & Kunasz (1987) Olson, G. L., & Kunasz, P. B. 1987, Journal of Quantitative Spectroscopy and Radiative Transfer, 38, 325
• Raskutti et al. (2016) Raskutti, S., Ostriker, E. C., & Skinner, M. A. 2016, MNRAS, submitted, arXiv:1608.04469
• Rosdahl et al. (2013) Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188
• Rosdahl et al. (2015) Rosdahl, J., Schaye, J., Teyssier, R., & Agertz, O. 2015, MNRAS, 451, 34
• Rosen et al. (2017) Rosen, A. L., Krumholz, M. R., Oishi, J. S., Lee, A. T., & Klein, R. I. 2017, Journal of Computational Physics, 330, 924
• Roth & Kasen (2015) Roth, N., & Kasen, D. 2015, ApJS, 217, 9
• Solomon et al. (1987) Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730
• Tominaga et al. (2015) Tominaga, N., Shibata, S., & Blinnikov, S. I. 2015, ApJS, 219, 38
• Wise et al. (2012) Wise, J. H., Abel, T., Turk, M. J., Norman, M. L., & Smith, B. D. 2012, MNRAS, 427, 311
• Zhang & Davis (2017) Zhang, D., & Davis, S. W. 2017, ApJ, 839, 54

## Appendix A Effective Face Construction and Coupling: Example Methods

As discussed in the text, our face-integrated methods require defining effective faces around any source and integrating the absorbed momentum in the direction of each face (Eq. 7). Regardless of the hydrodynamic grid, in moments-based methods (FLD/OTVET/M1) we must explicitly construct this “effective grid” on a sub-cell scale in order to obtain correct solutions (see § 3.3.1). And in hydrodynamic methods where the grid is anything but regular (e.g. AMR or moving-mesh or MFM/MFV methods), constructing faces and solving the relevant integrals for an arbitrary geometric distribution of mesh-generating points and source location is highly non-trivial. Moreover in SPH or finite-point methods (FPM) no explicit grid exists. We therefore present here a general method which can be used for any set of mesh-generating points (or “particles” in SPH/FPM), to construct a set of effective faces and solve Eq. 7 in the immediate vicinity of a radiation source.

1. Every timestep , for each source (at position ), determine and .

2. Identify “neighboring” gas surrounding the source. In regular-mesh methods this is straightforward; in mesh-free methods one identifies all gas elements within a “search radius” centered on (), and all elements for which is within their search radius ().777The most common approach to define is via a target “neighbor number” , where and is an appropriate kernel function which integrates (over volume) to unity.

3. Construct “effective faces” around the source. There are many possible choices for this, corresponding to different hydrodynamic methods. For example, a Voronoi tesselation (using and all as mesh-generating points). For the kernel-volume decomposition used in MFM/MFV methods (Hopkins, 2015): where and . In SPH, .

4. Solve the integral in Eq. 7, towards each face . Assume and are constant within each (sub-cell-scale) subdomain, and take to be the coordinate origin. Then for single-scattering, along each ray from the source (and we can add all sources independently). The integral can then be expressed as a set of vector weights :

 (˙pν)ab ≈fabsLaνc¯wba (8) ¯wba ≡wba∑c|wca| (9) wba ≡ωba∑+,−∑α(^x±ba)α(fα±)a (10) (fα±)a ≡⎧⎨⎩12⎡⎣1+(∑cωca|^x∓ca|α∑cωca|^x±ca|α)2⎤⎦⎫⎬⎭1/2 (11) ωba ≡12⎛⎜ ⎜⎝1−1√1+(Aba⋅^xba)/(π|xba|2)⎞⎟ ⎟⎠≈ΔΩba4π (12)

where the are the positive or negative (singly-signed) projection vectors:

 ^xba ≡xba|xba|=∑+,−^x±ba (13) (^x+ba)α ≡|xba|−1MAX(xαba,0)∣∣α=x,y,z (14) (^x−ba)α ≡|xba|−1MIN(xαba,0)∣∣α=x,y,z (15)

and, for single-scattering, . For multiple-scattering from a single source, the expression above is identical up to – if we assume a gray opacity and perfect re-radiation (within the single cell sub-domain) then .

These expressions are derived in Hopkins et al. (2017b) (up to the trivial addition here of calculating ); while non-trivial they have three key properties. (1) They maintain manifest conservation of linear momentum. (2) They give fluxes which are statistically isotropic in the rest frame of the source (i.e. the coupled momenta are not numerically biased in any particular direction, regardless of the position of the source within the cell or the position/distribution of mesh-generating points/particles). (3) They approximate, as closely as possible without a computational expensive exact numerical quadrature, the exact integral of Eq. 7.

5. Verify conservation: one should ensure total momentum and energy are conserved, regardless of how Eq. 7 is solved. For an isotropic source in a uniform density field, this means ensuring the total linear momentum coupled vanishes (the momentum-coupling is symmetric). If one solves Eq. 7 exactly at all points, and has an exact solution for at all points, and has a set of faces which form an exactly-closed convex hull, then this is guaranteed mathematically, but numerically these conditions are not usually all satisfied. Here, this step is included above: the terms in Eq. 11 are numerically evaluated in a first-loop over the neighbors, to calculate , and represent the “correction” (vector re-normalization) needed to guarantee conservation (one can verify that if the conditions above were met, then one would always have exactly).

6. Couple the momentum to the gas fully-conservatively, either as direct fluxes/forces into each neighbor , or in a Riemann problem.

The procedure above resolves the primary issues in the text, on arbitrary grid, for arbitrary source positions. Note that we have calculated the face-integrated coupling in the vicinity of the source. Far away from the source, one can extend a method like this if desired. However, far from a source, the direction of the flux is not strongly varying within a single cell – so simply using the cell-centered average flux (equivalently, assuming the net flux at all positions inside a cell , far from any sources, points in approximately the same direction) is not a significant source of error in this regime.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters