Radiation transport methods in star formation simulations

Abstract

Radiation transport plays a crucial role in star formation models, as certain questions within this field cannot be accurately addressed without taking it into account. Given the high complexity of the interstellar medium from which stars form, numerical simulations are frequently employed to model the star formation process. This study reviews recent methods for incorporating radiation transport into star formation simulations, discussing them in terms of the used algorithms, treatment of radiation frequency dependence, the interaction of radiation with the gas, and the parallelization of methods for deployment on supercomputers. Broadly, the algorithms fall into two categories: i) moment-based methods, encompassing the flux-limited diffusion approximation, M1 closure, and variable Eddington tensor methods, and ii) methods directly solving the radiation transport equation, including forward and reverse ray tracing, characteristics-based methods, and Monte Carlo techniques. Beyond discussing advantages and disadvantages of these methods, the review also lists recent radiation hydrodynamic codes implemented the described methods.

1 Introduction

Electromagnetic radiation plays a fundamental role in the process of star formation as it regulates the dynamics, thermodynamics, and chemistry of star-forming regions. The interstellar radiation field (ISFR) heats the interstellar dust, dissociates molecules and affects the ionization balance, which is critical for the chemistry of the interstellar medium and the formation of molecular clouds. During the collapse of molecular cores, the dust absorbs the cooling radiation and once the medium becomes optically thick, the first proto-stellar core is formed. When the first stars are formed, their radiation exerts a profound influence on the surrounding gas and dust shaping the morphology of molecular clouds. This set of processes is called the radiative feedback as it regulates the rate and efficiency of star formation. Extreme ultraviolet photons from young, massive stars can ionize the surrounding gas, creating HII regions. The ionizing radiation also influences the chemistry of the gas, affecting the abundance of molecules necessary for the cooling and fragmentation of the cloud. Furthermore, radiation pressure from newly formed stars can counteract gravitational collapse, slow down the gas inflow and drive the winds and outflows.

The inclusion of radiation transport in star formation simulations is therefore desirable. However, it is also very hard due to both the high computational costs and the complexity of the processes that have to be taken into account. The radiation transport equation (RTE) in its general form readswhere the radiation field is represented by the specific intensity I_ν. It is defined through the element of the radiation energy dE_ν flowing across an area element da located at position x in time dt in the solid angle dΩ about the direction n in the frequency interval ν + dν aswhere θ is the angle that n makes with a normal to the area element da. The first term on the r.h.s of Eq. 1 represents the radiation losses due to absorption (described by the absorption opacity κ_ν,a) and scattering (described by the scattering opacity κ_ν); ρ is the mass density of the medium through which the radiation propagates. The second term on the r.h.s of Eq. 1 represents the emission of radiation describe by the emission coefficient j_ν, the last r.h.s term represents contribution of the radiation scattered from other directions. Note that in the above equations, the specific intensity is expressed in the frame of coordinate system, however, emission coefficient and opacities are defined in the comoving frame.

The structure of Eqs 1, 2 reveals that the computational costs associated with solving them directly in their complete generality are exceedingly high. Specifically, in the case of three-dimensional problems, I_ν becomes a function of three components of x, two components of n, time t, and the frequency of the radiation ν (written as a subscript, as it is often grouped into several ranges known as wavebands)—amounting to seven independent variables. This issue becomes particularly pronounced when implementing the radiation transport algorithm within a star formation simulation that demands high resolution and, consequently, a substantial number of computational points along the aforementioned variables. To address this, various computational techniques have been developed that solve only a portion of Eq. 1 (e.g., incorporating only a subset of the processes it describes) or utilize diverse approximations, including those transforming it into a different (albeit related) equation that can be solved more computationally efficiently.

This paper reviews the most common techniques to solve the radiation transport equation within (magneto-)hydrodynamic simulations of star formation. The focus of this work is primarily on continuum radiation transport, rather than line transport (atomic or molecular). There is an extensive amount of work on solving radiation transport on static grids or as tools for synthetic observations, however, covering it is beyond the scope of this work. Additional information can be found, e.g., in classical works on radiation hydrodynamics Mihalas and Mihalas (1984); Castor (2004), in a review of dust radiative transfer by Steinacker et al. (2013) and in a review of numerical methods in simulations of star formation by Teyssier and Commerçon (2019). The paper is organised as follows. Section §2 describes different way how radiation transport codes can be classified. Section §3 reviews the most common radiation transport algorithms used in star formation, consisting of §3.1 discussing the moment-based methods and §3.2 describing the methods solving Eq. 1 directly. Finally, Section §4 summarizes this review.

2 Classifications of the radiation transport treatments

There are several criteria based on which radiation transport codes can be classified.

Firstly, the method can either explicitly include the first term of Eq. 1 (partial time derivative) or assume the infinite speed of light and solve the RTE equation implicitly. The former approach, also known as the evolutionary method, is typically employed by moment-based methods. In some cases, the speed of light is artificially reduced to avoid excessively small time steps mandated by the stability conditions of an explicit solver—a technique demonstrated to be useful in other fields as well (e.g., Cohen et al., 1999). The implicit approach, alternatively termed the instantaneous method, is often utilized by codes based on the method of characteristics, ray tracing, and Monte Carlo methods.

Secondly and related to the above, the method can either directly solve Eq. 1 (albeit with some approximations) or transform it into another type of equation. The former approach includes ray tracing, long and short characteristic methods, and Monte Carlo methods. While ray tracing and long characteristic methods are among the most accurate, they are also computationally expensive. The computational cost of ray tracing typically scales with the number of discrete radiation sources, making it particularly expensive for simulations with a high number of sources. In the latter approach, Eq. 1 is commonly transformed into a hyperbolic system (M1 closure methods) or a parabolic system (flux-limited diffusion method).

Thirdly, various methods concentrate on distinct wavelengths of radiation, incorporating different physical processes related to the interaction between radiation and matter. In the simplest approach, termed monochromatic, all photons are assumed to possess a single (representative) frequency. A slightly more complex method is the grey approximation, where a specific constant spectrum—often the black body spectrum—is assumed for all radiation. In certain scenarios, it becomes necessary to categorize radiation based on its frequency into multiple groups that exhibit qualitative differences (e.g., where only some frequencies ionize a specific species). This approach is referred to as multi-group frequency treatment. Lastly, certain codes have the capability to encompass full frequency dependence by sampling the frequency with a designated number of bins.

An essential consideration in implicit radiation transport methods is whether the opacity at a specific location is depends on the radiation intensity at that point. A notable instance of this behavior is the treatment of ionizing (EUV) radiation using the on-the-spot approximation. In this approach, ionizing photons are effectively destroyed by case B recombinations, as they do not result in the replication of ionizing photons. The number of case B recombinations occurring in a given gas parcel needs to be distributed among all the rays passing through it, particularly when using ray tracing or the long characteristic method. Consequently, such methods necessitate iteration, leading to an increase in computational cost.

Finally, different methods pose varying challenges in terms of parallelization, with particular emphasis on the parallelization using domain decomposition suitable for distributed memory machines. Many hydrodynamic codes opt for this type of parallelization, as distributed memory machines offer extensive memory and processor core capabilities. Unfortunately, this parallelization scheme is also the most complex from a coding perspective. Radiation transport, being primarily a long-distance interaction problem, presents significant difficulties in parallelization. This is especially true for ray tracing, the long characteristics method, and Monte Carlo, where numerous calculations over considerable distances must be executed within a limited timeframe.

3 Radiation transport algorithms

3.1 Moment-based methods

Instead of the specific intensity I_ν, the radiation field can be described by its moments obtained by integrating I_ν over all directions. The first three moments are defined as (see, e.g., Shu, 1991):where n is a unit vector. The mean radiation energy density E_ν, the radiation flux F_ν, and the radiation pressure tensor are the zeroth, first and second moments of the radiation field, respectively. The radiation pressure is often expressed as where is the so called Eddington tensor.

Moments of Eq. 1 can be obtained by multiplying it with a certain power of n and integrating it over all directions. The scattering is assumed to be isotropic. Then, the first two moment equations areandIn the above equations, the opacity and the emission coefficient are expressed in the frame of reference co-moving with the fluid (where they are typically isotropic), however, E_ν, F_ν and are expressed in the coordinate system of the simulation (the laboratory frame). This problem is usually solved by either transforming Eqs 4, 5 into the co-moving frame (CMF) or by transforming κ_ν,a and j_ν into the laboratory frame using the first order expansion (mixed frame formulation). The radiation moment equations in the co-moving frame have a form (see Mihalas and Mihalas, 1984)andwhere u is the velocity of the gas. For details on the mixed frame formulation see e.g., Mihalas and Klein (1982). The advantage of the mixed frame approach is that the radiation transfer equations are significantly simpler and that it conserves the total energy, however, it is not suitable for computing for instance line opacities in the supersonic flow due to their rapid changes (Moens et al., 2022).

