Pebble dynamics and accretion onto rocky planets.
II. Radiative models
Abstract
We investigate the effects of radiative energy transfer on a series of nestedgrid, highresolution hydrodynamic simulations of gas and particle dynamics in the vicinity of an Earthmass planetary embryo. We include heating due to the accretion of solids and the subsequent convective motions. Using a constant embryo surface temperature, we show that radiative energy transport results in a tendency to reduce the entropy in the primordial atmosphere, but this tendency is alleviated by an increase in the strength of convective energy transport, triggered by a correspondingly increased superadiabatic temperature gradient. As a consequence, the amplitude of the convective motions increase by roughly an order of magnitude in the vicinity of the embryo. In the cases investigated here, where the optical depth towards the disk surface is larger than unity, the reduction of the temperature in the outer parts of the Hill sphere relative to cases without radiative energy transport is only 100K, while the mass density increase is on the order of a factor of two in the inner parts of the Hill sphere. Our results demonstrate that, unless unrealistically low dust opacities are assumed, radiative cooling in the context of primordial rocky planet atmospheres can only become important after the disk surface density has dropped significantly below minimummasssolarnebula values.
keywords:
planets and satellites: formation, terrestrial planets – protoplanetary discs – hydrodynamics (HD) – radiative transfer1 Introduction
Planetary embryos of sufficient mass embedded in optically thick protoplanetary accretion disks are necessarily accompanied by extended primordial atmospheres in near hydrostatic equilibrium, merging smoothly with the disk background at distances on the order of the Hill radius, , where is the semimajor axis of the embryo’s orbit, and and are the masses of the embryo and the central star, respectively (Perri & Cameron, 1974; Bodenheimer & Pollack, 1986; Pollack et al., 1996; D’Angelo & Bodenheimer, 2013; Alibert, 2017).
The existence of dust and solid particles in the protoplanetary disk leads to accretion of solids onto such embryos, resulting in growth of the embryo, accretion heating at the embryo surface, and frictional heating of the surrounding primordial atmosphere (e.g. Brouwers et al., 2018). As long as the embryo atmosphere and surrounding disk remain optically thick, convective motions are the most effective means of transporting heat through the atmosphere, with efficient convective transport leading to a nearadiabatic stratification of the atmosphere (Stevenson, 1982; Wuchterl, 1993). Nearadiabaticity implies that the inner parts of the atmosphere are relatively hot, which, when combined with the reduction of the embryo’s gravity as the square of the distance, results in extended, relatively low mass atmospheres (e.g. Venturini et al., 2015).
If, on the other hand, the optical depth of the primordial atmosphere is not very large, radiative energy transport leads to cooling, resulting in possibly significant changes to the stratification of the atmosphere. Specifically, radiative cooling would lead to smaller pressures and density scale heights, thus increasing the mass of the atmosphere. Significant cooling would eventually lead to a situation where hydrostatic equilibrium is no longer sustainable, resulting in gravitational collapse and the formation of gas giants (e.g. Mizuno et al., 1978; Bodenheimer & Pollack, 1986; Hori & Ikoma, 2011). Opposing this effect is the replenishment of gas inside the Hill sphere with disk material (Ormel et al., 2015; Cimerman et al., 2017).
In cases where catastrophic collapse of the atmosphere does not occur, the ultimate formation outcome depends crucially on the balance between loss of atmospheric gas caused by the reduction of the density and pressure of the surrounding disk over time, and the opposing tendency caused by the increased radiative cooling—occurring for essentially the same reason. As shown by Ginzburg et al. (2016), this balance may, in the end, dictate if the final outcome is a gasrich planet, or a rocky planet with a thin remnant atmosphere.
Properly accounting for the radiative transfer (RT) of energy and the corresponding effects on the nearhydrostatic equilibrium of primordial atmospheres is thus crucial for realistic and accurate modelling of early planet formation. Such efforts are additionally complicated by effects of scattering, which are known to be important at the wavelengths and temperatures of relevance (e.g. Pinte et al., 2009).
In this Letter, we present firstofakind results, using ray tracing RT with scattering to demonstrate its effects on primordial atmospheres, and the resulting effects on particle dynamics and pebble accretion onto rocky planets.
2 Methods
As in Popovas et al. (2018, hereafter PNRO18), this study is carried out using the DISPATCH framework (Nordlund et al., 2018), employing a threedimensional, Cartesian (shearing box) domain, with a set of static, nested patches. We continue to use an ideal gas equation of state (EOS) with adiabatic index = 1.4 and molecular weight = 2, and embedment in a disk with a nominal surface mass density of 170 g cm, corresponding to 1/10 of the minimummasssolarnebula (MMSN). The grid set up is the same as in PNRO18 for a planet at 1 AU distance from the central star. We conduct a series of radiativeconvective simulations, for which the basic RT and accretion heating parameters are summarized in Table 1.
To investigate the effects of radiative cooling while maintaining a nearly adiabatic (i.e. convective) atmosphere, we deliberately choose a combination of disk surface density and opacity such that the disk and embryo atmosphere are slightly optically thick. To this end we adopt a total opacity (with 80% scattering and 20% absorption) of 0.1 cm g, and a disk surface density of 170 g cm. We thus obtain a midplane optical depth in the unperturbed disk of . The additional optical depth of the initial, adiabatic atmosphere is , and hence the optical properties are such that we expect to begin to see the effects of radiative cooling.
The optical properties of dust and gas in hot, primordial atmospheres are quite uncertain, mainly due to uncertainties related to dust coagulation and thermal processing. Nevertheless, if we consider the opacities from Semenov et al. (2003) at densities and temperatures relevant to this work (e.g. Fig. 1), then the vast majority of our domain has opacities on the order of 0.1 – 1 cm g.
2.1 Initial and boundary conditions
For initial conditions we take a fiducial, fully relaxed MMSN/10 run with adiabatic stratification (m095t10 in PNRO18). The external and spherical boundary conditions for density , entropy per unit mass and mass flux are also the same as in PNRO18. However, as we also consider radiative energy transport in the current work, appropriate boundary conditions must be considered. The disk is optically thick in the radial and azimuthal directions and is initially in radiative equilibrium in those directions, i.e., no heat exchange occurs radially or azimuthally. However, protoplanetary disks do cool radiatively in the vertical direction. Therefore, we adopt the following external boundary conditions for the radiative transfer: at external horizontal (radial and azimuthal) boundaries; at external vertical boundaries, where is the radiation intensity in a specific direction, is the local source function, and is the heating/cooling rate.
Run  [cm g]  [M yr]  

