Radiation transport methods in star formation simulations

Radiation transport plays a crucial role in star formation models, as certain questions within this ﬁeld 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 ﬂux-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.


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 reads where 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ν as where θ is the angle that n makes with a normal to the area element da.The first term on the r.h.s. of Equation 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 Equation 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 Equations 1 and 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 Equation 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 Equation 1 directly.Finally, Section §4 summarizes this review.

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 Equation 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 Equation 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, Equation 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.

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 P ν are the zeroth, first and second moments of the radiation field, respectively.The radiation pressure is often expressed as P ν = DE ν where D is the so called Eddington tensor.
Moments of Equation 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 are In 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 P ν are expressed in the coordinate system of the simulation (the laboratory frame).This problem is usually solved by either transforming Equations 4 and 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) where 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).

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 D = 1 3 I where I is the identity tensor.Assuming the radiation field is stationary, i.e. neglecting the first two terms in Equation 7, Wünsch and considering ∇ • P ν = 1 3 ∇E ν , 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 Equation 8 in the following way The flux limiter λ(R) is defined as a function of the normalized radiation energy density gradient in such a way that it ensures Equation 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 form which 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 form where 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 as and 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 codea 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  1. Radiation transport codes using the flux limited diffusion (FLD) approximation.The first column gives the reference to the work describing the code or algorithm; the second column give the name of the radiation transport code or RHD code (if they exist); the third column specifies the frame of reference in which the radiation energy equation is formulated (mixed, co-moving -CMF, or radiation energy advection terms are neglected); the fourth column tells whether the code is parallel using the domain decomposition; and the last column gives eventually an additional information.
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).

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 as where χ 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 Equations ( 6) and ( 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 Equations 6 and 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.

Variable Eddington Tensor methods
Another approach to solve moment Equations 4 and 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 Equation 1. Methods using it are referred to as Variable Eddington Tensor (VET) methods, and the Eddington tensor, D, is calculated by one of the methods described in §3.2 (ray   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 D 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 D and by separating the radiation and gas energy densities.The 1D RHD code TITAN (Gehmeyr and Mihalas, 1994) uses VET and calculates D by ray tracing.Gnedin and Abel (2001) developed an algorithm based on the Optically Thin Variable Eddington Tensor (OTVET) approximation where D is calculated in the optically thin regime.This is done by integrating contributions of radiation sources with attenuation factor 1/r 2 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 D and implemented it into 1D RHD code NIKE.This method was later implemented into code ATHENA by Jiang et al. (2012) where D 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 D 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.

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.

Ray tracing
The most straightforward method to solve the radiation transport Equation 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 Equation 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 the 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.1and 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 HARM 2 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 Equation 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).

Characteristics-based methods
A specific way of ray tracing are methods of characteristics, which solve the Radiative Transfer Equation (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 C 2 RAY (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 D (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 Wünsch 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.

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.

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 Carlobased 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).

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 4 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.Table 4. Radiation transport codes using ray tracing and methods of characteristics.The first column gives the reference to the work describing the code or algorithm; the second column give the name of the radiation transport code or RHD code (if they exist); the third column whether the code solves timedependent RTE; the fourth column tells whether the code is parallel using the domain decomposition; the fifth column gives how the code deals with the frequency dependence (grey, ion.-a single band of ionising photons, MG -multi-group approach, ν-dep -full frequency dependence); and the last column gives eventually an additional information.

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 selfinteraction, 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.

Wünsch
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 timedependent 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. (2006Iliev 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.

Table 2 .
Radiation transport codes using the M1 closure method.The meaning of columns is the same as in Table.1.

Table 3 .
Radiation transport codes using the Variable Eddington Tensor (VET) approach.The meaning of columns is the same as in Table. 1 tracing, long or short characteristics, or Monte Carlo).Subsequently, D is combined with the radiation moment equations to calculate the radiation energies and fluxes.