3.1.1 Flux limited diffusion (FLD)

The simplest moment-based method was derived as an extension of the radiation transport in a diffusion limit (Minerbo, 1978; Levermore and Pomraning, 1981). If the photon free mean path is small and the optical depth is high, the radiation is isotropic and the Eddington tensor has form where is the identity tensor. Assuming the radiation field is stationary, i.e., neglecting the first two terms in Eq. 7, and considering , the radiation flux is related to the radiation energy density as

and the radiation transport problem reduces to the Fick’s diffusion law.

Relation 8 is not valid if the optical depth is low, because F_ν can exceed the physical limit F_max = cE_ν there. Therefore, the FLD method introduces the so called flux limiterλ and modifies Eq. 8 in the following wayThe flux limiter λ(R) is defined as a function of the normalized radiation energy density gradient

in such a way that it ensures Eq. 9 reduces to the Fick’s law in the optically thick regime (i.e., λ = 1/3) and to the maximum allowed flux, F_max, in the optically thin regime (i.e., λ = 1/R). Various flux limiters fulfilling the above requirements have been suggested in the literature (see, e.g., Minerbo, 1978; Levermore and Pomraning, 1981; Levermore, 1984; Turner and Stone, 2001). As an example we give the flux limiter of Levermore and Pomraning (1981) in the form

and the corresponding approximate Eddington tensor

The FLD method then solves only the 0-th moment equation in the formwhich is, from a mathematical point of view, a parabolic equation, and hence the standard parabolic solvers can be used.

The above equations are monochromatic, i.e., they are written for a single wavelength of the radiation. However, it is often practical to use the FLD method (and other methods as well) for radiation with a certain range of wavelengths. Then we speak about the grey radiation transport if there is a single range covering the broad wavelength range important for a given problem, or about the multi-group/multi-band radiation transport if there are several wavelength ranges. Additionally, the grey methods often use the Planck black body law, B_ν(T), as the source function, and then the 0-th moment equation has formwhere E_r ≡∫E_νdν is the bolometric radiation energy density, T is the gas/dust temperature, a_r = 4σ_B/c is the radiation constant, σ_B is the Stefan-Boltzmann constant, and κ_P and κ_R are Planck and Rosseland mean opacities defined asand

RHD codes using the FLD method are discussed below and summarized in Table 1 giving a reference to the code paper, code name, whether it uses co-moving or mixed frame of reference, whether it is parallel with the domain decomposition and some additional information. The first usage of the FLD method is probably by Alme and Wilson (1973) for the problem of the X-ray emission resulting from the accretion of material onto a neutron star. In the field related to star formation, the grey FLD method was used by Kley (1989); Bodenheimer et al. (1990); Klahr et al. (1999); Boss (2001) for simulations of accretion discs around young low mass stars. It yielded reasonably correct temperatures in the optically thick and optically thing regions, however, it showed deviations from the correct solution in the transition regions (Boley et al., 2007). The treatment of radiation in such models has been significantly improved by Kuiper et al. (2010) who developed a hybrid scheme combining the grey FLD method with the frequency dependent ray tracing, and implemented it to the MHD code PLUTO. They at first calculate the stellar radiation up to the first absorption by a simple ray tracing algorithm and subsequently, they use FLD to calculate the radiation re-emitted by the dust. This approach turned out to be very successful and many other authors implemented this or similar schemes into their codes. Bitsch et al. (2013) used it with the NIRVANA code to study the influence of opacity and stellar irradiation on the disc structure and on the migration of planets. Flock et al. (2013) and Kolb et al. (2013) created two independent implementations of this hybrid method for the PLUTO code. Klassen et al. (2014) generalized this method for multiple sources and arbitrary geometry wrote it as a module for the MHD code FLASH. Similarly, Ramsey and Dullemond (2015) implemented a generalized version of this method for the AZEUS code—a ZEUS family AMR code with a fully staggered mesh.

The simplicity of the FLD method allowed it to be implemented into many general purpose (magneto-) hydrodynamic codes used in astrophysics. The widely used ZEUS2D code was originally released with the radiation module based on the tensor variable Eddington factor, however, the FLD method was written for it by Turner and Stone (2001). The ORION code used for many simulations of star formation, including the feedback from massive stars, was developed with FLD method being part of it (Krumholz et al., 2007). Gittings et al. (2008) developed the AMR code RAGE including the grey FLD, for which they introduce a new technique of solving the diffusion and material energy equations, allowing larger time steps and more robust behavior. Commerçon et al. (2011) implemented FLD into the AMR code RAMSES, and Commerçon et al. (2014) added adaptive time stepping and the implicit time integration. The block AMR code CRASH (van der Holst et al., 2011) was developed with grey FLD where electrons and ions are allowed to have different temperatures. The grey FLD method in its mixed-frame formulation is also part of the AMR code CASTRO (Zhang et al., 2011), one of the first 3D RHD codes using the unsplit PPM solver. TRHD (Sijoy and Chaturvedi, 2015) is two-dimensional unstructured mesh Lagrangian code using the grey FLD with three temperature approach (ions, electrons and radiation are allowed to have different temperatures). The FLASH code (Fryxell et al., 2000), a computational framework consisting of modules implementing various physical processes including hydrodynamics and radiation transport, was provided by the FLD module implemented by Chatzopoulos and Weide (2019) and tested on 1D simulations of supernova Ia explosions. The general fluid dynamics code GIZMO (Hopkins and Grudić, 2019) includes several radiation transport solvers including the one based on FLD. Recently, the grey FLD method was implemented into the finite-volume MHD code MPI-AMRVAC (Moens et al., 2022).

The FLD method was implemented also into the Smoothed Particle Hydrodynamics (SPH) codes. Basic methods for calculating radiative heat diffusion were described in seminal works on SPH by Lucy (1977) and Brookshaw (1985). The full grey FLD method was implemented into the SPH by Whitehouse and Bate (2004) and Whitehouse et al. (2005). Later, this code originating from code SPHNG (Benz, 1990) was used by Price and Bate (2009); Bate and Keto (2015); Bate (2019) and many other works for simulating star forming clouds, where the FLD method can be efficiently used to follow collapsing pre-stellar cores with high optical depths. Mayer et al. (2007) used implemented FLD into SPH code GASOLINE and used it to study the gravitational instability in protoplanetary discs. Fryer et al. (2006) developed code SNSPH including the FLD radiation transport and used it for simulating supernova explosions, where FLD can also be used efficiently due to high optical depths. Another implementations of FLD into SPH codes are by Viau et al. (2006) and recently by Bassett et al. (2021) (code SPHERAL).

3.1.2 M1 closure

A more complex and accurate two-moments method is the M1 closure (Dubroca and Feugeas, 1999). It assumes that the Eddington tensor is axially symmetric around the flux vector and hence can be expressed aswhere χ is a scalar function called the Eddington factor. The M1 closure method uses one of its simplest forms discussed in Levermore (1984), corresponding to the Lorentz boosted isotropic distribution of the radiation,where f = |F_ν|/cE_ν is the reduced flux. Then, the method solves Eqs 6, 7 or their equivalent in the mixed frame approach.

The advantage of the M1 closure is that it captures correctly the two extremes of the radiation field: i) the free streaming (or optically thin) limit with |f| = 1 and χ = 1, and ii) the diffusion (or optically thick) limit with |f| = 0 and χ = 1/3. The ability to maintain the original direction of the free streaming radiation is a pronounced improvement in comparison to the FLD method leading to the artificial diffusion of the free streaming radiation. On the other hand, the assumed simple form of the radiation pressure tensor cannot correctly describe more complex radiation fields as for instance the radiation from multiple sources in which case the nonphysical interaction would occur.

A specific feature of the M1 closure making is suitable for the implementation into (magneto-) hydrodynamic codes is that the system of Eqs 6, 7 is hyperbolic and hence it can be solved with the same algorithms as the equations of the fluid dynamics. Table 2 gives a list of star formation and related fields codes implementing the M1 closure method. The first such code was HERACLES (González et al., 2007) designed to study radiative shocks (and compare the results with laboratory experiments), jets of young stars, formation of pre-stellar cores and protoplanetary disks. This method was implemented into the cosmic reionisation code ATON (Aubert and Teyssier, 2008), and later ported to the general RHD code RAMSES (Rosdahl et al., 2013; Rosdahl and Teyssier, 2015). Sądowski et al. (2013) included the M1 closure method into the general relativistic code KORAL formulating the RHD equations in the covariant form. Skinner and Ostriker (2013) describe the M1 closure module for the ATHENA code, in which they use the reduced speed of light approximation to improve the performance. The general fluid dynamics code GIZMO includes the M1 closure as one of its radiation transport solvers. Kannan et al. (2019) describe the radiation transport module for the unstructured moving mesh code AREPO, based on a new higher order implementation of the M1 closure method. This method is also included into code FORNAX (Skinner et al., 2019) primarily intended for core collapse supernova simulations, general fuild dynamics code PLUTO (Melon Fuksman and Mignone, 2019), and hybrid OpenMP/MPI/GPU parallel code ARK-RT (Bloch et al., 2021). Chan et al. (2021) present module SPH-M1RT for task-based parallel SPH code SWIFT.