m095  –  –  – 
m095conv2e6  –  –  2 10 
m095conv2e6rt  0.1  0.2  2 10 
m095conv2e5rt  0.1  0.2  2 10 
m095conv2e6rt1  1.0  0.2  2 10 
2.1.1 Heating from the embryo
Accretion heating due to solids is included via a source term in the entropy equation (PNRO18), , where , is the distance from the embryo, is the gravitational constant and is the adopted accretion rate (Tab. 1).
As the envelope surrounding the embryo is optically thick for most of the disk evolution, the embryo and its atmosphere cannot radiate heat away effectively when solids are accreted. The atmosphere is thus expected to remain nearly adiabatic through most of the buildup of the embryo mass, and the embryo must be correspondingly hot, especially towards the end of the buildup period. Detailed predictions of the temperature and its history requires realistic equations of state for both the primordial atmosphere and the planet, in combination with modelling of the embryo thermal evolution over disk evolution time scales. This is a formidable task by itself, where sequences of simulations such as the ones reported here may be used to provide the required estimates of instantaneous cooling rates.
In any case, since the heat capacity of a planet is large, the planet is expected to remain hot for a considerable time, even after the envelope becomes optically thin (Ginzburg et al., 2016). For the relatively brief periods of time covered by the current simulations, it is thus appropriate to consider the surface temperature of the embryo to be fixed, and we therefore adopt the surface temperature predicted by the adiabatic atmosphere model as a fixed lower boundary condition in all models.
2.2 Radiative energy transport
We use a hybridcharacteristics ray tracing scheme (Nordlund et al., 2018, Appendix A, available online only) with 26 raydirections (forward and reverse directions along 3 axes, 6 face diagonals, and 4 space diagonals) and a single frequency bin (constant opacity ). Since scattering is an important effect at the temperatures and wavelengths of relevance here, we split the opacity into an absorption part, , and a scattering part, . At the levels of scattering typically assumed—from 50% to 80%—and with the relatively slow time evolution dictated by convective motions, scattering can be handled with what is essentially “ iteration”, i.e., iteratively feeding the mean intensity from previous time steps into the source function of the next time step (Hubeny, 2003). Here, the time slices stored in DISPATCH may be used to predict the mean intensity in the next time step, thus efficiently reducing the time lag otherwise resulting from the use of the previous mean intensity.
In order to avoid that the intense radiation from the embryo surface results in narrow, parallel beams of strong radiation in the 26 angular directions used in the raybased radiation solver, we replace the solution from the raybased solver with a diffusion approximation in the region near the embryo. In the optically thick limit, the integral solution of the radiative transfer equation along any direction may be approximated by
(1) 
where and are the specific radiation intensities in the forward and reverse directions, and is the optical depth along the ray direction. A term proportional to the first derivative of the source function is omitted, since it cancels out when averaging forward and reverse solutions along a given ray direction. Given the availability of the source function and opacity values on the Cartesian mesh, it is trivial to evaluate this expression in the three axis directions, and estimate the full space angle, integrated heat exchange rate per unit mass:
(2) 
The heat exchange rate for the raybased radiation solver, after averaging over the 26 forward and reverse ray directions, is meanwhile:
(3) 
We apply the raybased solver everywhere, but we weight the resultant heat exchange rate by
(4) 
and the diffusionbased value of from Eq. 2 by , where is 10 times the radius of the embryo.
2.3 Particles
We use 12 million macroparticles, each representing a swarm of real particles with a given size and mass. The initial spatial distribution of macroparticles is proportional to the local gas density, with particle sizes ranging from 10 m to 1 cm, and a constant number of macroparticles per logarithmic size bin. Rather than making assumptions about the settling and actual size distribution, we instead analyse subpopulations of our initial distribution. Specifically, to measure accretion rate as a function of particle size, we tag and follow only the particles that initially reside within one Hill radius of the midplane. See PNRO18 for more details about the particle distribution, their motion, particle size distribution within macroparticles, and selection.
3 Results and Discussion
3.1 Gas dynamics
Figure 1 shows the thermal and density structures of the envelopes for the five cases considered.
Convective motions effectively transport the excess heating from accretion of solids, even when RT is not considered (or is sufficiently high), and therefore the differences in thermal structure relative to the purely convective case are small. When radiative energy transport is used, radiative cooling effects become important in the outer parts of the atmosphere and the temperature decreases. However, cooling effects do not penetrate all the way to the embryo, and the thermal structure near the embryo remains close to the adiabatic one. As the envelope adjusts its nearhydrostatic equilibrium to the lower temperature, it becomes on average denser. Since this is an effect that accumulates with depth—pressure scale height is proportional to the local temperature—the effect on the density is larger than on the temperature. Indeed, in the radiativelycooled models, the total mass of the atmosphere is increased by about 4% and 27% for = 2 10 and 2 10 M yr, respectively (with = 0.1 cm g), relative to adiabatic and purely convective models.
Figure 2 (middle) shows an example of the convective velocity flow for a simulation where radiative energy transfer is included (m095conv2e6rt). The fluid motions are similar to the case with pure convection (left panel of Fig. 2), but the velocity amplitudes are larger—especially at small radii—when radiative energy transfer is included. The reason for the increased velocity amplitudes is illustrated in the right panel of Fig. 2. The radiative cooling tends to lower the entropy in the neighbourhood of the embryo relative to runs without RT, while accretion heating simultaneously keeps adding entropy, and the temperature of the embryo surface is held fixed. The net effect is a tendency to increase the temperature gradient, thus making it more superadiabatic than in the case with only accretion heating, resulting in stronger driving of the convection. Near the embryo, the increase of the convective energy transport is largely able to compensate for the increased radiative cooling. Further out, where the convection motions are weaker and the vertical optical depth is smaller, the effects of radiative cooling become more visible.
The increase of convective velocity amplitudes is illustrated in Fig. 3 against a background of particle drift speeds. The figure shows that velocity amplitudes near the embryo for simulations that include RT are increased by nearly an order of magnitude over the simulation that does not, becoming comparable to the vertical drift speed of 23 mm particles. Conversely, however, increasing the accretion heating further does not significantly increase the dispersion. Further out, where the horseshoe and shear flow patterns start to contribute to the velocity dispersion, the effects of radiative cooling on the velocity dispersion diminish.
3.2 Solid accretion rates
In our previous paper (PNRO18), we determined the accretion rates of solids and showed that they scale linearly with the particle size. In Figure 4, we affirm that the previously determined accretion rates are robust. Indeed, the accretion rates in simulation runs with convective motions and radiative cooling do not differ significantly from runs without; deviations are mainly due to increased noise as the convective motions of the gas stir the solids. Moreover, we find that stronger convective motions due to a larger accretion heating (e.g. that could result from changing the particle size distribution) do not significantly affect the accretion rates. Although the amplitude of the convective motions increase with accretion heating, they remain confined to a region (Fig. 3). Whether a particle will be accreted or not is, meanwhile, determined at larger radii (cf. Fig. 19 of PNRO18), and thus the strength of the convection does not significantly affect the accretion rate of particles.
4 Conclusions and Outlook
In these firstofakind, threedimensional, radiativeconvective models of hot and extended primordial atmospheres of Earthmass embryos, with realistic—albeit schematic—scattering included, we find that radiative cooling leads to a local increase of the shell averaged mass density on the order of a factor of two , but only a slight decrease in the temperature. For the conditions modelled here, cooling effects do not penetrate to the bottom of the atmosphere where convection is the dominant mechanism of energy transport. In response to the tendency from cooling to increase the radial temperature gradient, however, the amplitude of convective motions increase by nearly an order of magnitude near the embryo surface relative to runs without RT.
Even though radiative cooling modifies the atmosphere structure, the accretion heating and resulting convective energy transport still dominate near the planetary embryo, where the atmosphere temperature remains “anchored” to the surface temperature of the embryo. This is encouraging, since it implies that planetary embryos embedded in protoplanetary disks can retain hot atmospheres by adjusting their nearly adiabatic stratifications only slightly, throughout much of the evolution of the disk. For cooling to become important at earlier times, the opacities would have to be much lower than assumed here (0.1  1 g cm).
When the disk and primordial atmosphere finally become genuinely optically thin, the atmospheres that remain around low mass planetary embryos are therefore relatively light, with their future fate depending on the relative balance between the slow cooling of the still hot embryo, and the secular loss of the remaining atmosphere. As argued by Ginzburg et al. (2016), it is the outcome of this balance that ultimately determines if a planet ends up as a gasrich planet with a significant H+He atmosphere, or becomes a rocky planet, with an atmosphere consisting of heavier gas molecules, possibly even dominated by outgassing.
The current modelling should be seen as a pilot effort, exploring effects and probing what is currently possible. As briefly discussed in Appendix B (available online only), the current results were obtained with a moderate computational effort. Obvious factors that could use improvement are the EOS, the constant and grey opacity, and the spatial resolution, which could stand to be increased in order to better resolve the convective motions. Improving all these factors is certainly possible, and will be the subject of future work.
Acknowledgements
We are grateful to the anonymous reviewer for a constructive report that helped to significantly improve the quality of this manuscript. The work of AP and ÅN was supported by grant 132300199B from the Danish Council for Independent Research (DFF). The Centre for Star and Planet Formation is funded by the Danish National Research Foundation (DNRF97). Storage and computing resources at the University of Copenhagen HPC centre, funded in part by the Villum Foundation (VKR023406), were used to carry out the simulations.
References
 Alibert (2017) Alibert Y., 2017, A&A, 606, A69
 Amanatides & Woo (1987) Amanatides J., Woo A., 1987, in In Eurographics ’87. pp 3–10
 Bodenheimer & Pollack (1986) Bodenheimer P., Pollack J. B., 1986, Icarus, 67, 391
 Brouwers et al. (2018) Brouwers M. G., Vazan A., Ormel C. W., 2018, A&A, 611, A65
 Cen (2002) Cen R., 2002, ApJS, 141, 211
 Cimerman et al. (2017) Cimerman N. P., Kuiper R., Ormel C. W., 2017, MNRAS, 471, 4662
 D’Angelo & Bodenheimer (2013) D’Angelo G., Bodenheimer P., 2013, ApJ, 778, 77
 Dullemond & Turolla (2000) Dullemond C. P., Turolla R., 2000, A&A, 360, 1187
 Ginzburg et al. (2016) Ginzburg S., Schlichting H. E., Sari R., 2016, ApJ, 825, 29
 González et al. (2015) González M., Vaytet N., Commerçon B., Masson J., 2015, A&A, 578, A12
 Heinemann et al. (2006) Heinemann T., Dobler W., Nordlund Å., Brandenburg A., 2006, A&A, 448, 731
 Hori & Ikoma (2011) Hori Y., Ikoma M., 2011, MNRAS, 416, 1419
 Hubeny (2003) Hubeny I., 2003, in Hubeny I., Mihalas D., Werner K., eds, Astronomical Society of the Pacific Conference Series Vol. 288, Stellar Atmosphere Modeling. p. 17
 Levermore & Pomraning (1981) Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321
 Mihalas & Mihalas (1984) Mihalas D., Mihalas B. W., 1984, Foundations of radiation hydrodynamics
 Mizuno et al. (1978) Mizuno H., Nakazawa K., Hayashi C., 1978, Progress of Theoretical Physics, 60, 699
 Nordlund (1982) Nordlund A., 1982, A&A, 107, 1
 Nordlund et al. (2018) Nordlund Å., Ramsey J. P., Popovas A., Küffmeier M., 2018, MNRAS, 477, 624
 Ormel et al. (2015) Ormel C. W., Shi J.M., Kuiper R., 2015, MNRAS, 447, 3512
 Perri & Cameron (1974) Perri F., Cameron A. G. W., 1974, Icarus, 22, 416
 Pinte et al. (2009) Pinte C., Harries T. J., Min M., Watson A. M., Dullemond C. P., Woitke P., Ménard F., DuránRojas M. C., 2009, A&A, 498, 967
 Pollack et al. (1996) Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62
 Popovas et al. (2018) Popovas A., Nordlund Å., Ramsey J. P., Ormel C. W., 2018, MNRAS,
 Rijkhorst et al. (2006) Rijkhorst E.J., Plewa T., Dubey A., Mellema G., 2006, A&A, 452, 907
 Robitaille (2011) Robitaille T. P., 2011, A&A, 536, A79
 Rosdahl & Teyssier (2015) Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380
 Semenov et al. (2003) Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611
 Stevenson (1982) Stevenson D. J., 1982, Planet. Space Sci., 30, 755
 Stone et al. (1992) Stone J. M., Mihalas D., Norman M. L., 1992, ApJS, 80, 819
 Venturini et al. (2015) Venturini J., Alibert Y., Benz W., Ikoma M., 2015, A&A, 576, A114
 Wuchterl (1993) Wuchterl G., 1993, Icarus, 106, 323
