Numerical Problems in Coupling Photon Momentum (Radiation Pressure) to Gas
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 radiationhydrodynamics simulations (simply integrating over cell volumes or evaluating at cell/particle centers) can severely underestimate the correct RP force around sources, unless photon meanfreepaths are highlyresolved in the fluid grid. The existence of this error is independent of the numerical radiation transfer (RT) method (even in exact raytracing/MonteCarlo methods), because it depends on how the RT solution is interpolated back onto fluid elements. Bruteforce convergence (resolving meanfree paths) is impossible in many cases (especially where highopacity UV/ionizing photons are involved). Instead, we show a “faceintegrated” 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 raytracing, MonteCarlo, and moments RT methods in both grid and meshfree fluid schemes. We consider an example problem of star formation in a molecular cloud with singlescattering UV/ionizing RP. At stateoftheart resolution, cellintegrated methods underestimate the effects of RP by an order of magnitude, leading (incorrectly) to the conclusion that RP is unimportant, while faceintegrated methods predict strong selfregulation of star formation and rapid cloud destruction via RP.
keywords:
galaxies: formation — galaxies: evolution — galaxies: active — stars: formation — cosmology: theory1 Introduction
Radiationhydrodynamics (RHD) is fundamental to the behavior of a huge range of astrophysical systems, including galaxies; starforming giant molecular clouds (GMCs) and star clusters; the interstellar (ISM) circumgalactic (CGM) and intergalactic (IGM) medium; protoplanetary 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, nonlinear, turbulent, multiphysics 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) onthefly 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 ordersofmagnitude, unless the photon mean free paths are very wellresolved. 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 meanfreepaths of especially ionizing and UV photons in the dense gas around sources (where radiation pressure is most important) can be many ordersofmagnitude smaller than stateoftheart 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 singlescattering limit, and shows how the typical “cellintegrated” 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, faceintegrated 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 (raytracing, MonteCarlo, momentsbased) RTE methods, and various (fixed grid, meshfree) hydrodynamic methods (additional details in App. A). § 4 generalizes to multiplescattering, 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 higherorder relativistic & Hubbleflow effects unimportant here), is
(1)  
(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
(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:
(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 (RayTracing & MonteCarlo) Methods
Raytracing with “long characteristics” (RTLC) 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 infinitespeedoflight/equilibrium solution (dropping the term) or “reduced speed of light” () in explicit nonequilbrium solutions, both of which allow for larger timesteps; (3) neglecting scattering (or assuming it is isotropic); (4) neglecting velocitydependence in (for continuum transfer) and/or higherorder 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 “shortcharacteristics” (RTSC; 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 nonlocal shadowing/extinction to a local isotropic region around each source (LEBRON; Hopkins et al. 2014, 2017a); (8) treating absorption and reemission 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 nexthigher moment), so one truncates the series by adopting an adhoc closure ansatz.
Fluxlimited diffusion (FLD; Levermore & Pomraning 1981) closes the hierarchy at , giving , with the ansatz with Eddington tensor ( is the identity matrix) and the “fluxlimiter” (with ) chosen to interpolate between freestreaming and isotropicdiffusionlike 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 “momentone” (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 freestreaming 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 singlevalued at a given – one has made the fluid approximation for the radiation. This gives the wellknown result that moments methods do not, in fact, converge to the correct RTE solutions (e.g. antiparallel 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: SingleScattering
First consider “singlescattering”: photons freestream until they are absorbed (no reemission 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 sidelength (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: CellIntegrated 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) “”:^{1}^{1}1Even simpler than the “cellintegrated” Eq. 4, some RHD implementations (in e.g. moments methods or some SPH/finitepoint methods) adopt a “cellcentered” 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 cellcentered (likewise for , ), it produces the same errors as cellintegrated approaches. In addition it is noisier than cell or faceintegrated methods, and can violate linear momentum conservation, so we will not discuss this particular case in more detail.
(4) 
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 unresolved 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 meanfree paths are all wellresolved, i.e. .^{2}^{2}2Another 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. “undershooting” 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 freestreaming 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 singlescattering absorption comes primarily from neutral hydrogen absorption of ionizing photons, with (for neutral gas), and/or nearUV 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
(5) 
(where ) for ionizing photons in neutral gas. Even for just nearUV photons this gives , so a grid of elements is needed for a pcradius 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
(6) 
(for ionizing photons in neutral gas) or (for nearUV). This is wildly beyond the resolution of stateoftheart simulations.
3.3 Solution 2: FaceIntegrated Coupling
Consider instead a faceintegrated momentum coupling. Instead of considering only the volume element , note that is surrounded by a set of faces^{3}^{3}3In “meshfree” methods (e.g. SPH, FPM) one can always define “effective” faces by reference to equationofmotion 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 :
(7) 
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) raytracing methods where absorption is calculated along rays with fixed , it is straightforward to evaluate Eq. 7 exactly.
In the opticallythick limit of particular interest, Eq. 7 is approximately . In discrete (MonteCarlo) methods, this becomes .
Since the momentum flux is defined now at faces, rather than cell centers, it can be incorporated either as a simple operatorsplit 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 MomentsBased Methods: Additional Complications and Errors
In momentsbased methods, two additional, independent and very important errors arise.

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 cellcentered (if , the solution is converged, if , the flux outside the origin cell is negligible).

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 higherorder gradient estimators in regular Cartesianmesh codes (see Fig. B1 in Rosdahl et al. 2015). This means a mass has its momentum underestimated, which makes the resolution criterion in Eqs. 56 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 raytrace” 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 singlesource, opticallythin limit; however if there are multiple sources with intersecting rays this can overestimate 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 lightcrossing 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 cellintegrated or faceintegrated approach. For the moments methods our “faceintegrated” 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 cellintegrated case; for the faceintegrated case we sum over where points to the facecenter. 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 “cellintegrated” fashion, even exact raytracing or MC methods severely suppress the momentum unless . The “cellintegrated,” “default” moments (M1/FLD) methods fare even worse, requiring . With “faceintegrated” 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 faceintegrated 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.^{4}^{4}4In fact the rectilinear Cartesian mesh we adopt gives nearly the worstcase scenario for this geometric effect. If we instead assume a glass configuration of meshgenerating points with a Voronoi tesselation or the meshfree 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 faceintegrated momentum as is . Imagine a perfectlyresolved 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 cellintegrated methods have decreasing exponentially as ).
For this simple test problem, approximateray (LEBRON/ART) methods give the same result as exact raytracing/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 faceintegrated method and fixes in § 3.3.1, in moments methods are not perfect: at they converge (much faster than with cellintegrated 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 MultipleScattering
Briefly, consider multiplescattering: photons are absorbed and reemitted (or scattered) many times before escaping/being destroyed. Where this is of interest, opacities are usually low : e.g. in IR dust reprocessing (), freeelectron scattering (), Xray metalline absorption ( at the Fe edge). So resolution requirements are less extreme: or . But even these are not always met; moreover multiplescattering can also be important in resonance lines where is much larger.
It is trivial to see that in this limit once again, a “cellintegrated” 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 “faceintegrated” coupling (Eq. 7) gives the correct solution in this limit.
We can repeat the numerical experiment in Fig. 2, for example, but instead of singlescattering assume multiplescattering with a gray (constant) opacity , take compared to resolved scales, and assume perfect reemission (all absorbed energy is reradiated). 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 singlescattering.
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),^{5}^{5}5http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html, with the meshlessfinitemass (MFM) Godunov MHD solver, and fullyadaptive/Lagrangian forcesoftening. We initialize a solarmetallicity 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 selfgravity, ideal MHD, radiative cooling (with the cooling curves from Hopkins et al. 2017a from K, assuming a universal MilkyWay like radiation field), and star formation into sink particles in locally selfgravitating, Jeansunstable, converging flows (all details in Grudić et al., 2018). While the simulations in Grudić et al. (2018) included stellar massloss, SNe (Ia & II), and multiband RHD with single and multiscattering RP & photoheating, here we disable these terms and consider only monochromatic singlescattering RP (to clearly isolate the effects of interest). To further simplify, each star particle (once formed) is assigned a constant lighttomass ratio and we adopt a constant (approximately appropriate for nearUV luminosities & dust opacities around young populations).
We solve RHD using one of two approximate methods: (1) the raybased LEBRON (for implementation details, see Hopkins et al. 2017a), or (2) momentsbased M1 (following Rosdahl et al. 2013). For each, we consider both “cellintegrated” and “faceintegrated” implementations as Figs. 12 (with the “faceintegrated” 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 opticallythin limit, M1 in the opticallythick multiplescattering limit, so neither solves the RTE exactly in this problem where there are both thin & thick singlescattering 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 cellcentered runs, the RP is strongly suppressed at this resolution and has no effect on the SFE. In both facecentered runs, the SFE is suppressed by a significant factor – i.e. RP has a dramatic effect on the cluster evolution.^{6}^{6}6The 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 overwhelm 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 “faceintegrated” LEBRON method, while a similar study by Raskutti et al. (2016) used the “cellintegrated” (and otherwise “uncorrected”) M1 method. As expected from our test here, Raskutti et al. (2016) find an orderofmagnitude 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. photoionization or stellar winds or SNe). Consider: if some other (nonRP) physics first creates a lowopacity “bubble” around the source (e.g. stellar winds sweeping gas into a shell surrounding a very lowdensity cavity, or SNe destroying dust, or photoionization 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 highopacity “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 wellresolved, 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 unimportant.
6 Conclusions
We have shown that “cellintegrated” coupling of absorbed photon momentum (radiation pressure) in radiationhydrodynamics treatments – the most common approach used in the literature – severely underestimates the true momentum flux around sources, unless the photon meanfreepaths are wellresolved 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 singlescattering radiation pressure) formally requires mass resolution (or in FLD/M1 “moments” methods)!
Fortunately, we show that adopting a “faceintegrated” 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 cellcentered 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 MonteCarlo or raytracing 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 momentsbased 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 gridbased codes and meshfree (e.g. SPH) methods.
We show that the erroneous formulation can fail entirely to capture the effects of radiation pressure (underestimating the true photon momentum by ordersofmagnitude) when . We illustrate the effects of this in a fully nonlinear example of stateoftheart simulations of star cluster formation, where we show an improper cellintegrated treatment of the singlescattering photonmomentum leads to incorrect conclusions. Specifically, with a correct treatment of the radiation pressure (in either M1 or raybased methods), radiation pressure from massive stars rapidly disrupts the starforming GMC and greatly suppresses the star formation efficiency. With an incorrect (cellintegrated) treatment (again in both M1 and raybased 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 cellintegrated approaches, because the behavior will change sharply (and converge to the correct solution) only when the mass resolution becomes smaller than the ludicrouslysmall value above (which would require particles in our example starcluster 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 TGAST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD167592.
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., FaucherGiguè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., FaucherGiguere, 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 eprints
 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 faceintegrated 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 momentsbased methods (FLD/OTVET/M1) we must explicitly construct this “effective grid” on a subcell 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 movingmesh or MFM/MFV methods), constructing faces and solving the relevant integrals for an arbitrary geometric distribution of meshgenerating points and source location is highly nontrivial. Moreover in SPH or finitepoint methods (FPM) no explicit grid exists. We therefore present here a general method which can be used for any set of meshgenerating 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.

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

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

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 meshgenerating points). For the kernelvolume decomposition used in MFM/MFV methods (Hopkins, 2015): where and . In SPH, .

Solve the integral in Eq. 7, towards each face . Assume and are constant within each (subcellscale) subdomain, and take to be the coordinate origin. Then for singlescattering, 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 :
(8) (9) (10) (11) (12) where the are the positive or negative (singlysigned) projection vectors:
(13) (14) (15) and, for singlescattering, . For multiplescattering from a single source, the expression above is identical up to – if we assume a gray opacity and perfect reradiation (within the single cell subdomain) then .
These expressions are derived in Hopkins et al. (2017b) (up to the trivial addition here of calculating ); while nontrivial 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 meshgenerating points/particles). (3) They approximate, as closely as possible without a computational expensive exact numerical quadrature, the exact integral of Eq. 7.

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 momentumcoupling 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 exactlyclosed 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 firstloop over the neighbors, to calculate , and represent the “correction” (vector renormalization) needed to guarantee conservation (one can verify that if the conditions above were met, then one would always have exactly).

Couple the momentum to the gas fullyconservatively, 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 faceintegrated 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 cellcentered 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.