3.1.3 Variable Eddington Tensor methods

Another approach to solve moment Eqs 4, 5, as opposed to the closure assuming an approximate form of the Eddington tensor, is to calculate the Eddington tensor self consistently by directly solving the radiation transport Eq. 1. Methods using it (Table 3) are referred to as Variable Eddington Tensor (VET) methods, and the Eddington tensor, , is calculated by one of the methods described in §3.2 (ray tracing, long or short characteristics, or Monte Carlo). Subsequently, is combined with the radiation moment equations to calculate the radiation energies and fluxes.

The VET method was developed for 1D calculations of plane-parallel stellar atmospheres by Auer and Mihalas (1970), however, the idea of iterating the ratio of second and zeroth radiation moments dates back to Eddington (1926). The first multi-dimensional astrophysical code using this method was ZEUS2D (Stone et al., 1992), where was calculated by the short characteristics. The algorithm was later updated in code ZEUS-MP (Hayes and Norman, 2003) by adding more methods to calculate and by separating the radiation and gas energy densities. The 1D RHD code TITAN (Gehmeyr and Mihalas, 1994) uses VET and calculates by ray tracing. Gnedin and Abel (2001) developed an algorithm based on the Optically Thin Variable Eddington Tensor (OTVET) approximation where is calculated in the optically thin regime. This is done by integrating contributions of radiation sources with attenuation factor 1/r² over the computational domain taking advantage of the solver used for the self-gravity. They implemented this algorithm into the Soften Lagrangian Hydrodynamic code (SLH) following all physical quantities on a moving deformed mesh. This algorithm was later implemented by Petkova and Springel (2009) into the SPH code GADGET. Sekora and Stone (2010) developed a hybrid Godunov method based on VET with user defined and implemented it into 1D RHD code NIKE. This method was later implemented into code ATHENA by Jiang et al. (2012) where is calculated by the short characteristics method described by Davis et al. (2012). The general relativistic radiation magnetohydrodynamic (GR-RMHD) code INAZUMA (Asahina et al., 2020) also uses the VET based method to treat the radiation, and it calculates integrating the local radiation intensity over the angular directions. Recently, the RHD code QUOKKA (Wibking and Krumholz, 2022) implements the general two-moment radiation transport method. The code uses the block structured adaptive mesh refinement handled by the AmREX library and it is parallelized for both MPI and usage on graphic processing units (GPUs) using the CUDA library. Its radiation solver is written to use a general Eddington tensor, the default option is the M1 closure.

3.2 Methods solving RTE directly

Table 4 lists codes that employ the direct solution of the RTE: forward and reverse ray tracing, methods of characteristics, and Monte Carlo methods.

3.2.1 Ray tracing

The most straightforward method to solve the radiation transport Eq. 1 is to follow the radiation from its sources along the lines of its propagation, a technique called ray tracing. It is widely used in many other fields, in particular in the computer image generation. The original idea of ray tracing dates back to 16th century to Albrecht Dürer who described multiple techniques for projecting 3D scenes onto an image plane (Hofmann, 1990).

As previously mentioned, the primary challenge in directly solving Eq. 1 is its high dimensionality, resulting in high computational costs. Therefore, a key characteristic of a ray tracing algorithm is a set of approximations designed to reduce the complexity of problem. A crucial aspect is the method by which the radiation field is discretized in both space and its propagation directions. The simplest approach, employed by many early and recent codes, involves calculating the radiation on the same grid as the gas dynamics. This confines the direction of radiation propagation to the coordinate axes. Examples of this approach include Wolfire and Cassinelli (1986) in a 1D spherical grid, Murray et al. (1994); Garcia-Segura and Franco (1996) in a 2D spherical grid, and numerous other works, the review of which exceeds the scope of this study. This method is particularly well-suited for treating radiation emanating from a single source (e.g., a central star) located at the center of the coordinate system. It can be effectively combined with other methods that describe the diffuse radiation component (e.g., FLD, as mentioned in Section 3.1.1; Table 1).

In many cases, confining the radiation propagation solely along the hydrocode grid axes is overly restrictive. Therefore, general radiation transport codes typically employ an independent method to define rays, making ray tracing algorithms particularly suitable for implementation into unstructured grid or smooth particle hydrodynamic codes. One approach is to cast rays from a source in spherical coordinates, associating each ray with a spherical segment of the same size (Abel et al., 1999; Razoumov and Scott, 1999). This method was later refined in the HEAlPix library (Górski et al., 2005), which has become a de facto standard. HEALPix facilitates the tessellation of the sphere surface into equal area regions in a hierarchical manner, a feature utilized by adaptive ray tracing algorithms that split rays in regions requiring finer angular resolution (Abel and Wandelt, 2002). This approach has proven highly successful and has been employed by various authors, including Pawlik and Schaye (2008) (TRAPHIC code linked to SPH code GADGET2), Bisbas et al. (2009) (SPH code SEREN), Wise and Abel (2011) (MORAY code linked to AMR code ENZO), Baczynski et al. (2015) (FERVENT code within the AMR code FLASH), and others.

Depending on the number of cast rays, ray tracing algorithms can achieve high accuracy, although with correspondingly higher computational costs that typically scale with the number of radiation sources. As a result, they are occasionally employed in conjunction with other methods where ray tracing specifically handles radiation from a relatively low number of discrete sources, such as massive stars. An example is the HARM2 code (Rosen et al., 2017), which calculates ionizing stellar radiation using adaptive ray tracing and the diffuse component using the FLD method in the ORION code. Similarly, the TORUS-3DPDR code (Bisbas et al., 2015a) utilizes adaptive ray tracing for FUV radiation and the Monte Carlo method for the ionizing and diffuse components.

A slightly different approach to spatially adaptive ray tracing was used by Dale et al. (2007) who integrated the RTE on the line-of-sight between the source and the target particle using evaluation points defined by SPH particles close to the line-of-sight selected using method by Kessel-Deynet and Burkert (2000).

Ritzerveld and Icke (2006) developed algorithm SimpleX to solve Eq. 1 on an unstructured grid calculated from the properties of the medium in which the radiation propagates, using the Delaunay triangulation (Delaunay, 1934). The computational costs of this method do not depend on the number of sources, similarly to moment-based methods or some reverse ray tracing schemes. It is typically less diffusive that moment based methods and it is more suitable for parallelisation with domain decomposition than most of the ray tracing methods. It is particularly suitable for being used with SPH or unstructured grid codes, Clementel et al. (2014) implemented it into an SPH code. The SimpleX algorithm was improved and parallelized by Paardekooper et al. (2010) and implemented as code SPRAI integrated into code AREPO by Jaura et al. (2018).

3.2.2 Characteristics-based methods

A specific way of ray tracing are methods of characteristics, which solve the RTE along predefined rays (characteristics) covering the entire computational domain rather than casting rays from discrete sources. In this manner, characteristics-based methods can describe the diffuse radiation component. They can be categorized into long, short and hybrid characteristics based on the end points of the rays. Long characteristics extend across the entire computational domain, with their separation comparable to the size of the computational elements and defined in specific directions. While long characteristics are among the most accurate methods, they are also the most expensive, depending on the number of rays. Another disadvantage, applicable to some extent to all ray tracing methods, is that they are difficult to get parallelized within codes that use domain decomposition on distributed memory machines. This challenge arises because rays crossing multiple domains typically represent a substantial amount of information that needs to be communicated over the (relatively) slow network. An example of a code utilizing long characteristics is LAMPRAY (Frostholm et al., 2018), implemented into the RAMSES code. LAMPRAY calculates both diffuse and discrete sources radiation on an adaptive oct-tree mesh, allowing the definition of several frequency bins for both ionizing and non-ionizing radiation. It is coupled with the non-equilibrium chemistry code KROME (Grassi et al., 2014), providing the chemical networks to evolve species abundances.

The short characteristics methods create rays only among the neighbour cells (or generally computational elements) and interpolate the radiation intensity at the cell borders. This makes this class of methods more diffusive, but also computationally cheaper and much easier to parallelize with domain decomposition codes. Short characteristics are used, e.g., by code C2RAY (Mellema et al., 2006), code ATHENA (Davis et al., 2012), or code PION (Mackey et al., 2021). They have been also used together with the moment-based VET method to calculate the local radiation pressure tensor (see §3.1.3).

Rijkhorst et al. (2006) developed the hybrid characteristics method, combining the advantages of both long and short characteristics, particularly in reducing the necessary parallel communication. This method is implemented within the block-based daptive AMR code FLASH. It divides the radiation field into two components: i) the local component within the blocks, calculated using the long characteristics method, and ii) the global component, transporting radiation among blocks (and parallel sub-domains) using interpolation similar to the short characteristics method. This algorithm was further enhanced by Peters et al. (2010) and later re-implemented by Buntemeyer et al. (2016) to also handle diffuse radiation and operate on GPU-accelerated architectures.