Appendix A The raybased radiation solver
Being a 7dimensional problem (three space + time + wavelength + unit direction vector (polar and azimuthal angles), RT can be a daunting physical process to solve. A number of approaches have been developed over the last several decades in an attempt to tackle the problem. One common approach is to reduce the dimensionality of the problem. These approximations are known as moment methods, since they take the moments of radiative transfer equation (i.e. the angular dependence is averaged out). It can be a good approximation for smooth, diffuse radiation fields in optically thick media when the radiation is tightly coupled to the gas. Common approximations of various complexity include flux limited diffusion (FLD) (e.g. Levermore & Pomraning 1981; González et al. 2015), the M1 method (e.g. Rosdahl & Teyssier 2015), or variable tensor Eddington factors (Dullemond & Turolla, 2000). Although protoplanetary disks and the deep interiors of the planetary atmospheres are generally optically thick, there are always radiative zones where moment methods are not reliable.
An alternative to moment methods are raytracing techniques. In this case, the accuracy of the general solution depends on the sampling of ray directions and the density of sampling points. Naturally, more directions and greater sampling density means the RT solver becomes more accurate but also more computationally expensive. In a shortcharacteristics raytracing scheme (e.g. Stone et al. 1992), only the neighbouring grid cells are used to interpolate the intensities. It is faster than other raybased approaches, but is also more diffusive. In a longcharacteristics scheme (e.g. Nordlund 1982; Heinemann et al. 2006), rays are traced cellbycell across long distances. This scheme provides a maximum possible accuracy at the cost of a large number of redundant calculations. Hybridcharacteristics raytracing schemes (Rijkhorst et al., 2006) are a combination of long and short characteristics methods. In addition to moment methods and raytracing, it is of course also possible to solve radiative transfer using Fourier transforms (e.g. Cen 2002) or Monte Carlo techniques (e.g. Robitaille 2011).
Although briefly described in (Nordlund et al., 2018), here we present additional details about the hybridcharacteristics raytracing RT scheme that we employ in the simulations presented here. In DISPATCH, the RT scheme uses long rays inside patches and short rays inbetween. The RT solver consists of a ray geometry component, an initialization component, and runtime schemes. In addition, since the RT depends on the values of, for example, density and temperature, it relies on the coexistence of instances of (magneto)hydrodynamical patches that can provide these quantities.
a.1 The equation of radiative transfer
The radiative transfer equation describes the change of the specific intensity during its propagation through a medium and is determined by emission and absorption processes. Solving the equation of radiative energy transport (e.g. Mihalas & Mihalas 1984) gives the net heating or cooling,
(5) 
which can then be coupled to the fluid dynamics. In Eq. 5, is a solid angle, is frequency, is the mass density, is the absorption coefficient per unit frequency, and are the specific intensity and the source function per unit area, frequency, time, and solid angle. It is advantageous to write the radiative transfer equation directly in the net difference rather than in terms of , as is typical, to circumvent numerical roundoff problems where the optical depth is very high and approaches closely to (Nordlund, 1982). At large optical depth, the difference goes as the 2 order derivative of the source function, and therefore formulations that assume constant or linear within an interval are not accurate enough. In this work, we therefore use an integral formulation which represents to 2 order. A modified Feautrier solver following Nordlund (1982) is, however, also available in the code.
a.1.1 Integral formulation
Given an optical depth scale , where is the optical distance along the propagation direction of the ray, Eq. 5 can be rewritten as
(6) 
in the direction of increasing and
(7) 
in the direction of decreasing . For simplicity, the explicit reference to the frequency has been dropped here. Equations 6 and 7 have the same form as Eq. 5, but with a term rather than . For an assumed quadratic source function , the source function derivative is a linear function, and it is therefore easy to establish an analytical solution from one point to the next:
(8)  
(9) 
where the coefficients , and are:
(10)  
and .
a.2 Ray geometry
When the RT module is initialized, it first creates a ray geometry (RG) for existing patches. The only information the RG component needs is the simulation geometry type (i.e. Cartesian/cylindrical/spherical), dimensions of the patch (number of cells per direction), the desired number of raydirections, and their angular separation.
a.2.1 Ray tracing through a patch
The ray casting starts by selecting the main direction; this is the global coordinate that forms the smallest angle with the ray. If the main direction is along, e.g., the axis, the rays are cast from the lower plane with a ray spacing in the  or direction equal to the patch cell size in that direction. The 45 angles are the exception – if they are used, the selection of the main direction is hardcoded to  and directions, thus avoiding repetition. Ray points, where the required hydrodynamic quantities are interpolated to the ray, are cellcentred along the direction of the ray. With such positioning, one needs to use 2D interpolations for a 3D space. Ideally, the ray points are coincident with patch cell centres. In this case, the required hydrodynamic quantities (e.g. density and temperature) do not need to be interpolated to the RG, and subsequently the resulting heating rates neither need to be interpolated back to the mesh, saving a significant amount of computational time. Figure 5 presents a 2D representation of the RG. Rays with three different inclinations are shown; circles denote where the hydrodynamic data is interpolated from.
To trace a ray through a patch we use a scheme that is similar to the “fast voxel traversal algorithm” (Amanatides & Woo, 1987), but is slightly simpler, and therefore a bit faster:

Consider a ray with a main direction along the axis. First, calculate the distance the ray travels to cross one cell along the main axis:
(11) where is the polar angle in the zenith direction and is the cell size in the same direction.

Next, calculate a step size in the  and directions that the ray must take in order to travel the distance :
(12) and
(13) The step size in the direction is simply .

Step back to the nearest guard cell. Rays need to know exactly where in an upstream patch or boundary the incoming heating rates should be taken from. For example, in Figure 5, if we consider the green ray, which starts in cell , the upstream point is the guard cell , as indicated by the dashed green line. This is stored as the upstream point for this ray in the forward direction.

Starting from the first point (e.g. in Figure 5), we enter a loop over the position and simply add the step size:
(14) 
After each step, check if the ray is still inside the patch (i.e. it has not crossed any of the patch boundaries). If not, continue the loop.

However, if the ray has left the patch after the last step, store the last position as the downstream point for this ray in the forward direction (or as an upstream point in the reverse direction). Depending on which wall was crossed, we either:

stop the ray if the boundary was crossed, or

restart the ray. This is done by resetting the index, which is currently outside the patch, to the lower inner value (e.g. in Figure 5, the yellow rays are restarted from the green rays while the orange rays are restarted from the red ones) and going back to step (iii).

The calculations in steps (i) and (ii) are done once per direction. This procedure requires only three floating point additions and three floating point comparisons per iteration.
a.2.2 Hierarchy
Once the casting of all the rays inside a patch is finished, we have a large stack of rays with different lengths, that have crossed different walls, and are in different directions. The results are now sorted following a particular hierarchy:

Ray directions – in principle, as long as scattering is not considered, rays in one direction do not need to know about rays in another direction, and they can be updated as independent tasks by different threads. To solve RT inside a patch, boundary conditions from the patch wall that a ray originates from are required (i.e. the upstream/downstream points that were stored during the ray tracing procedure). Note that not all slanted rays with one direction originate and terminate at the same wall (see Figure 5).

The complete set of rays is further subdivided into ray bundles; these are defined as sets of rays that originate and terminate on the same pair of patch walls. This helps avoid waiting time originating from sets of ray directions that end on more than one patch wall. Ray bundles also hold metadata about the rays it contains: their step sizes , their spatial orientation (the main direction and inclination angles), and how the rays are sorted within the bundle. Raybased RT is, in general, a repetitive task, particularly when a large number of rays is considered. It is thus highly advantageous to use schemes that maximize vectorization and therefore computational efficiency. As such, all rays should, preferably, be the same length. This is a natural feature of rays that are parallel to a patch coordinate axis.

Slanted rays, meanwhile, are further rearranged into ray packets; these are sets of rays in a ray bundle that have the same length. By organising the ray packet as a single array, we promote faster data lookup (e.g. interpolating ray coordinates to patch coordinates, origin/termination points, etc.) and vectorization.
In this work, we use 26 ray directions with 45 angular separation between rays. These rays have the advantage that they will generally terminate at locations coincident with the origins of rays in a ‘downstream’ (with respect to RT) patch. This feature is demonstrated in Figure 5, where red rays terminate at the wall and then restart at wall (yellow rays), with exactly the same position.
Finally, we note that “bundle chains” are not used in the RT in this work (in contrast to Nordlund et al. 2018). We, instead, run the RT solver in socalled “nochains” mode.
a.2.3 Initialisation and execution
After the ray geometry has been generated, the RT initialization component iterates through all of the hydrodynamical patches in the task list and creates a corresponding and overlapping RT patch. A single RT patch keeps track of all the raydirections crossing the patch. For each ray bundle, the RT patch maintains a corresponding array (with size equal to the total number of rays in this bundle) which points to the appropriate locations in the guard cells where rays must pick up the necessary upstream/downstream values.
Now that all the RT patches and rays are set up, the initialization procedure generates neighbour relations between the RT patches similarly to how it is done for hydrodynamical patches (see Nordlund et al., 2018).
Finally, RT is ready for execution. Courtesy of the connection to an underling hydrodynamical patch, the RT patch is updated immediately after its “parent” hydrodynamical patch. After the mean intensities everywhere within the patch are updated, the heating rates (Eq. 3) are calculated and stored for use by the underlying hydrodynamical patch.
a.2.4 Applying heating/cooling rates to hydrodynamics
Once the RT solver returns the frequencydependent heating rates, they are stored, via reverse volume mapping, on a mesh similar to the one belonging to the parent hydrodynamical patch. is then integrated over frequency bins and the different directions. At this point, the heating rate is ready to be applied to the hydrodynamics. Whenever the hydrodynamical patch needs the heating rates, it picks them up directly without need for further interpolation.
Appendix B Performance considerations
The current results in this work were obtained with only moderate computational effort, with each model requiring of the order of a few thousand core hours. The DISPATCH code framework (Nordlund et al., 2018) can handle arbitrary equations of state and opacities using an optimized table lookup without significantly increasing the computational cost. The RT solver is also ready for multifrequency experiments, using either representative frequencies or opacity distribution functions. Increasing the spatial resolution results, as always, in a cost increase scaling with the inverse fourth order power of the smallest resolution element. Although this can be partly mitigated by the use of local time steps (e.g. by splitting a spatial region into a number of patches each with their own time step), the scaling of the computing cost will remain essentially , where is the number of opacity bins. Taking advantage of the linear weak scaling of DISPATCH (cf. Nordlund et al., 2018, Fig. 5), it is thus possible to perform simulations with a tabular EOS and opacities, several opacity bins, and resolution increased by a factor of 8, at costs per model of the order 1020 million core hours. Alternatively one could run, using similar amounts of computing resources, a dozen or more models with a factor of four better resolution than the current endeavour.