Jiang (2021) developed an implicit radiation transport algorithm for MHD code ATHENA++. It represents the radiation by its specific intensities along discrete rays, and evolves them using a conservative finite volume method. It avoids using a closure as in moment-based methods and the related artificial diffusion. On the other hand, it solves time-dependent radiation transport equation, contrary to majority of ray or characteristics based methods. The original implementation was grey, however, Jiang (2022) extended it to include the frequency dependence using the multi-group approach.

3.2.3 Reverse ray tracing

An alternative approach is the so-called reverse ray tracing. In contrast to the “standard” forward ray tracing, which casts rays from the source outward, reverse ray tracing schemes cast rays from the target point where the radiation interacts with the gas (e.g., grid cell or SPH particle) and transport the radiation inward to the target point. This approach offers several advantages, including the highest density of rays (i.e., the highest spatial resolution) near the location where radiation interacts with the gas. Another advantage is its suitability for handling external diffuse radiation. Hasegawa and Umemura (2010) implemented this algorithm into the SPH code START, incorporating tree-based acceleration to group distant sources into oct-tree nodes. A grid-based version of this algorithm was developed by Okamoto et al. (2012) for the code ARGOT. Altay and Theuns (2013) created the reverse ray tracing code URCHIN and integrated it into an SPH code focused on cosmological simulations. Another implementation of this approach is TREVR for the SPH code GASOLINE by Grond et al. (2019), who introduced adaptive refinement along rays and directionally dependent absorption coefficients for tree nodes. Wünsch et al. (2021) developed the code TREERAY, implementing tree-accelerated reverse ray tracing for the Adaptive Mesh Refinement (AMR) code FLASH. It utilizes a fast parallel tree solver described in Wünsch et al. (2018) and is coupled with the chemical network developed by Glover and Mac Low (2007). Later, it was extended by Klepitko et al. (2023), who added diffuse grey infrared radiation and related radiation pressure, and by Gaches et al. (2023), who added frequency-dependent X-ray radiation. TREERAY has been used in the SILCC project (Walch et al., 2015, and follow-up papers) using RMHD simulations to model the full cycle of molecular clouds in a part of a galactic disc.

3.2.4 Monte Carlo methods

Monte Carlo Radiation Transport (MCRT) algorithms belong to a broad category of statistical sampling methods that share this name. They simulate the movement of individual photons through a medium by stochastically sampling random events, including interactions with matter. Given that tracking individual photons is prohibitively expensive, MCRT algorithms introduce “photon packets” that carry information about the position, direction, frequency distribution, etc., of the group of photons they represent. This approach closely mirrors the reality of propagating radiation, making it highly versatile, and facilitating the inclusion of various complex processes such as scattering.

However, the statistical nature of this method introduces Poisson noise, and maintaining it below a required level often necessitates a very large number of photon packets. Consequently, MCRT codes have historically mostly been used for static problems (e.g., the code MOCASSIN by Ercolano et al., 2003, the code HYPERION by Robitaille, 2011, or the code RADMC3D by Dullemond et al., 2012). Due to the potential long distances photon packets can travel, Monte Carlo methods are relatively challenging to parallelize using domain decomposition. On the other hand, thread-based parallelization on shared memory architectures is typically straightforward.

As this method does not rely on any grid, it is natural to use it with SPH or unstructured grid codes. An algorithm to incorporate MCRT into SPH codes was proposed by Lucy (1999). A Monte Carlo-based ray-tracing algorithm was developed by Altay et al. (2008) and implemented into the publicly available code SPHRAY, primarily focused on cosmological simulations. Another implementation of MCRT into an SPH code was done by Nayakshin et al. (2009), who used it to calculate the radiation pressure force. Roth and Kasen (2015) developed a 1D MCRT code and coupled it with both Eulerian and Lagrangian hydrodynamic codes. The first star formation simulations with MCRT were performed by Lomax and Whitworth (2016), who developed the code SPAMCART coupled with SPH hydrodynamics. They used it to calculate synthetic intensity maps and spectra of an embedded protostellar multiple system. Vandenbroucke and Wood (2018) developed the radiation hydrodynamic code CMACIONIZE, which works on an unstructured moving grid and is based on the Monte Carlo method calculating ν-dependent transport of ionizing photons. The code is parallelized using a task-based scheme and has an option to work in a distributed memory configuration. Code TORUS, developed by Harries et al. (2019), is an MCRT that includes native AMR grid-based hydrodynamics and can also be coupled with the SPH code SPHNG. It implements numerous microphysics modules, including atomic and molecular line transport in moving media, dust radiative equilibrium, photoionization equilibrium, and time-dependent radiative transfer, and has been used, e.g., for modelling massive stars feedback in turbulent clouds (Ali et al., 2018). It is also designed for postprocessing the output of several other (magneto-)hydrodynamic codes. Recently, MCRT has also been implemented into the unstructured moving mesh code AREPO by Smith et al. (2020).

3.3 Other approximations

Many studies conducting hydrodynamic simulations of star formation further simplify the radiation transport problem. For example, Vázquez-Semadeni et al. (2017) include stellar ionizing radiation feedback by calculating the Strömgren sphere radius from the EUV photon production rate of a star and the characteristic density obtained as the geometric mean of the density at the star and at the target cell. Subsequently, they set the gas temperature within the sphere to 10⁴ K.

A specific case involves non-ionizing diffuse ambient interstellar radiation fields, crucial for the physics of molecular clouds as they provide heating and affect molecular gas chemistry. In molecular clouds, this radiation intensity strongly depends on shielding, i.e., the optical depth between a given location and the edge of the cloud. This consideration led Clark et al. (2012) to develop the TREECOL algorithm, implemented into the code GADGET2. It efficiently calculates optical depths in all directions (discretized by HealPix) for each SPH particle using a tree, enabling the determination of heating and dissociation rates. This algorithm was later implemented by Valdivia and Hennebelle (2014) into RAMSES and by Wünsch et al. (2018) into FLASH.

A different approach to solving a very similar problem was suggested by Stamatellos et al. (2007), who estimate the mean optical depth from local density, temperature, and gravitational potential using a set of precalculated models of collapsing clouds. This method provides reasonably good results for almost zero computational costs.

The simplest treatment of radiation is to modify the gas equation of state to calculate the equilibrium gas temperature from local gas density, thus accounting for the most probable result of the heating-cooling balance. A popular choice is the barotropic equation of state (e.g., Bonnell, 1994; Whitworth et al., 1995; Bate, 1998), a simple analytical formula designed to mimic the thermodynamics of spherically symmetric collapse of a single, isolated protostar.

4 Summary

This work reviews methods to solve radiation transport equation (RTE) in simulations of star formation and related fields. According to the algorithm, they can be roughly divided into two groups: moment-based methods and methods solving the RTE directly.

Moment-based methods approximate the radiation specific intensity by integrating it over all directions and expressing it in terms of several moments. The order of the moments in the sequence closure determines the level of approximation, leading to artificial diffusion of the radiation, unphysical self-interaction, and other artificial effects. These methods typically include the time-dependent term of the RTE and sometimes artificially decrease the speed of light to reduce computational costs. Their computational expenses are independent of the number of sources, making them well-suited for calculating diffuse radiation. Frequency-wise, moment-based methods often use the grey approximation or, occasionally, the multi-group frequency approach. They are sometimes combined with ray tracing, which calculates radiation from one or several discrete sources, while the moment-based method addresses the diffuse radiation component. The most commonly used methods of this type include flux-limited diffusion, the M1 closure, and the variable Eddington tensor method.

Methods that directly solve the RTE integrate it along specific lines, commonly referred to as rays. These methods include forward and reverse ray tracing, long, hybrid, and short characteristics, as well as Monte Carlo methods. Ray tracing and characteristics-based methods typically do not include the time-dependent RTE term, utilizing the infinite speed of light approximation; however, there are exceptions. Monte Carlo methods are usually time-dependent. Depending on angular discretization, ray tracing and long/hybrid characteristics can achieve high accuracy but at the expense of high computational costs. Short characteristics methods tend to be more diffusive. While Monte Carlo methods can be highly accurate, achieving this requires a substantial number of photon packets, leading to increased computational costs. Monte Carlo methods are also versatile, suitable, e.g., for treating radiation scattering. Computational costs for most of these methods scale with the number of sources, with exceptions like reverse ray tracing methods that group distant sources and methods combining ray tracing with advection schemes. Various treatments of radiation frequency dependence are employed; generally, methods designed for a small number of sources can incorporate detailed frequency dependence. Forward ray tracing and long characteristics methods are typically challenging to parallelize using domain decomposition, whereas short characteristics are relatively straightforward. Hybrid characteristics, reverse ray tracing, and Monte Carlo methods fall in the middle ground regarding parallelization difficulty.

Radiation transport has become a well-established component of star formation simulations, with many methods and codes available, continuously advancing in both maturity and complexity. Consequently, it becomes imperative to not only validate these codes against standard tests but also to benchmark them against each other and, ideally, against realistic problems inherent in the realm of star formation. In this regard, initiatives such as the code comparison projects, such as those conducted by Iliev et al. (2006), Iliev et al. (2009) for cosmological radiation transport codes and the STARBENCH project led by Bisbas et al. (2015b) that compares radiation transport codes for star formation, prove to be very valuable. Moreover, there is a growing need for more comparison projects to enhance our understanding and confidence in the capabilities of these codes.

References

• 1

  AbelT.NormanM. L.MadauP. (1999). Photon-conserving radiative transfer around point sources in multidimensional numerical cosmology. ApJ523, 66–71. 10.1086/307739

  • CrossRef
  • Google Scholar

• 2

  AbelT.WandeltB. D. (2002). Adaptive ray tracing for radiative transfer around point sources. MNRAS330, L53–L56. 10.1046/j.1365-8711.2002.05206.x

  • CrossRef
  • Google Scholar

• 3

  AliA.HarriesT. J.DouglasT. A. (2018). Modelling massive star feedback with Monte Carlo radiation hydrodynamics: photoionization and radiation pressure in a turbulent cloud. MNRAS477, 5422–5436. 10.1093/mnras/sty1001

  • CrossRef
  • Google Scholar

• 4

  AlmeM. L.WilsonJ. R. (1973). X-ray emission from a neutron star accreting material. ApJ186, 1015–1026. 10.1086/152566

  • CrossRef
  • Google Scholar

• 5

  AltayG.CroftR. A. C.PelupessyI. (2008). SPHRAY: a smoothed particle hydrodynamics ray tracer for radiative transfer. MNRAS386, 1931–1946. 10.1111/j.1365-2966.2008.13212.x

  • CrossRef
  • Google Scholar

• 6

  AltayG.TheunsT. (2013). URCHIN: a reverse ray tracer for astrophysical applications. MNRAS434, 748–764. 10.1093/mnras/stt1067

  • CrossRef
  • Google Scholar

• 7

  AsahinaY.TakahashiH. R.OhsugaK. (2020). A numerical scheme for general relativistic radiation magnetohydrodynamics based on solving a grid-based Boltzmann equation. ApJ901, 96. 10.3847/1538-4357/abaf51

  • CrossRef
  • Google Scholar

• 8

  AubertD.TeyssierR. (2008). A radiative transfer scheme for cosmological reionization based on a local Eddington tensor. MNRAS387, 295–307. 10.1111/j.1365-2966.2008.13223.x

  • CrossRef
  • Google Scholar

• 9

  AuerL. H.MihalasD. (1970). On the use of variable Eddington factors in non-LTE stellar atmospheres computations. MNRAS149, 65–74. 10.1093/mnras/149.1.65

  • CrossRef
  • Google Scholar

• 10

  BaczynskiC.GloverS. C. O.KlessenR. S. (2015). Fervent: chemistry-coupled, ionizing and non-ionizing radiative feedback in hydrodynamical simulations. MNRAS454, 380–411. 10.1093/mnras/stv1906

  • CrossRef
  • Google Scholar

• 11

  BassettB. R.OwenJ. M.BrunnerT. A. (2021). Efficient smoothed particle radiation hydrodynamics I: thermal radiative transfer. J. Comput. Phys.429, 109996. 10.1016/j.jcp.2020.109996

  • CrossRef
  • Google Scholar

• 12

  BateM. R. (1998). Collapse of a molecular cloud core to stellar densities: the first three-dimensional calculations. ApJ Let.508, L95–L98. 10.1086/311719

  • CrossRef
  • Google Scholar

• 13

  BateM. R. (2019). The statistical properties of stars and their dependence on metallicity. MNRAS484, 2341–2361. 10.1093/mnras/stz103

  • CrossRef
  • Google Scholar

• 14

  BateM. R.KetoE. R. (2015). Combining radiative transfer and diffuse interstellar medium physics to model star formation. MNRAS449, 2643–2667. 10.1093/mnras/stv451

  • CrossRef
  • Google Scholar

• 15

  BenzW. (1990). “Smooth particle hydrodynamics - a review,” in Numerical modelling of nonlinear stellar pulsations problems and prospects. Editor BuchlerJ. R. (Netherlands: Springer), 269.

  • Google Scholar

• 16

  BisbasT. G.HaworthT. J.BarlowM. J.VitiS.HarriesT. J.BellT.et al (2015a). TORUS-3DPDR: a self-consistent code treating three-dimensional photoionization and photodissociation regions. MNRAS454, 2828–2843. 10.1093/mnras/stv2156

  • CrossRef
  • Google Scholar

• 17

  BisbasT. G.HaworthT. J.WilliamsR. J. R.MackeyJ.TremblinP.RagaA. C.et al (2015b). STARBENCH: the D-type expansion of an H II region. MNRAS453, 1324–1343. 10.1093/mnras/stv1659

  • CrossRef
  • Google Scholar

• 18

  BisbasT. G.WünschR.WhitworthA. P.HubberD. A. (2009). Smoothed particle hydrodynamics simulations of expanding H II regions. I. Numerical method and applications. A&A497, 649–659. 10.1051/0004-6361/200811522

  • CrossRef
  • Google Scholar

• 19

  BitschB.CridaA.MorbidelliA.KleyW.Dobbs-DixonI. (2013). Stellar irradiated discs and implications on migration of embedded planets. I. Equilibrium discs. A&A549, A124. 10.1051/0004-6361/201220159

  • CrossRef
  • Google Scholar

• 20

  BlochH.TremblinP.GonzálezM.PadioleauT.AuditE. (2021). A high-performance and portable asymptotic preserving radiation hydrodynamics code with the M₁ model. A&A646, A123. 10.1051/0004-6361/202038579

  • CrossRef
  • Google Scholar

• 21

  BodenheimerP.YorkeH. W.RozyczkaM.TohlineJ. E. (1990). The Formation phase of the solar nebula. ApJ355, 651. 10.1086/168798

  • CrossRef
  • Google Scholar

• 22

  BoleyA. C.DurisenR. H.NordlundÅ.LordJ. (2007). Three-dimensional radiative hydrodynamics for disk stability simulations: a proposed testing standard and new results. ApJ665, 1254–1267. 10.1086/519767

  • CrossRef
  • Google Scholar

• 23

  BonnellI. A. (1994). A new binary formation mechanism. MNRAS269, 837–848. 10.1093/mnras/269.3.837

  • CrossRef
  • Google Scholar

• 24

  BossA. P. (2001). Gas giant protoplanet formation: disk instability models with thermodynamics and radiative transfer. ApJ563, 367–373. 10.1086/323694

  • CrossRef
  • Google Scholar

• 25

  BrookshawL. (1985). A method of calculating radiative heat diffusion in particle simulations. PASA6, 207–210. 10.1017/S1323358000018117

  • CrossRef
  • Google Scholar

• 26

  BuntemeyerL.BanerjeeR.PetersT.KlassenM.PudritzR. E. (2016). Radiation hydrodynamics using characteristics on adaptive decomposed domains for massively parallel star formation simulations. New Astron.43, 49–69. 10.1016/j.newast.2015.07.002

  • CrossRef
  • Google Scholar

• 27

  CastorJ. I. (2004). Radiation hydrodynamics. Cambridge University Press.

  • Google Scholar

• 28

  ChanT. K.TheunsT.BowerR.FrenkC. (2021). Smoothed particle radiation hydrodynamics: two-moment method with local Eddington tensor closure. MNRAS505, 5784–5814. 10.1093/mnras/stab1686

  • CrossRef
  • Google Scholar

• 29

  ChatzopoulosE.WeideK. (2019). Gray radiation hydrodynamics with the FLASH code for astrophysical applications. ApJ876, 148. 10.3847/1538-4357/ab18f9

  • CrossRef
  • Google Scholar

• 30

  ClarkP. C.GloverS. C. O.KlessenR. S. (2012). TreeCol: a novel approach to estimating column densities in astrophysical simulations. MNRAS420, 745–756. 10.1111/j.1365-2966.2011.20087.x

  • CrossRef
  • Google Scholar

• 31

  ClementelN.MaduraT. I.KruipC. J. H.IckeV.GullT. R. (2014). 3D radiative transfer in η Carinae: application of the SIMPLEX algorithm to 3D SPH simulations of binary colliding winds. MNRAS443, 2475–2491. 10.1093/mnras/stu1287

  • CrossRef
  • Google Scholar

• 32

  CohenJ.PratchettT.StewardI. (1999). The science of discworld. London: Ebury Press.

  • Google Scholar

• 33

  CommerçonB.DeboutV.TeyssierR. (2014). A fast, robust, and simple implicit method for adaptive time-stepping on adaptive mesh-refinement grids. A&A563, A11. 10.1051/0004-6361/201322858

  • CrossRef
  • Google Scholar

• 34

  CommerçonB.TeyssierR.AuditE.HennebelleP.ChabrierG. (2011). Radiation hydrodynamics with adaptive mesh refinement and application to prestellar core collapse. I. Methods. A&A529, A35. 10.1051/0004-6361/201015880

  • CrossRef
  • Google Scholar

• 35

  DaleJ. E.ErcolanoB.ClarkeC. J. (2007). A new algorithm for modelling photoionizing radiation in smoothed particle hydrodynamics. MNRAS382, 1759–1767. 10.1111/j.1365-2966.2007.12486.x

  • CrossRef
  • Google Scholar

• 36

  DavisS. W.StoneJ. M.JiangY.-F. (2012). A radiation transfer solver for athena using short characteristics. ApJS199, 9. 10.1088/0067-0049/199/1/9

  • CrossRef
  • Google Scholar

• 37

  DelaunayB. (1934). Sur la sphère vide. A la mémoire de Georges Voronoï. Bull. de l’Académie des Sci. de l’URSS. Cl. des Sci. mathématiques naturelles6, 793–800.

  • Google Scholar

• 38

  DubrocaB.FeugeasJ. (1999). Etude théorique et numérique d’une hiérarchie de modèles aux moments pour le transfert radiatif. Acad. Des. Sci. Paris Comptes Rendus Ser. Sci. Math.329, 915–920. 10.1016/S0764-4442(00)87499-6

  • CrossRef
  • Google Scholar

• 39

  DullemondC. P.JuhaszA.PohlA.SereshtiF.ShettyR.PetersT.et al (2012). RADMC-3D: a multi-purpose radiative transfer tool. N/A ascl (Astrophysics Source Code Library). record ascl:1202.015.

  • Google Scholar

• 40

  EddingtonA. S. (1926). The internal constitution of the stars. Cambridge University Press.

  • Google Scholar

• 41

  ErcolanoB.BarlowM. J.StoreyP. J.LiuX. W. (2003). MOCASSIN: a fully three-dimensional Monte Carlo photoionization code. MNRAS340, 1136–1152. 10.1046/j.1365-8711.2003.06371.x

  • CrossRef
  • Google Scholar

• 42

  FlockM.FromangS.GonzálezM.CommerçonB. (2013). Radiation magnetohydrodynamics in global simulations of protoplanetary discs. A&A560, A43. 10.1051/0004-6361/201322451

  • CrossRef
  • Google Scholar

• 43

  FrostholmT.HaugbølleT.GrassiT. (2018). Lampray: multi-group long characteristics ray tracing for adaptive mesh radiation hydrodynamics. arXiv e-prints , arXiv:1809.05541. 10.48550/arXiv.1809.05541

  • CrossRef
  • Google Scholar

• 44

  FryerC. L.RockefellerG.WarrenM. S. (2006). SNSPH: a parallel three-dimensional smoothed particle radiation hydrodynamics code. ApJ643, 292–305. 10.1086/501493

  • CrossRef
  • Google Scholar

• 45

  FryxellB.OlsonK.RickerP.TimmesF. X.ZingaleM.LambD. Q.et al (2000). FLASH: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes. ApJS131, 273–334. 10.1086/317361

  • CrossRef
  • Google Scholar

• 46

  GachesB. A. L.WalchS.WünschR.MackeyJ. (2023). Tree-based solvers for adaptive mesh refinement code FLASH - IV. An X-ray radiation scheme to couple discrete and diffuse X-ray emission sources to the thermochemistry of the interstellar medium. MNRAS522, 4674–4690. 10.1093/mnras/stad1206

  • CrossRef
  • Google Scholar

• 47

  Garcia-SeguraG.FrancoJ. (1996). From ultracompact to extended H II regions. ApJ469, 171. 10.1086/177769

  • CrossRef
  • Google Scholar

• 48

  GehmeyrM.MihalasD. (1994). Adaptive grid radiation hydrodynamics with TITAN. Phys. D. Nonlinear Phenom.77, 320–341. 10.1016/0167-2789(94)90143-0

  • CrossRef
  • Google Scholar

• 49

  GittingsM.WeaverR.CloverM.BetlachT.ByrneN.CokerR.et al (2008). The RAGE radiation-hydrodynamic code. Comput. Sci. Discov.1, 015005. 10.1088/1749-4699/1/1/015005

  • CrossRef
  • Google Scholar

• 50

  GloverS. C. O.Mac LowM.-M. (2007). Simulating the formation of molecular clouds. I. Slow formation by gravitational collapse from static initial conditions. ApJS169, 239–268. 10.1086/512238

  • CrossRef
  • Google Scholar

• 51

  GnedinN. Y.AbelT. (2001). Multi-dimensional cosmological radiative transfer with a Variable Eddington Tensor formalism. New Astron.6, 437–455. 10.1016/S1384-1076(01)00068-9

  • CrossRef
  • Google Scholar

• 52

  GonzálezM.AuditE.HuynhP. (2007). HERACLES: a three-dimensional radiation hydrodynamics code. A&A464, 429–435. 10.1051/0004-6361:20065486

  • CrossRef
  • Google Scholar

• 53

  GórskiK. M.HivonE.BandayA. J.WandeltB. D.HansenF. K.ReineckeM.et al (2005). HEALPix: a framework for high-resolution discretization and fast analysis of data distributed on the sphere. ApJ622, 759–771. 10.1086/427976

  • CrossRef
  • Google Scholar

• 54

  GrassiT.BovinoS.SchleicherD. R. G.PrietoJ.SeifriedD.SimonciniE.et al (2014). KROME - a package to embed chemistry in astrophysical simulations. MNRAS439, 2386–2419. 10.1093/mnras/stu114

  • CrossRef
  • Google Scholar

• 55

  GrondJ. J.WoodsR. M.WadsleyJ. W.CouchmanH. M. P. (2019). TREVR: A general N log²N radiative transfer algorithm. MNRAS485, 3681–3695. 10.1093/mnras/stz525

  • CrossRef
  • Google Scholar

• 56

  HarriesT. J.HaworthT. J.AcremanD.AliA.DouglasT. (2019). The TORUS radiation transfer code. Astronomy Comput.27, 63–95. 10.1016/j.ascom.2019.03.002

  • CrossRef
  • Google Scholar

• 57

  HasegawaK.UmemuraM. (2010). START: smoothed particle hydrodynamics with tree-based accelerated radiative transfer. MNRAS407, 2632–2644. 10.1111/j.1365-2966.2010.17100.x

  • CrossRef
  • Google Scholar

• 58

  HayesJ. C.NormanM. L. (2003). Beyond flux-limited diffusion: parallel algorithms for multidimensional radiation hydrodynamics. ApJS147, 197–220. 10.1086/374658

  • CrossRef
  • Google Scholar

• 59

  HofmannG. R. (1990). Who invented ray tracing?Vis. Comput.6, 120–124. 10.1007/bf01911003

  • CrossRef
  • Google Scholar

• 60

  HopkinsP. F.GrudićM. Y. (2019). Numerical problems in coupling photon momentum (radiation pressure) to gas. MNRAS483, 4187–4196. 10.1093/mnras/sty3089

  • CrossRef
  • Google Scholar

• 61

  IlievI. T.CiardiB.AlvarezM. A.MaselliA.FerraraA.GnedinN. Y.et al (2006). Cosmological radiative transfer codes comparison project - I. The static density field tests. MNRAS371, 1057–1086. 10.1111/j.1365-2966.2006.10775.x

  • CrossRef
  • Google Scholar

• 62

  IlievI. T.WhalenD.MellemaG.AhnK.BaekS.GnedinN. Y.et al (2009). Cosmological radiative transfer comparison project - II. The radiation-hydrodynamic tests. MNRAS400, 1283–1316. 10.1111/j.1365-2966.2009.15558.x

  • CrossRef
  • Google Scholar

• 63

  JauraO.GloverS. C. O.KlessenR. S.PaardekooperJ. P. (2018). SPRAI: coupling of radiative feedback and primordial chemistry in moving mesh hydrodynamics. MNRAS475, 2822–2834. 10.1093/mnras/stx3356

  • CrossRef
  • Google Scholar

• 64

  JiangY.-F. (2021). An implicit finite volume scheme to solve the time-dependent radiation transport equation based on discrete ordinates. ApJS253, 49. 10.3847/1538-4365/abe303

  • CrossRef
  • Google Scholar

• 65

  JiangY.-F. (2022). Multigroup radiation magnetohydrodynamics based on discrete ordinates including compton scattering. ApJS263, 4. 10.3847/1538-4365/ac9231

  • CrossRef
  • Google Scholar

• 66

  JiangY.-F.StoneJ. M.DavisS. W. (2012). A Godunov method for multidimensional radiation magnetohydrodynamics based on a variable Eddington tensor. ApJS199, 14. 10.1088/0067-0049/199/1/14

  • CrossRef
  • Google Scholar

• 67

  KannanR.VogelsbergerM.MarinacciF.McKinnonR.PakmorR.SpringelV. (2019). AREPO-RT: radiation hydrodynamics on a moving mesh. MNRAS485, 117–149. 10.1093/mnras/stz287

  • CrossRef
  • Google Scholar

• 68

  Kessel-DeynetO.BurkertA. (2000). Ionizing radiation in smoothed particle hydrodynamics. MNRAS315, 713–721. 10.1046/j.1365-8711.2000.03451.x

  • CrossRef
  • Google Scholar

• 69

  KlahrH. H.HenningT.KleyW. (1999). On the azimuthal structure of thermal convection in circumstellar disks. ApJ514, 325–343. 10.1086/306926

  • CrossRef
  • Google Scholar

• 70

  KlassenM.KuiperR.PudritzR. E.PetersT.BanerjeeR.BuntemeyerL. (2014). A general hybrid radiation transport scheme for star formation simulations on an adaptive grid. ApJ797, 4. 10.1088/0004-637X/797/1/4

  • CrossRef
  • Google Scholar

• 71

  KlepitkoA.WalchS.WünschR.SeifriedD.DinnbierF.HaidS. (2023). Tree-based solvers for adaptive mesh refinement code FLASH - III: a novel scheme for radiation pressure on dust and gas and radiative transfer from diffuse sources. MNRAS521, 160–184. 10.1093/mnras/stad385

  • CrossRef
  • Google Scholar

• 72

  KleyW. (1989). Radiation hydrodynamics of the boundary layer in accretion disks. I - numerical methods. A&A208, 98–110.

  • Google Scholar

• 73

  KolbS. M.StuteM.KleyW.MignoneA. (2013). Radiation hydrodynamics integrated in the PLUTO code. A&A559, A80. 10.1051/0004-6361/201321499

  • CrossRef
  • Google Scholar

• 74

  KrumholzM. R.KleinR. I.McKeeC. F.BolstadJ. (2007). Equations and algorithms for mixed-frame flux-limited diffusion radiation hydrodynamics. ApJ667, 626–643. 10.1086/520791

  • CrossRef
  • Google Scholar

• 75

  KuiperR.KlahrH.DullemondC.KleyW.HenningT. (2010). Fast and accurate frequency-dependent radiation transport for hydrodynamics simulations in massive star formation. A&A511, A81. 10.1051/0004-6361/200912355

  • CrossRef
  • Google Scholar

• 76

  LevermoreC. D. (1984). Relating Eddington factors to flux limiters. J. Quant. Spec. Radiat. Transf.31, 149–160. 10.1016/0022-4073(84)90112-2

  • CrossRef
  • Google Scholar

• 77

  LevermoreC. D.PomraningG. C. (1981). A flux-limited diffusion theory. ApJ248, 321–334. 10.1086/159157

  • CrossRef
  • Google Scholar

• 78

  LomaxO.WhitworthA. P. (2016). SPAMCART: a code for smoothed particle Monte Carlo radiative transfer. MNRAS461, 3542–3551. 10.1093/mnras/stw1568

  • CrossRef
  • Google Scholar

• 79

  LucyL. B. (1977). A numerical approach to the testing of the fission hypothesis. AJ82, 1013–1024. 10.1086/112164

  • CrossRef
  • Google Scholar

• 80

  LucyL. B. (1999). Computing radiative equilibria with Monte Carlo techniques. A&A344, 282–288.

  • Google Scholar

• 81

  MackeyJ.GreenS.MoutzouriM.HaworthT. J.KavanaghR. D.ZargaryanD.et al (2021). PION: simulating bow shocks and circumstellar nebulae. MNRAS504, 983–1008. 10.1093/mnras/stab781

  • CrossRef
  • Google Scholar

• 82

  MayerL.LufkinG.QuinnT.WadsleyJ. (2007). Fragmentation of gravitationally unstable gaseous protoplanetary disks with radiative transfer. ApJ Let.661, L77–L80. 10.1086/518433

  • CrossRef
  • Google Scholar

• 83

  MellemaG.IlievI. T.AlvarezM. A.ShapiroP. R. (2006). C ²-ray: a new method for photon-conserving transport of ionizing radiation. New Astron.11, 374–395. 10.1016/j.newast.2005.09.004

  • CrossRef
  • Google Scholar

• 84

  Melon FuksmanJ. D.MignoneA. (2019). A radiative transfer module for relativistic magnetohydrodynamics in the PLUTO code. ApJS242, 20. 10.3847/1538-4365/ab18ff

  • CrossRef
  • Google Scholar

• 85

  MihalasD.KleinR. I. (1982). On the solution of the time-dependent inertial-frame equation of radiative transfer in moving media to. J. Comput. Phys.46, 97–137. 10.1016/0021-9991(82)90007-9

  • CrossRef
  • Google Scholar

• 86

  MihalasD.MihalasB. W. (1984). Foundations of radiation hydrodynamics. Mineola, New York: Courier Corporation.

  • Google Scholar

• 87

  MinerboG. N. (1978). Maximum entropy Eddington factors. J. Quant. Spec. Radiat. Transf.20, 541–545. 10.1016/0022-4073(78)90024-9

  • CrossRef
  • Google Scholar

• 88

  MoensN.SundqvistJ. O.El MellahI.PoniatowskiL.TeunissenJ.KeppensR. (2022). Radiation-hydrodynamics with MPI-AMRVAC. Flux-limited diffusion. A&A657, A81. 10.1051/0004-6361/202141023

  • CrossRef
  • Google Scholar

• 89

  MurrayS. D.CastorJ. I.KleinR. I.McKeeC. F. (1994). Accretion disk coronae in high-luminosity systems. ApJ435, 631. 10.1086/174842

  • CrossRef
  • Google Scholar

• 90

  NayakshinS.ChaS.-H.HobbsA. (2009). Dynamic Monte Carlo radiation transfer in SPH: radiation pressure force implementation. MNRAS397, 1314–1325. 10.1111/j.1365-2966.2009.15091.x

  • CrossRef
  • Google Scholar

• 91

  OkamotoT.YoshikawaK.UmemuraM. (2012). ARGOT: accelerated radiative transfer on grids using oct-tree. MNRAS419, 2855–2866. 10.1111/j.1365-2966.2011.19927.x

  • CrossRef
  • Google Scholar

• 92

  PaardekooperJ. P.KruipC. J. H.IckeV. (2010). SimpleX2: radiative transfer on an unstructured, dynamic grid. A&A515, A79. 10.1051/0004-6361/200913821

  • CrossRef
  • Google Scholar

• 93

  PawlikA. H.SchayeJ. (2008). TRAPHIC - radiative transfer for smoothed particle hydrodynamics simulations. MNRAS389, 651–677. 10.1111/j.1365-2966.2008.13601.x

  • CrossRef
  • Google Scholar

• 94

  PetersT.BanerjeeR.KlessenR. S.Mac LowM.-M.Galván-MadridR.KetoE. R. (2010). H II regions: witnesses to massive star formation. ApJ711, 1017–1028. 10.1088/0004-637X/711/2/1017

  • CrossRef
  • Google Scholar

• 95

  PetkovaM.SpringelV. (2009). An implementation of radiative transfer in the cosmological simulation code GADGET. MNRAS396, 1383–1403. 10.1111/j.1365-2966.2009.14843.x

  • CrossRef
  • Google Scholar

• 96

  PriceD. J.BateM. R. (2009). Inefficient star formation: the combined effects of magnetic fields and radiative feedback. MNRAS398, 33–46. 10.1111/j.1365-2966.2009.14969.x

  • CrossRef
  • Google Scholar

• 97

  RamseyJ. P.DullemondC. P. (2015). Radiation hydrodynamics including irradiation and adaptive mesh refinement with AZEuS. I. Methods. A&A574, A81. 10.1051/0004-6361/201424954

  • CrossRef
  • Google Scholar

• 98

  RazoumovA. O.ScottD. (1999). Three-dimensional numerical cosmological radiative transfer in an inhomogeneous medium. MNRAS309, 287–298. 10.1046/j.1365-8711.1999.02775.x

  • CrossRef
  • Google Scholar

• 99

  RijkhorstE. J.PlewaT.DubeyA.MellemaG. (2006). Hybrid characteristics: 3D radiative transfer for parallel adaptive mesh refinement hydrodynamics. A&A452, 907–920. 10.1051/0004-6361:20053401

  • CrossRef
  • Google Scholar

• 100

  RitzerveldJ.IckeV. (2006). Transport on adaptive random lattices. Phys. Rev. E74, 026704. 10.1103/PhysRevE.74.026704

  • CrossRef
  • Google Scholar

• 101

  RobitailleT. P. (2011). HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code. A&A536, A79. 10.1051/0004-6361/201117150

  • CrossRef
  • Google Scholar

• 102

  RosdahlJ.BlaizotJ.AubertD.StranexT.TeyssierR. (2013). RAMSES-RT: radiation hydrodynamics in the cosmological context. MNRAS436, 2188–2231. 10.1093/mnras/stt1722

  • CrossRef
  • Google Scholar

• 103

  RosdahlJ.TeyssierR. (2015). A scheme for radiation pressure and photon diffusion with the M1 closure in RAMSES-RT. MNRAS449, 4380–4403. 10.1093/mnras/stv567

  • CrossRef
  • Google Scholar

• 104

  RosenA. L.KrumholzM. R.OishiJ. S.LeeA. T.KleinR. I. (2017). Hybrid Adaptive Ray-Moment Method (HARM²): a highly parallel method for radiation hydrodynamics on adaptive grids. J. Comput. Phys.330, 924–942. 10.1016/j.jcp.2016.10.048

  • CrossRef
  • Google Scholar

• 105

  RothN.KasenD. (2015). Monte Carlo radiation-hydrodynamics with implicit methods. ApJS217, 9. 10.1088/0067-0049/217/1/9

  • CrossRef
  • Google Scholar

• 106

  RozyczkaM. (1985). Two-dimensional models of stellar wind bubbles. I - numerical methods and their application to the investigation of outer shell instabilities. A&A143, 59–71.

  • Google Scholar

• 107

  SądowskiA.NarayanR.TchekhovskoyA.ZhuY. (2013). Semi-implicit scheme for treating radiation under M1 closure in general relativistic conservative fluid dynamics codes. MNRAS429, 3533–3550. 10.1093/mnras/sts632

  • CrossRef
  • Google Scholar

• 108

  SekoraM. D.StoneJ. M. (2010). A hybrid Godunov method for radiation hydrodynamics. J. Comput. Phys.229, 6819–6852. 10.1016/j.jcp.2010.05.024

  • CrossRef
  • Google Scholar

• 109

  ShuF. H. (1991). The physics of astrophysics. Volume 1: radiation. Mill Valley, CA: University Science Books.

  • Google Scholar

• 110

  SijoyC. D.ChaturvediS. (2015). TRHD: three-temperature radiation-hydrodynamics code with an implicit non-equilibrium radiation transport using a cell-centered monotonic finite volume scheme on unstructured-grids. Comput. Phys. Commun.190, 98–119. 10.1016/j.cpc.2015.01.019

  • CrossRef
  • Google Scholar

• 111

  SkinnerM. A.DolenceJ. C.BurrowsA.RadiceD.VartanyanD. (2019). FORNAX: a flexible code for multiphysics astrophysical simulations. ApJS241, 7. 10.3847/1538-4365/ab007f

  • CrossRef
  • Google Scholar

• 112

  SkinnerM. A.OstrikerE. C. (2013). A two-moment radiation hydrodynamics module in athena using a time-explicit Godunov method. ApJS206, 21. 10.1088/0067-0049/206/2/21

  • CrossRef
  • Google Scholar

• 113

  SmithA.KannanR.TsangB. T. H.VogelsbergerM.PakmorR. (2020). AREPO-MCRT: Monte Carlo radiation hydrodynamics on a moving mesh. ApJ905, 27. 10.3847/1538-4357/abc47e

  • CrossRef
  • Google Scholar

• 114

  StamatellosD.WhitworthA. P.BisbasT.GoodwinS. (2007). Radiative transfer and the energy equation in SPH simulations of star formation. A&A475, 37–49. 10.1051/0004-6361:20077373

  • CrossRef
  • Google Scholar

• 115

  SteinackerJ.BaesM.GordonK. D. (2013). Three-dimensional dust radiative transfer. Annu. Rev. Astronomy Astrophysics51, 63–104. 10.1146/annurev-astro-082812-141042

  • CrossRef
  • Google Scholar

• 116

  StoneJ. M.MihalasD.NormanM. L. (1992). ZEUS-2D: a radiation magnetohydrodynamics code for astrophysical flows in two space dimensions. III. The radiation hydrodynamic algorithms and tests. ApJS80, 819. 10.1086/191682

  • CrossRef
  • Google Scholar

• 117

  TeyssierR.CommerçonB. (2019). Numerical methods for simulating star formation. Front. Astronomy Space Sci.6, 51. 10.3389/fspas.2019.00051

  • CrossRef
  • Google Scholar

• 118

  TurnerN. J.StoneJ. M. (2001). A module for radiation hydrodynamic calculations with ZEUS-2D using flux-limited diffusion. ApJS135, 95–107. 10.1086/321779

  • CrossRef
  • Google Scholar

• 119

  ValdiviaV.HennebelleP. (2014). A fast tree-based method for estimating column densities in adaptive mesh refinement codes. Influence of UV radiation field on the structure of molecular clouds. A&A571, A46. 10.1051/0004-6361/201423720

  • CrossRef
  • Google Scholar

• 120

  VandenbrouckeB.WoodK. (2018). The Monte Carlo photoionization and moving-mesh radiation hydrodynamics code CMACIONIZE. Astronomy Comput.23, 40–59. 10.1016/j.ascom.2018.02.005

  • CrossRef
  • Google Scholar

• 121

  van der HolstB.TóthG.SokolovI. V.PowellK. G.HollowayJ. P.MyraE. S.et al (2011). CRASH: a block-adaptive-mesh code for radiative shock hydrodynamics—implementation and verification. ApJS194, 23. 10.1088/0067-0049/194/2/23

  • CrossRef
  • Google Scholar

• 122

  Vázquez-SemadeniE.González-SamaniegoA.ColínP. (2017). Hierarchical star cluster assembly in globally collapsing molecular clouds. MNRAS467, stw3229–1328. 10.1093/mnras/stw3229

  • CrossRef
  • Google Scholar

• 123

  ViauS.BastienP.ChaS.-H. (2006). An implicit method for radiative transfer with the diffusion approximation in smooth particle hydrodynamics. ApJ639, 559–570. 10.1086/499328

  • CrossRef
  • Google Scholar

• 124

  WalchS.GirichidisP.NaabT.GattoA.GloverS. C. O.WünschR.et al (2015). The SILCC (SImulating the LifeCycle of molecular Clouds) project - I. Chemical evolution of the supernova-driven ISM. MNRAS454, 246–276. 10.1093/mnras/stv1975

  • CrossRef
  • Google Scholar

• 125

  WhitehouseS. C.BateM. R. (2004). Smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation. MNRAS353, 1078–1094. 10.1111/j.1365-2966.2004.08131.x

  • CrossRef
  • Google Scholar

• 126

  WhitehouseS. C.BateM. R.MonaghanJ. J. (2005). A faster algorithm for smoothed particle hydrodynamics with radiative transfer in the flux-limited diffusion approximation. MNRAS364, 1367–1377. 10.1111/j.1365-2966.2005.09683.x

  • CrossRef
  • Google Scholar

• 127

  WhitworthA. P.ChapmanS. J.BhattalA. S.DisneyM. J.PongracicH.TurnerJ. A. (1995). Binary star formation: accretion-induced rotational fragmentation. MNRAS277, 727–746. 10.1093/mnras/277.2.727

  • CrossRef
  • Google Scholar

• 128

  WibkingB. D.KrumholzM. R. (2022). QUOKKA: a code for two-moment AMR radiation hydrodynamics on GPUs. MNRAS512, 1430–1449. 10.1093/mnras/stac439

  • CrossRef
  • Google Scholar

• 129

  WiseJ. H.AbelT. (2011). ENZO+MORAY: radiation hydrodynamics adaptive mesh refinement simulations with adaptive ray tracing. MNRAS414, 3458–3491. 10.1111/j.1365-2966.2011.18646.x

  • CrossRef
  • Google Scholar

• 130

  WolfireM. G.CassinelliJ. P. (1986). The temperature structure in accretion flows onto massive protostars. ApJ310, 207. 10.1086/164676

  • CrossRef
  • Google Scholar

• 131

  WünschR.WalchS.DinnbierF.SeifriedD.HaidS.KlepitkoA.et al (2021). Tree-based solvers for adaptive mesh refinement code FLASH - II: radiation transport module TreeRay. MNRAS505, 3730–3754. 10.1093/mnras/stab1482

  • CrossRef
  • Google Scholar

• 132

  WünschR.WalchS.DinnbierF.WhitworthA. (2018). Tree-based solvers for adaptive mesh refinement code FLASH - I: gravity and optical depths. MNRAS475, 3393–3418. 10.1093/mnras/sty015

  • CrossRef
  • Google Scholar

• 133

  ZhangW.HowellL.AlmgrenA.BurrowsA.BellJ. (2011). CASTRO: a new compressible astrophysical solver. II. Gray radiation hydrodynamics. ApJS196, 20. 10.1088/0067-0049/196/2/20

  • CrossRef
  • Google Scholar