© The Authors 2025

1 Introduction

Astrophysical models often require radiative transfer to properly account for heating and cooling processes. This is generally challenging, as the radiative transfer equation (RTE) depends not only on space and time but also on the direction and energy of the transported photons. Moreover, it is often necessary to model the radiation and matter distributions simultaneously. In this context, many radiative transfer methods rely on simplified moment-based methods, which involve averaging the radiative intensity over all angles and often over frequency as well. Widely used examples include flux-limited diffusion (FLD; Levermore & Pomraning 1981), which involves a single equation for the total radiation energy density, and the two-moment M1 closure (Levermore 1984), which introduces an extra equation for the radiation flux.

While these methods often suffice in the diffusion regime, both the angle- and frequency-averaging strategies have drawbacks that can become significant depending on the problem. Frequency averaging may lead to inaccuracies when dealing with opacities that vary significantly with frequency (see, e.g., Dullemond et al. 2002). Furthermore, when angle-averaging the radiation intensity, physically motivated ad hoc assumptions must be made in order to close the resulting system of equations, which also introduces inaccuracies. In particular, FLD becomes excessively diffusive in the free-streaming regime, where it is unable to cast shadows (Hayes & Norman 2003). While the M1 closure does not suffer from this issue, it remains a hydrodynamical method for photon transport, and thus cannot propagate radiation in multiple directions at the same point in space. Therefore, both FLD and M1 lead to artificial interactions between crossing beams (e.g., Weih et al. 2020b), causing inaccuracies wherever such interactions become physically relevant.

These issues can be avoided by applying more refined methods that better capture the angular dependence of the radiation field. This includes approaches relying on an angular discretization of the radiation specific intensity, such as discrete ordinates methods (Jiang 2021, 2022; Ma et al. 2025), including long-(Magic et al. 2013; Frostholm et al. 2018; Stein et al. 2024) and short-characteristics (Davis et al. 2012) methods-sometimes integrated within variable Eddington tensor approaches (Jiang et al. 2012, 2014) – lattice Boltzmann methods (Asinari et al. 2010; Weih et al. 2020a), and reverse ray tracing (Wünsch et al. 2021), to name a few. Recent reviews of angle-discretized solvers can be found in Wünsch (2024) and Ma et al. (2025). Due to either the additional angular dimension or the iterative nature of some of these approaches, these are generally more computationally expensive than moment methods (Ma et al. 2025), and their applicability depends on the required spatial, angular, and frequency resolution. Some examples of applications to frequency-averaged radiative transfer in the context of planet-forming disks can be found, for example, in Bailey et al. (2024) (circumplanetary disks) and Zhang et al. (2024) (circumstellar disks).

Another alternative is provided by Monte Carlo (MC) methods, which involve the propagation of a large number of photon packets through an interacting medium. These schemes naturally account for the frequency and angular dependence of the radiation field and are highly accurate at computing temperature distributions in stationary cases. They are also ideal tools for reproducing scattering-dependent effects, as well as generating synthetic observations from simulations (e.g., Barraza-Alfaro et al. 2024). Coupling MC temperature computations with hydrodynamics (HD) or magnetohydrodynamics (MHD) via operator splitting is not trivial, and may introduce unphysical phenomena if the timescales of radiative and hydrodynamical processes are altered (Melon Fuksman & Klahr 2022). Instead, a viable option is to use MC methods to compute moments of the radiative intensity, which can then be naturally coupled with the HD/MHD equations (Foucart 2018; Smith et al. 2020).

Generally, the high cost of MC methods compared to the previously described approaches (Ma et al. 2025), especially in problems of high optical depth, makes them typically best suited for post-processing computations rather than dynamical simulations (see also Wünsch 2024). Yet, unlike dynamical simulations, such computations do not account for heating and cooling terms that depend on the velocity field, such as expansion and compression work. These processes cannot simply be included as heat sources in MC simulations, as this would leave out local cooling processes. As a result, even in stationary systems, MC post-processing temperature computations can be inaccurate near sources of significant compression and expansion, such as strong spiral waves in circumstellar disks (Krieger et al. 2025).

In standard benchmarks of starlight-irradiated circumstellar disks, several studies show minimal differences between the temperatures obtained with moment-based and MC methods (Pascucci et al. 2004; Flock et al. 2013; Melon Fuksman et al. 2021). However, in largely optically thick disk regions, Krieger et al. (2025) obtained larger midplane temperatures with FLD than with MC (see also Mignon-Risse et al. 2020). Besides the frequency averaging in FLD, a major reason behind this discrepancy is the fact that the flux entering the disk from the starlight-heated upper layers crosses the flux leaving the disk and cooling it down (Melon Fuksman & Klahr 2022). In moment methods, this results in an artificial interaction between the outward and inward fluxes, reducing the disk’s cooling efficiency. This inspires the goal of this paper: to prevent this interaction and develop a method suitable for dynamical simulations that produces accurate disk temperatures with minimal angular and frequency discretization.

To remove the interaction between photons entering and leaving the disk, it is immediate to think of the two-stream approximation, originally developed for stellar atmospheres (Schuster 1905; Eddington 1916) and widely used in 1D planetary atmosphere models (Malik et al. 2017; Mukherjee et al. 2023). This method replaces the radiative intensity with two fields corresponding to its integral in opposite halves of a plane; in other words, separate “upward” and “downward” fluxes are considered. In the context of circumstellar disks, we can expect that using two different fields for radiation entering and leaving the disk should prevent the artificial interaction produced in moment-based methods and yield more accurate temperature distributions.

A class of two-stream methods is given by half-moment (HM) closures, first introduced for 1D transport in Dubroca & Klar (2002) and extended to multigroup radiative transfer in Turpault et al. (2004). These methods are similar to moment-based closures, with the difference that the solid-angle integration of the radiation specific intensity is replaced with two integrals over hemispheres (half moments). The resulting system of equations is closed by imposing a maximum entropy constraint on the resulting fields, just as the M1 closure can be obtained using a closed angle integration. As in the M1 case, the obtained closure can be expressed in terms of simple relations between the involved half moments. The resulting transport equations are hyperbolic and can be solved using appropriate standard methods.

A generalization of HM methods to multidimensional transport is not straightforward, and few attempts have been made. The case of 2D transport was approached in Frank et al. (2006) by using partial moments integrated over four quadrants. The resulting closure requires a numerical inversion between Lagrange multipliers and partial moments, leading to obvious limitations in computational efficiency. Alternatively, Ripoll & Wray (2005) derived a 3D HM closure by assuming that (I) the partial fluxes are always parallel to the total flux and (II) the radiation specific intensity is isotropic in the plane perpendicular to the total flux. Though suitable for the diffusion regime, these approximations fail in cases involving obliquely crossing beams.

In this work, we introduce an HM closure suitable for multidimensional radiative transport. This closure can be derived by replacing the exact calculation of half moments via entropy maximization with an approximate functional form that reproduces the exact maximum-entropy closure with high accuracy, coinciding with it in the free-streaming and diffusion regimes. The half moments are defined by integrating over two halves of a plane whose normal, henceforth the “splitting direction”, is fixed in time. We also extend this closure to multigroup transport (see, e.g., Rosdahl et al. 2013; González et al. 2015); specifically, we divide the frequency domain into a series of discrete bins and solve separate HM equations for each of them, coupled only by radiation-matter interaction. As is shown later, this extension is essential to reproduce vertical temperature variations in irradiated disks.

We implemented these methods as extensions of the M1 radiative transfer code by Melon Fuksman et al. (2021), integrated into the open-source PLUTO code (Mignone et al. 2007) for HD and MHD. We employed a finite-volume scheme to compute the transport of the half moments, integrating the resulting system of equations via implicit-explicit (IMEX) methods (Pareschi & Russo 2005). This approach removes the need for globally implicit solvers, which is generally advantageous for parallelization. The closure is simple enough to compute the eigenvalues of the HM block exactly, which we used to derive and implement Riemann solvers.

For nonrelativistic problems, the light-crossing time of one grid cell is several orders of magnitude smaller than the fluid-crossing time, leading to a Courant-Friedrichs-Lewy (CFL) time step criterion that makes explicit integration prohibitively expensive. One way around this problem, which we employ in the present work, is to replace the speed of light, c, with an artificially reduced value, $\hat{c}$ (e.g., Gnedin & Abel 2001), which reproduces the correct physical limit provided that$\hat{c}$is nevertheless significantly larger than all relevant velocities in the problem. Alternatively, one could use a globally implicit approach incorporating both transport and source terms, avoiding the CFL limit.

This article is organized as follows. In Sect. 2, we introduce general maximum-entropy closures and define the closure implemented in this work. In Sect. 3 we summarize the systems of equations considered for HM and M1 radiative transfer. In Sect. 4 we describe the numerical implementation of our methods, and in Sect. 5 we test the code’s performance on a series of numerical benchmarks. In Sect. 6, we test the performance of the HM and M1 methods in modeling the temperature of passive irradiated protoplanetary disks in comparison with MC simulations. In Sect. 7, we discuss the accuracy and caveats of these disk models and characterize our methods’ computational efficiency. Finally, in Sect. 8 we summarize our conclusions. Additional details, benchmarks, and calculations, including a comparison between the exact and approximate HM closures, are included in the appendices.

2 Maximum-entropy half-moment closures

2.1 Full and half moments

Our goal is to derive a closed system of equations for half moments of the radiation specific intensity, I_v. This field depends on the photon frequency, v, the transport direction, $\hat{\mathit{n}}$, and the location in space and time, (t, r), where r denotes the spatial position. Our starting point is the radiative transfer equation, which determines the evolution of I_v as $$\frac{1}{c}{\mathrm{\partial }}_t{I}_v+\hat{\mathit{n}}\cdot \mathrm{\nabla }{I}_v={\eta }_v-\rho {\chi }_v{I}_v,$$(1)

where η_v and χ_v are, respectively, the medium’s emissivity and total opacity, while ρ is its mass density and c is the speed of light. The total opacity equals the sum of κ_v, the absorption opacity, and σ_v, the scattering opacity. Assuming coherent, isotropic scattering and thermal emission according to Kirchhoff’s law (Mihalas & Mihalas 1984), we assume an emissivity of the form $${\eta }_v(T)=\rho {\kappa }_v{B}_v(T)+\rho {\sigma }_v{J}_v,$$(2)

where J_v is the angle-averaged value of I_v, while B_v(T)= (2 h v³/c²)(e^hv/k_BT −1)⁻¹ is the Planck spectral radiance, with h and k_B the Planck and Boltzmann constants, respectively. Moment approaches rely on an angle integration of this equation multiplied by successive products of $\hat{\mathit{n}}$, resulting in a hierarchy of equations for radiative moments in which the flux of the n-th moment is the evolved field in the (n+1)-th equation. This procedure needs to be truncated at a finite number of equations in which the highest-order moment is determined as a function of the previous ones via a physically motivated closure relation.

Two-moment methods, such as the M1 closure by Levermore (1984), consist of closed systems of equations for the radiation energy density, its flux, and the radiation pressure tensor, which can be expressed as $$\begin{array}{rl}E& =\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v\oint \mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})\\ F^i& =\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v\oint \mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i\\ P^{ij}& =\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v\oint \mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i{n}^j,\end{array}$$(3)

respectively. Here we have rescaled the radiation flux by c in such a way that all of these fields are defined in energy density units.

Instead, here we aim to construct an HM scheme, for which we choose a splitting direction $\hat{\mathit{d}}$ and define + and – half moments by integrating Eq. (1) separately over the solid angle regions ${\mathcal{S}}_{+}$ and ${\mathcal{S}}_{-}$ verifying $\hat{\mathit{n}}\cdot \hat{\mathit{d}}>0$ and < 0, respectively. Specifically, we define $$\begin{array}{rl}& E_{\pm }=\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v{\int }_{{\mathcal{S}}_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})\\ & F_{\pm }^i=\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v{\int }_{{\mathcal{S}}_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i\\ & P_{\pm }^{ij}=\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v{\int }_{{\mathcal{S}}_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i{n}^j,\end{array}$$(4)

and so the full moments are recovered as $$\begin{array}{rl}E& =E_{+}+E_{-}\\ F^i& =F_{+}^i+F_{-}^i\\ P^{ij}& =P_{+}^{ij}+P_{-}^{ij}.\end{array}$$(5)

With these definitions, the + and − radiation moments are independent of each other. Thus, the separate integration of Eq. (1) over ${\mathcal{S}}_{+}$ and ${\mathcal{S}}_{-}$ yields subsystems for the + and − moments that remain independent in vacuum and are only coupled through the radiation-matter interaction terms. In this work, we do not consider moments beyond the pressure tensors, meaning that we need a closure for $P_{\pm }^{ij}$ as functions of (E_±, F_±).

2.2 Maximum-entropy closures

In the maximum-entropy¹ approach, reviewed in Müller & Ruggeri (1998), a closure relation is obtained by finding the distribution function that maximizes the entropy density of the radiation field, $$s\left(I_v\right)={\int }_0^{\mathrm{\infty }}\mathrm{d}v\oint \mathrm{d}\mathrm{\Omega }{h}_s\left(I_v\right)$$(6)

under the constraint that I_v reproduces a series of moments via defining relations such as Eq. (3) or (4). Here we use the expression for the entropy of a photon gas, with $$h_s\left(I_v\right)=\frac{2k_B{v}^2}{c^3}((\psi +1)\log (\psi +1)-\psi \log \psi ),$$(7)

where ψ = c² I_v/2hv³. This procedure, which can more generally be applied to classical or degenerate gases, yields the exact same closures as those obtained by enforcing the entropy inequality, ∂_ts+∇ · Φ_s ≥ 0 (Dreyer 1987; Müller & Ruggeri 1998), where Φ_s is the entropy flux. This provides a clear physical motivation for the entropy maximization approach.

2.3 The M1 closure

It is illustrative at this point to showcase how the M1 closure can be obtained by maximizing s under the constraints E = ⟨I_v⟩ and $\mathit{F}=⟨I_{\nu }\hat{\mathit{n}}⟩$, where ⟨·⟩ denotes $\frac{1}{c}{\int }_0^{\mathrm{\infty }}\mathrm{d}v\oint \mathrm{d}\mathrm{\Omega }$. To do this, we can apply the method of Lagrange multipliers in calculus of variations. Namely, we can introduce Lagrange multipliers a and b and replace s with the Lagrangian functional $$\mathcal{L}=s\left(I_v\right)+a(⟨I_v⟩-E)+\mathit{b}\cdot (⟨I_v\hat{\mathit{n}}⟩-\mathit{F}).$$(8)

We then zero the variation of $\mathcal{L}$ for small variations of I_v around the function that extremizes s, resulting in the equation ${\mathrm{\partial }}_{I_v}\mathcal{L}=0$, which is equivalent to $${\mathrm{\partial }}_{I_v}{h}_s+a+\mathit{b}\cdot \hat{\mathit{n}}=0.$$(9)

After rescaling a and b, it is straightforward to obtain that I_v is of the form $$I_v\propto \frac{v^3}{e^{av(1-b\cdot \hat{n})}-1},$$(10)

and therefore $${\int }_0^{\mathrm{\infty }}\mathrm{d}v{I}_v\propto \frac{1}{(1-\mathit{b}\cdot \hat{\mathit{n}}{)}^4}$$(11)

By equating E = ⟨I_ν⟩ and $\mathit{F}=⟨I_{\nu }\hat{\mathit{n}}⟩$, we can then invert the Lagrange multipliers as functions of E and F, after which the pressure tensor can be obtained as $\mathbb{P}=⟨I_V\hat{\mathit{n}}\hat{\mathit{n}}⟩$.

The resulting M1 closure can be expressed as $$P^{ij}=D^{ij}E,$$(12)

where the components of the Eddington tensor $\mathbb{D}$ are defined as follows: $$D^{ij}=\frac{1-\xi }{2}{\delta }^{ij}+\frac{3\xi -1}{2}{n}^i{n}^j.$$(13)

As derived by Levermore & Pomraning (1981), the Eddington factor ξ in this expression has the form $${\xi }_{\mathrm{M}1}(f)=\frac{5-2\sqrt{4-3f^2}}{3},$$(14)

where $f=\parallel \mathit{F}\parallel /E$ and $n^i=F^i/\parallel \mathit{F}\parallel$, while δ^ij is the Kronecker delta. The diffusion and free-streaming limits, given by f ≪ 1 and f → 1, correspond to D^ij = δ^ij/3 and D^ij = nⁱ n^j, respectively.

2.4 Half-moment closures for 1D transport

The same procedure can be used to find closures for the defined half moments. In this case, we need to maximize Eq. (6) for fixed$E_{\pm }={⟨I_{\nu }⟩}_{{\mathcal{S}}_{\pm }}$ and ${\mathit{F}}_{\pm }={⟨I_{\nu }\hat{\mathit{n}}⟩}_{{\mathcal{S}}_{\pm }}$ where the brackets’ subscripts indicate the angular region of integration. Since we have twice the number of constraints as before, we also need to duplicate the number of Lagrange multipliers, thereby defining the Lagrangian as $$\begin{array}{rl}\mathcal{L}=& s\left(I_v\right)+a_{+}({⟨I_v⟩}_{{\mathcal{S}}_{+}}-E_{+})+{\mathit{b}}_{+}\cdot ({⟨I_v\hat{\mathit{n}}⟩}_{{\mathcal{S}}_{+}}-{\mathit{F}}_{+})\\ & +a_{-}({⟨I_v⟩}_{{\mathcal{S}}_{-}}-E_{-})+{\mathit{b}}_{-}\cdot ({⟨I_v\hat{\mathit{n}}⟩}_{{\mathcal{S}}_{-}}-{\mathit{F}}_{-}).\end{array}$$(15)

Following the same steps as before, we obtain $${\mathrm{\partial }}_{I_v}{h}_s+a_{\pm }+{\mathit{b}}_{\pm }\cdot \hat{\mathit{n}}=0,\hat{\mathit{n}}\cdot \hat{\mathit{d}}\gtrless 0,$$(16)

which results in $${\int }_0^{\mathrm{\infty }}\mathrm{d}v{I}_v\propto \frac{1}{{(1-{\mathit{b}}_{\pm }\cdot \hat{\mathit{n}})}^4},\hat{\mathit{n}}\cdot \hat{\mathit{d}}\gtrless 0.$$(17)

Although this expression is almost identical to Eq. (11), its angle integration is substantially more complicated unless $\hat{\mathit{d}}$ and b_± can be assumed to be parallel, in which case integration is straightforward in spherical coordinates. Such is the case for slab and spherically symmetric geometries (1D transport). In these cases, the HM closure by Dubroca & Klar (2002) is obtained, where the Eddington factor ξ, relating the pressure tensor and energy as $P_{\pm }^{xx}=\xi E_{\pm }$, is computed as $${\xi }_{\mathrm{H}\mathrm{M}}(f)=\frac{8f^2}{1+6f+\sqrt{1+12f(1-f)}},$$(18)

where f = |F_±|/E_±.

As in the M1 case, this expression tends to 1/3 in the diffusion regime and 1 in the free-streaming regime (see Fig. 1). However, given that we have integrated over a hemisphere, now the diffusion regime corresponds to f ≈ 1/2; that is, $F_{\pm }^x\approx \pm \frac{E_{\pm }}{2}$, and not to f = 0. Moreover, f must satisfy $$f\ge 2/7,$$(19)

which results from the exact computation of the half moments from Eq. (17).

Fig. 1 Comparison of the Eddington factor in this work’s half-moment closure (ξHM) and the M1 closure (ξM1) as functions of $f=\parallel {\mathit{F}}_{\pm }\parallel /E$ and $\parallel \mathit{F}\parallel /E$, respectively. Vertical dotted lines indicate f = 2/7 (HM minimum flux) and 1/2 (isotropic limit in the HM case), while the horizontal dotted line indicates ξ = 1/3 (isotropic limit).

2.5 Half-moment closure used in this work

The HM approach can be immediately generalized to multidimensional transport. Once $\hat{\mathit{d}}$ is fixed, the fluxes defined by Eq. (4) can have arbitrary directions provided that ${\mathit{F}}_{\pm }\cdot \hat{\mathit{d}}\gtrless 0$. In this case, the maximum-entropy HM closure requires a numerical inversion of the Lagrange multiplier in Eq. (17). Instead, here we seek an expression that resembles the exact closure and matches it exactly in the free-streaming and diffusion regimes. With this goal in mind, we define a multidimensional HM closure as follows:

1. We assume the Eddington tensor depends on ξ as in Eq. (13) for f ≥ 1/2, with $\hat{\mathit{n}}={\mathit{F}}_{\pm }/\parallel {\mathit{F}}_{\pm }\parallel$.

2. We compute ξ as in Eq. (18), defining f = ∥F_±∥/E_±.

3. In a basis $\{\hat{\mathit{x}},\hat{\mathit{y}},\hat{\mathit{d}}\}$, with $\hat{\mathit{y}}$ perpendicular to F_±, it follows from Eqs. (4) and (17) that all components of $\mathbb{D}$ are nonnegative. We enforce this property by flooring the term (3ξ–1)/2 in Eq. (13) to zero for i ≠ j in such a basis. Expressing this in coordinate-independent form, we compute the Eddington tensor for f < 1/2 as $$\begin{array}{rl}{D}^{ij}=& \frac{1-\xi }{2}{\delta }^{ij}+\frac{3\xi -1}{2}[n^i{n}^j\\ & -(\hat{\mathit{n}}\cdot \hat{\mathit{d}})(n^i{d}^j+n^j{d}^i)+2(\hat{\mathit{n}}\cdot \hat{\mathit{d}}{)}^2{d}^i{d}^j],\end{array}$$(20)

where dⁱ are the components of $\hat{\mathit{d}}$. Details leading to this expression can be found in Appendix A.

This closure guarantees that both desired limits are correctly recovered as in the M1 closure, with the difference that now the isotropic limit D^ij = δ^ij/3 is obtained for f = 1/2, corresponding to ${\mathit{F}}_{\pm }=\pm \frac{E_{\pm }}{2}\hat{\mathit{d}}$. Moreover, the exact HM closure, coinciding with that in Dubroca & Klar (2002), is recovered when F^± and $\hat{\mathit{d}}$ are parallel. As is shown in Appendix A, the resulting Eddington tensor closely approximates its functional form in the exact closure with typical differences under ∼10%, vanishing near the mentioned limits.

Both in the exact and approximate HM closures, the radiation fields must verify a few physicality constraints. The first one is the usual limit preventing energy transport faster than the speed of light; namely, $$\parallel {\mathit{F}}_{\pm }\parallel \le E_{\pm },$$(21)

which stems from Eqs. (4). The second one is a generalization of Eq. (19) given by $${\left({\mathit{F}}_{\pm ,\perp }\right)}^2+{\left(\frac{7}{2}\right)}^2{\left({\mathit{F}}_{\pm ,\parallel }\right)}^2\ge E_{\pm }^2,$$(22)

where we have defined ${\mathit{F}}_{\pm ,\parallel }=({\mathit{F}}_{\pm }\cdot \hat{\mathit{d}})\hat{\mathit{d}}$ as the projection of F_± onto $\hat{\mathit{d}}$ and F_{±, ⊥} = F_±–F_{±, ∥}. We obtained this relation from a numerical computation of the half moments using Eq. (17) (Appendix A). Lastly, the fluxes must verify ${\mathit{F}}_{+}\cdot \hat{\mathit{d}}>0,{\mathit{F}}_{-}\cdot \hat{\mathit{d}}<0$, which stems from their definition. Although these constraints are naturally enforced in the diffusion regime, this is not necessarily the case for free-streaming transport. We investigate the consequences of enforcing them in Sections 5.2 and 6.

Computing the pressure tensor in these expressions is significantly less expensive than in the exact maximum-entropy HM closure since it can be done explicitly, as in the M1 closure. This makes the approximate HM closure ideal for implementation in numerical solvers for hyperbolic systems, such as those summarized in Section 4.

2.6 Multigroup extension

We now turn to a generalization to multigroup moments. We divide the frequency domain in N_g intervals or frequency groups, [$v_{\ell }^{\text{min}},v_{\ell }^{\text{max}}$], where ℓ = 1, …, N_g indicates the group index (this notation is kept throughout this work). For each group, we define half moments as $$\begin{array}{rl}{E}_{\ell ,\pm }& =\frac{1}{c}{\int }_{v_{\ell }^{\text{min}}}^{v_{\ell }^{\text{max}}}\mathrm{d}v{\int }_{S_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})\\ F_{\ell ,\pm }^i& =\frac{1}{c}{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\mathrm{d}v{\int }_{S_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i\\ P_{\ell ,\pm }^{ij}& =\frac{1}{c}{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\mathrm{d}v{\int }_{S_{\pm }}\mathrm{d}\mathrm{\Omega }{I}_v(t,\mathit{r},\hat{\mathit{n}})n^i{n}^j,\end{array}$$(23)

in such a way that the total moments are recovered by summing over groups and signs as $$E=\sum _{\begin{array}{c}\ell \\ s=+,-\end{array}}{E}_{\ell ,s},\mathit{F}=\sum _{\begin{array}{c}\ell \\ s=+,-\end{array}}{\mathit{F}}_{\ell ,s},\mathbb{P}=\sum _{\begin{array}{c}\ell \\ s=+,-\end{array}}{\mathbb{P}}_{\ell ,s}.$$(24)

To obtain a closure for the defined half moments, we would need to repeat the entropy maximization process and integrate the radiative intensity within each frequency group to obtain similar expressions to Eqs. (11) and (17), using different ± Lagrange multipliers per group. Unfortunately, as for the multigroup M1 closure, this frequency integration cannot be achieved analytically for finite integration intervals, which means that the resulting Eddington factor can only be obtained either via numerical inversion of the Lagrange multipliers in each group. This approach has been followed in Turpault (2002, 2003, 2005); Ripoll & Wray (2008) for the M1 closure and in Turpault et al. (2004) for the 1D HM closure.

Instead, here we follow a simplified approach often used in multigroup M1 methods (Vaytet et al. 2011; Rosdahl et al. 2013; Anninos & Fragile 2020), which consists of using the same closure for each group as in the frequency-integrated case. For each group, we compute ${\mathbb{P}}_{\ell ,\pm }$ in the same way as in Section 2.5, this time using f = ∥F_ℓ,±∥/E_ℓ,± and $\hat{\mathit{n}}={\mathit{F}}_{\ell ,\pm }/\parallel {\mathit{F}}_{\ell ,\pm }\parallel$ to obtain the Eddington factors via Eq. (18). This approach is justified by the property that the Eddington tensor must still have the form of Eq. (13) in the diffusion and free-streaming limits. We adopt this method as a first approximation and defer a precise comparison with the exact maximum-entropy closure to future studies.

3 Half-moment and M1 radiation MHD

3.1 Multigroup half-moment radiative transfer

Here we summarize the evolution equations used in this work, beginning with HM radiative transfer in the general multigroup case. These can be obtained by integrating Eq. (1) in frequency and solid angle in the hemispheres ${\mathcal{S}}_{\pm }$, resulting in the following system: $$\begin{array}{rl}& \frac{1}{\hat{c}}{\mathrm{\partial }}_t{E}_{\ell ,\pm }+\mathrm{\nabla }\cdot {\mathit{F}}_{\ell ,\pm }=-G_{\ell ,\pm }^0\\ & \frac{1}{\hat{c}}{\mathrm{\partial }}_t{\mathit{F}}_{\ell ,\pm }+\mathrm{\nabla }\cdot {\mathbb{P}}_{\ell ,\pm }=-{\mathit{G}}_{\ell ,\pm },\end{array}$$(25)

where the components of ${\mathbb{P}}_{\ell ,\pm }$ are obtained as functions of (E_{ℓ, ±}, F_{ℓ, ±}) as specified in Section 2.6. To increase the maximum time step allowed by the explicit integration of transport terms, here we have replaced c with the reduced speed of light $\hat{c}<c$.

When coupled to HD or MHD, the energy and momentum equations read $$\begin{array}{rl}{\mathrm{\partial }}_t\mathcal{E}+S_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^{\mathcal{E}}& =c\sum _{\ell }(G_{\ell ,+}^0+G_{\ell ,-}^0)-\mathrm{\nabla }\cdot {\mathit{F}}_{\mathrm{I}\mathrm{r}\mathrm{r}}\\ {\mathrm{\partial }}_t\mathit{m}+{\mathit{S}}_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^{\mathit{m}}& =\sum _{\ell }({\mathit{G}}_{\ell ,+}+{\mathit{G}}_{\ell ,-}),\end{array}$$(26)

where $\mathcal{E}$ and m are the gas total energy and momentum densities, respectively. Here we have only explicitly written terms corresponding to radiative processes, summarizing the rest (e.g., advection and pressure terms) as $S_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^{\mathcal{E}}$ and $S_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^m$. We have also added a source term for external irradiation sources producing a flux F_Irr, which will be used when modeling irradiated disks in Section (6).

Radiation-matter interaction terms, obtained by integrating the right-hand side of Eq. (1) (Appendix B), are computed for constant opacities as $$G_{\ell ,\pm }^0=\rho {\kappa }_{\ell }(E_{\ell ,\pm }-\frac{B_{\ell }(T)}{2})+\rho {\sigma }_{\ell }(E_{\ell ,\pm }-\frac{E_{\ell }}{2})$$(27a) $${\mathit{G}}_{\ell ,\pm }=\rho {\kappa }_{\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{B_{\ell }(T)}{4}\hat{\mathit{d}})+\rho {\sigma }_{\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{E_{\ell }}{4}\hat{\mathit{d}})$$(27b)

where ρ is the matter density and κ_ℓ, σ_ℓ, and χ_ℓ are the group absorption, scattering, and total opacities, respectively, while E_ℓ = E_ℓ,++E_ℓ,− and $$B_{\ell }(T)=\frac{4\pi }{c}{\int }_{v_{\ell }^{\text{min}}}^{v_{\ell }^{\text{max}}}{B}_v(T)\mathrm{d}v$$(28)

is the integrated Planck energy density, where T is the temperature. In the single-group case with [v^min, v^max) = [0, ∞), this function becomes$B(T)=a_R{T}^4$, where a_R is the radiation constant.

For frequency-dependent opacities, the coefficients in Eq. (27a) are replaced with their Planck-averaged values, $${\kappa }_{P\ell }=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\kappa }_v{B}_v(T)\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{B}_v(T)\mathrm{d}v},{\sigma }_{P\ell }=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\sigma }_v{B}_v(T)\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{B}_v(T)\mathrm{d}v},$$(29)

while Eq. (27b) is replaced with $${\mathit{G}}_{\ell ,\pm }=\rho {\kappa }_{R\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{B_{\ell }(T)}{4}\hat{\mathit{d}})$$(30)

for pure absorption opacities and $${\mathit{G}}_{\ell ,\pm }=\rho {\sigma }_{R\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{E_{\ell }}{4}\hat{\mathit{d}})$$(31)

for pure scattering. Here the subscript R denotes Rosseland averaging, defined as $${\kappa }_{R\ell }^{-1}=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\kappa }_v^{-1}\frac{\mathrm{\partial }{B}_v(T)}{\mathrm{\partial }T}\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\frac{\mathrm{\partial }{B}_v(T)}{\mathrm{\partial }T}\mathrm{d}v},{\sigma }_{R\ell }^{-1}=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\sigma }_v^{-1}\frac{\mathrm{\partial }{B}_v(T)}{\mathrm{\partial }T}\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\frac{\mathrm{\partial }{B}_v(T)}{\mathrm{\partial }T}\mathrm{d}v}.$$(32)

The case of frequency-dependent combined absorption and scattering is addressed in Appendix B.

3.2 Multigroup M1 radiative transfer

We now introduce the multigroup M1 equations solved in this work, which generalize the single-group system considered in Melon Fuksman et al. (2021). These are expressed as $$\begin{array}{r}\frac{1}{\hat{c}}{\mathrm{\partial }}_t{E}_{\ell }+\mathrm{\nabla }\cdot {\mathit{F}}_{\ell }=-G_{\ell }^0]]\\ \frac{1}{\hat{c}}{\mathrm{\partial }}_t{\mathit{F}}_{\ell }+\mathrm{\nabla }\cdot {\mathbb{P}}_{\ell }=-{\mathit{G}}_{\ell },\end{array}$$(33)

where the angle-averaged moments are defined for each group analogously to Eq. (23), this time using closed solid angle integrations. As in our HM closure and the M1 methods by Vaytet et al. (2011); Rosdahl et al. (2013); Anninos & Fragile (2020), we compute ${\mathbb{P}}_{\ell }(E_{\ell },{\mathit{F}}_{\ell })$ using the same M1 closure as in the frequency-integrated case (Section 2.3). The radiation fields are then coupled with HD and MHD as follows: $$\begin{array}{rl}{\mathrm{\partial }}_t\epsilon +S_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^{\epsilon }& =c\sum _{\ell }{G}_{\ell }^0-\mathrm{\nabla }\cdot {\mathit{F}}_{\mathrm{I}\mathrm{r}\mathrm{r}}\\ {\mathrm{\partial }}_t\mathit{m}+{\mathit{S}}_{\mathrm{H}\mathrm{D}/\mathrm{M}\mathrm{H}\mathrm{D}}^{\mathit{m}}& =\sum _{\ell }{\mathit{G}}_{\ell },\end{array}$$(34)

where the source terms are computed as $$\begin{array}{rl}{G}_{\ell }^0& =\rho {\kappa }_{P\ell }(E_{\ell }-B_{\ell }(T))\\ {\mathit{G}}_{\ell }& =\rho {\chi }_{R\ell }{\mathit{F}}_{\ell }.\end{array}$$(35)

4 Numerical implementation

To solve Equations (25) and (26), we implemented a numerical scheme that extends the gray M1 radiative transfer code by Melon Fuksman et al. (2021), which is freely distributed² as part of PLUTO 4.4-patch3. Our solution strategy is based on an explicit integration of advection terms via Godunov solvers and a locally implicit integration of the potentially stiff radiation-matter interaction terms. Operator splitting is used to integrate the transport equation of the radiation and MHD fields separately, each with their corresponding maximum time steps allowed by the CFL stability criterion. The integration of radiation transport and interaction terms is carried out simultaneously using the same IMEX schemes implemented for the M1 closure in Melon Fuksman & Mignone (2019). While here we only describe our implementation of the multigroup HM method, we use analogous methods to integrate the equations of both single-group and multigroup HM and M1 radiative transfer.

4.1 Explicit step

The explicit integration of radiation transport terms involves the solution of 2 N_g independent systems of hyperbolic PDEs of the form $${\mathrm{\partial }}_t{\mathcal{U}}_{\ell ,s}+\mathrm{\nabla }\cdot {\mathcal{F}}_{\ell ,s}=0$$(36)

corresponding to the left-hand side of Eq. (25), where ${\mathcal{U}}_{\ell ,s}=(E_{\ell ,s},{\mathit{F}}_{\ell ,s})$, with s = + or −. We solve these systems independently using the same reconstruct-solve-update strategy as in PLUTO, following a finite volume approach in which the discretized values of ${\mathcal{U}}_{\ell ,s}$ are interpreted as cell averages. In this implementation, $\hat{\mathit{d}}$ must be parallel to one of the coordinate axes in the chosen grid geometry; that is, it can only align with a Cartesian, cylindrical, or spherical unit vector.

To compute the fluxes, a set ${\mathcal{V}}_{\ell ,s}$ of primitive fields mapped from ${\mathcal{U}}_{\ell ,s}$ must be interpolated into the inter-cell interfaces. Here we chose ${\mathcal{V}}_{\ell ,s}={\mathcal{U}}_{\ell ,s}$. We also explored definitions such as ${\mathcal{V}}_{\ell ,\pm }=(E_{\ell ,\pm },{\mathit{F}}_{\ell ,\pm }\mp \left(E_{\ell ,\pm }/2\right)\hat{\mathit{d}})$ to properly compute the flux in diffusion problems where ∥F_ℓ∥ ≪ E_ℓ but ${\mathit{F}}_{\ell ,\pm }\approx \pm \left(E_{\ell ,\pm }/2\right)\hat{\mathit{d}}$. However, we found that the former definition of ${\mathcal{V}}_{\ell ,s}$ suffices to reproduce the correct diffusion flux in all cases we considered.

The reconstructed fields are then used to compute inter-cell fluxes via Riemann solvers. We implemented a Lax-Friedrichs (LF) solver and a Harten-Lax-van Leer (HLL) solver identical in form to those in Melon Fuksman & Mignone (2019). These solvers require knowledge of the maximum and minimum signal speeds along each direction, which we compute as the eigenvalues of the Jacobian matrices of Eq. (36). These can be explicitly obtained as functions of the radiation fields, as detailed in Appendix C.

The obtained fluxes are finally used to update the cell-averaged ${\mathcal{U}}_{\ell ,s}$, while signal speeds are used to update the time step according to the CFL condition. The latter step motivates our use of a reduced speed of light, which is not a necessary feature of the HM closure.

4.2 Implicit step

The implicit step handles the integration of radiation-matter interaction terms. By adding Equations (25) and (26) zeroing all transport and MHD terms, we can see that this step conserves the following modified total energy and momentum densities: $$\begin{array}{rl}{E}_{\mathrm{t}\mathrm{o}\mathrm{t}}& =\mathcal{E}+\frac{c}{\hat{c}}\sum _{\ell ,s}{E}_{\ell ,s}\\ {\mathit{m}}_{\mathrm{t}\mathrm{o}\mathrm{t}}& =\mathit{m}+\frac{1}{\hat{c}}\sum _{\ell ,s}{\mathit{F}}_{\ell ,s},\end{array}$$(37)

which coincide with the total energy and momentum only if $\hat{c}=c$. As extensively tested in various physical contexts (e.g., Skinner & Ostriker 2013; Melon Fuksman et al. 2021, 2024a; Melon Fuksman & Klahr 2022), these modified conservation laws do not alter the solutions for large enough $\hat{c}$. If external irradiation is included, −∇ · F_Irr is simply added to E_tot (see Appendix D), while if magnetic fields are included, the magnetic energy density B²/2 is subtracted from $\mathcal{E}$ as it does not intervene in the radiative energy exchange.

We use these conservation laws to simplify the equations solved during the implicit step, which become $$\begin{array}{rl}{\mathrm{\partial }}_t{E}_{\ell ,\pm }& =-\hat{c}{G}_{\ell ,\pm }^0\\ {\mathrm{\partial }}_t{\mathit{F}}_{\ell ,\pm }& =-\hat{c}{\mathit{G}}_{\ell ,\pm },\end{array}$$(38)

together with the conservation equtions of E_tot and m_tot. Unlike in the explicit step, all of these equations must be solved simultaneously due to the coupling introduced by the interaction terms. To speed up this process, we use the fact that momentum and flux variations resulting from radiative processes are much smaller than energy variations, allowing us to update all energies via either Newton or fixed-point methods while updating F_{ℓ, ±} in a fixed-point fashion. We describe these methods in Appendix D.

5 Numerical benchmarks

In this section, we test our implementation of HM radiative transfer in single-group problems designed to evaluate the main features of our methods. All tests were performed using the IMEX1 method in Melon Fuksman & Mignone (2019) to integrate the radiation fields, switching off HD advection. Unless otherwise stated, zero-gradient conditions were imposed on all fields, and fluxes were computed using the HLL solver and linear reconstruction with the second-order limiter by van Leer (1974). Additional tests involving absorption and scattering in the diffusion and free-streaming limits are presented in Appendix E.

5.1 Freely streaming beams

As a starting point, we tested the multidimensional transport scheme in the free-streaming regime. We studied the propagation of freely streaming beams oblique to the grid to evaluate their spread in the most diffusive case possible (González et al. 2007). We also performed crossing beam tests to demonstrate how the HM approach prevents the artificial interaction between converging beams in specific directions.

We defined the computational domain as [0, 2] ×[0, 2] in 2D Cartesian coordinates, with 300 cells in each direction. Fixing $\hat{\mathit{d}}=\hat{\mathit{x}}$ as the splitting direction, we initially set E₊ = E₋ = 10⁻⁸ in the isotropic transport regime; namely, ${\mathit{F}}_{\pm }=\pm \left(E_{\pm }/2\right)\hat{\mathit{d}}$. At the left and right boundaries, we injected two beams in the freestreaming regime with ∥F_±∥ = E_± = 1 as follows: $${\kappa }_{P\ell }=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\kappa }_v{B}_v(T)\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{B}_v(T)\mathrm{d}v},{\sigma }_{P\ell }=\frac{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{\sigma }_v{B}_v(T)\mathrm{d}v}{{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}{B}_v(T)\mathrm{d}v},$$(39)

The produced beams, shown in the bottom left panel of Fig. 2, are correctly advected at the speed of light maintaining their original direction. Trivially, the beams cross each other without interacting, as they correspond to separate, non-interacting fields. This is always the case when the beams have fluxes with opposite signs along the splitting direction. Conversely, if a third beam converging at the same point were introduced, it would inevitably lead to beam-crossing interaction, as its flux would necessarily have the same sign along the splitting direction as one of the other two beams.

In contrast, we also considered a case in which the crossing beams propagate with positive x velocity. We injected these beams at the left boundary with ∥F_±∥ = E_± = 1, as $${\mathit{F}}_{+}=\{\begin{array}{ll}(\sqrt{5/6},1/\sqrt{6})& \text{if\hspace{0.17em}}0.5\le y\le 0.6\\ (\sqrt{5/6},-1/\sqrt{6})& \text{if\hspace{0.17em}}1.4\le y\le 1.5\end{array}.$$(40)

In this case, as with the M1 closure, the beams artificially interact and merge at their crossing point (Fig. 2, right panel), as they are both described by the same field (E₊in this case).

For comparison, we repeated both tests with the M1 closure. The resulting energy distributions are shown in Fig. 2. The M1 closure always produces an interaction between crossing beams, while the HM closure does so only when both beams have positive x flux. Due to the chosen beam angles, the upward-transported M1 beams produce an overpressurized region at their colliding point that emits beams with both positive and negative y flux. Before any interaction is produced, both HM and M1 beams propagate across the domain, maintaining their initial directions with similar numerical spread. In both cases, numerical broadening decreases when the resolution is doubled. The merged beams in the right panels are similar for both closures because these coincide in the free-streaming limit. In both cases, sharp gradients, combined with grid noise, produce ripples in the energy distribution as radiation propagates away from the convergence point.

Fig. 2 Crossing beams produced with the M1 and HM closures. This test highlights how the HM closure removes the artificial interaction between crossing fluxes in selected directions. In the HM case, $\hat{\mathit{d}}=\hat{\mathit{x}}$ has been chosen as the splitting direction. Arrows indicate the directions of the injected beams.

5.2 Shadows

The ability to produce shadows is a key feature of the M1 closure, absent in the more diffusive FLD method (Hayes & Norman 2003). We tested this property in our HM code by reproducing a shadow test akin to those in Hayes & Norman (2003) and González et al. (2007), which also serves to consider a case of transition between the diffusion and free-streaming regimes.

We defined a 2D Cartesian domain as [−0.5, 0.5] cm ×[−0.12, 0.12] cm, with a resolution of 400 × 96 cells. We set the density such that it transitioned between a minimum value of ρ₀ = 10⁻⁵ g cm⁻³ and a maximum value ρ₁ = 10² g cm⁻³ with the functional form $$\rho (x,y)={\rho }_0+\frac{{\rho }_1-{\rho }_0}{1+e^{\mathrm{\Delta }}},$$(41)

where Δ = 10[(x/x₀)²+(y/y₀)²−1], with (x₀, y₀)=(0.1, 0.06) cm. This defines an elliptic overdense region at the center of the domain, acting as an obstacle for the incoming radiation. To this end, we defined a uniform absorption opacity κ = 10 cm² g⁻¹, resulting in a maximum horizontal optical depth of ∼200 across the obstacle and ∼10⁻⁴ through the background medium.

With the system initially in local thermal equilibrium (LTE) at T = 10 K, we injected a freely streaming radiation flux at the left boundary as $\mathit{F}=E_1\hat{\mathit{x}}$, with $E_1=2a_R{T}_1^4$, where T₁ = 10⁴ K. This parameterization of E₁ is unrelated to the actual temperature distribution of the injected photons, which is undefined.

We ran this simulation with the M1 and HM closures, switching off the update of the density and internal energy in both cases to focus solely on radiation transport. In the HM case, we chose $\hat{\mathit{d}}=\hat{\mathit{y}}$ as the splitting direction, considering equal + and − fluxes. As is shown in the top two panels of Fig. 3, a shadow is correctly formed in both cases behind the obstacle, and the energy distributions in both simulations are almost indistinguishable from each other. This is because, except for the shadowed region, radiation is either in the free streaming regime or close to LTE in the diffusion regime, and the Eddington tensor in both closures coincides in both regimes. In the HM case, this test also shows that radiation is correctly transported perpendicularly to the splitting direction $\hat{\mathit{d}}$.

We also used this test to evaluate the effect of enforcing the minimum-flux constraint given by Eq. (22). After each IMEX integration, we achieve this by modifying the flux in one of the following two ways wherever F_±does not satisfy the condition:

$${\mathit{F}}_{\pm }\to \alpha {\mathit{F}}_{\pm }\left({\mathrm{F}}_{min}\mathrm{v}1\right)$$(42a) $${\mathit{F}}_{\pm ,\parallel }\to \pm \left(E_{\pm }/2\right)\hat{\mathit{d}}\left({\mathrm{F}}_{min}\mathrm{v}2\right),$$(42b)

where $\alpha \in {\mathbb{R}}_{>0}$ is the minimum number that ensures the constraint is satisfied.

In Fig. 3, we show a comparison of the energy distributions obtained when this limiting is not enforced versus when each of these transformations is applied. When using Eq. (42a), the shadow looks almost the same as before, but the transition to the diffusion regime where radiation encounters the obstacle is smoother in the radial direction. This is because, in that region, interaction with matter preferentially reduces $F_{\pm }^x$ as radiation propagates through the obstacle. We can see this in both the M1 and HM cases in the (|f_y|, f_x) histograms shown in Fig. 4 when the constraint is not applied. This reduction is not consistent with Eq. (22), and so the flux is modified in that region. Once in the diffusion regime, with f ∼1/2, such a flux limiting is no longer required.

Conversely, only increasing $\left|F_{\pm }^y\right|$ as in Eq. (42b) leads to a similar transition in the x direction as before, but smooths out the transition in the y direction. Moreover, in this case, the shadow is lost because the method favors transport parallel to $\hat{\mathit{d}}$. We conclude that care must be taken in general when enforcing Eq. (22). While this is generally not necessary, it can be useful in beamcrossing scenarios that lead to deviations from this constraint (Section 6).

We lastly showcase the HM method’s ability to reproduce shadows cast by radiation emitted in different directions, either by separate sources or a single extended source. We considered the latter case by increasing the vertical domain to [− 0.2, 0.2] with a resolution of 160 × 160, setting the left boundary fluxes as

Fig. 4 Frequency of (|fy|, fx) values in each computational cell in the shadow tests shown in Fig. 3, normalized to a maximum of 1. The reduced fluxes are defined as (fx, fy)=F/E and F±/E±in the M1 and HM closures, respectively. Solid lines indicate the limits given by Eqs. (21) and (22), while dotted lines correspond to f = 1/2.

Fig. 5 Oblique shadows produced with the HM closure.

$${\mathit{F}}_{\pm }=\frac{E_1}{2}(\sqrt{19/20},\pm \sqrt{1/20}).$$(43)

In this way, radiation was simultaneously injected at the left boundary in two different directions. The resulting radiation energy distribution, shown in Fig. 5, is the sum of the E₊ and E₋distributions, each casting a shadow in its respective transport direction. As a result, the radiation field forms a shadow pattern with umbra, penumbra, and antumbra regions, characteristic of spatially extended light sources. These features cannot be reproduced with the M1 closure, which does not allow transport in multiple directions at the same location.

Fig. 3 Total energy density in the horizontal shadow test obtained with the M1 and HM closures. For the HM closure, distributions are shown for the cases where Eq. (22) is not enforced, enforced via Eq. (42a), and enforced via Eq. (42b), in the second, third, and fourth panels from the top, respectively.

6 Thermal modeling of protoplanetary disks

We now investigate the accuracy of multigroup M1 and HM radiative transfer in an astrophysical scenario involving the computation of the temperature in a protoplanetary disk irradiated by a central T Tauri star. In particular, we evaluate whether our HM method can improve the agreement with MC models. We modeled the disk as a stationary axisymmetric gas and dust distribution with a fixed dust-to-gas mass ratio. In the M1 and HM runs, we used tabulated frequency-dependent dust opacities to compute both the stellar irradiation absorption term via ray tracing and Planck and Rosseland opacities for reprocessed radiation, using N_g = 1, 3, or 22 frequency groups for the latter. At the same time, we adopted the same dust distribution, stellar parameters, and tabulated opacities to compute the temperature distribution with the MC radiative transfer code RADMC-3D (Dullemond et al. 2012). Implementation details and simulation parameters for all methods are provided in Appendix F.

6.1 M1 simulations

In Figs. 6 and 7, we show a comparison of the temperature distributions obtained with the described methods in 2D maps and at fixed radii, respectively. In the single-group M1 case, the midplane temperature is consistently overestimated in the optically thick region up to 44% with respect to the MC distribution. There are two main reasons for this: the frequency averaging, addressed later in this section, and the interaction between incoming and outgoing photons. The latter phenomenon, discussed in Section 3.3 of Melon Fuksman & Klahr (2022), can be summarized as follows. Radiation is absorbed near the Planck-averaged radial τ = 1 surface for stellar photons (henceforth the irradiation surface) and re-emitted at a lower temperature both away and toward the disk. The flux emitted toward the disk encounters the flux escaping in the vertical direction, merging with it as exemplified in Section 5.1. This reduces the flux leaving the disk, resulting in an overestimation of the midplane temperature. Additionally, this merger increases the radial flux produced by the radial temperature gradient in the optically thin region above the vertical τ = 1 surface. The resulting energy concentration in that region causes the deviation from vertical isothermality shown in Fig. 7 at large radii.

Above 10 au, we obtained a significant underestimation of the disk temperature by up to 37%. This discrepancy is caused by an underestimation of the absorption opacity resulting from its Planck averaging. We verified this by computing the ratio between the single-group Planck opacity and its energy-averaged value for N_g = 22, computed as κ_E(N_g)=∑_ℓ κ_{P ℓ} E_ℓ/∑_ℓ E_ℓ. In the regions where the temperature is underestimated, we obtained κ_P(N_g = 1)/κ_E(N_g = 22) < 1, with minimum values of ∼1/5. Consequently, for small N_g, the optically thin absorption of reprocessed radiation is underestimated, explaining why the temperature increases at large radii when increasing N_g.

In the single-group case, both the M1 and HM single-group runs produce vertically isothermal temperature distributions up to the irradiation surface. This is because, in stationary LTE, the flux source terms tend to zero, and therefore so does the divergence of the pressure tensor (Eqs. (25) and (33)). Given that ξ = 1/3 in this regime, the E and E_±gradients must also tend to zero. Since these energies are proportional to T⁴ in LTE, the resulting temperature distribution is vertically isothermal.

Conversely, the multigroup runs result in U-shaped vertical temperature distributions in the optically thick regions (see also Pavlyuchenkov & Akimkin 2025). An explanation for this is provided in Dullemond et al. (2002), and coincides with what we observe in our simulations. In these, the locations of the τ = 1 surfaces for vertically transported photons depend on opacity, resulting in departures from LTE at decreasing height for decreasing frequency. The disk cools via negative $G_{\ell }^0$ terms for the lowest-frequency groups, compensating for the escaping flux (Eq. (33)). To maintain thermal balance, the sum of all $G_{\ell }^0$ must be zero below the irradiation surface (Eq. (34)), meaning that the higher-frequency $G_{\ell }^0$ terms must increase. In other words, absorption must overcome emission at the higherfrequency groups in order to compensate for the energy loss at the low-frequency groups. This occurs through a steepening of the temperature gradient, which enhances the inward diffusive flux of the higher-frequency groups and consequently increases the inward-transported energy by exactly the amount needed to cancel the sum of all $G_{\ell }^0$.

The steepening of the temperature gradient in the multigroup runs results in lower midplane temperatures in optically thick regions compared to the single-group runs, leading to a better agreement with the MC distribution. With M1, we obtained maximum overestimations of the midplane temperature inside of 10 au of 44%, 23%, and 21% for N_g = 1, 3, and 22, respectively.

This agreement is aided by the fact that the frequency splitting prevents some of the interaction between photons entering and leaving the disk. Specifically, high-frequency inward-transported photons are emitted at the irradiation surface, while low-frequency outward-transported photons are emitted at the disk. However, the temperature difference between the midplane and upper layers is not large enough to prevent middle-range frequency groups from containing both photons entering and leaving the disk (Fig. F.2). Thus, regardless of the groups’ definitions, frequency splitting alone is not sufficient to prevent the interaction between incoming and outgoing fluxes for all groups. As a result, the midplane temperature is still significantly overestimated in the optically thick regions.

Fig. 6 Relative deviation of the disk temperatures from the MC distribution in the M1 and HM runs.

Fig. 7 Disk temperatures at fixed r values in the MC, HM, and M1 runs.

Fig. 8 Fluxes in the HM disk simulation with Ng = 3. Top: total radial flux. Bottom: θ component of the + flux for the highest-frequency group (ℓ = 3), also showing the direction of the + flux for the same group.

6.2 Half-moment simulations

In the HM case, the artificial interaction between outward and inward fluxes in the θ direction is absent, which improves the disk’s cooling efficiency. This significantly improves the match with the MC distribution, leading to maximum temperature overestimations of 6% and 8% inside of 10 au for N_g = 22 and 3, respectively. In the single-group case, the deviations are still large (up to 37%) since, as already discussed, the U shape cannot be recovered with a single frequency group. In conclusion, using 3 frequency groups suffices to obtain similar errors as with a much finer frequency discretization (Fig. 6).

At large radii, the agreement between HM and MC is better for N_g = 3 than for N_g = 22. This results from the underestimation of the Planck opacity in that region for small N_g, which leads to a temperature decrease that compensates for the energy accumulation at r = 100 au caused by boundary effects.

Despite the improved agreement with the MC distribution, the HM runs still overestimate the midplane temperatures in optically thick regions. Moreover, we also obtained a consistent temperature underestimation at the transition layers to the surface temperature.

To understand these differences, we need to look at the flux distribution in the HM runs. In Fig. 8, we show the distribution of the total radial flux for N_g = 3, together with the $F_{\ell ,+}^{\theta }$ distribution for the highest-frequency group (ℓ = 3). Also shown is the direction of the vector field F_3,+. This figure evidences that, as radiation escapes the disk, the vertical flux merges with the radial flux resulting from the radial temperature gradient. This results in a localized increase in the energy density above its LTE value below the irradiation surface, shown in Fig. 9 near θ−π/2 ∼0.07. As a result, according to Eq. (25), the vertical flux gradient must become negative. This produces a significant decrease in the outward flux, shown in Fig. 9, which reduces the cooling efficiency. This effect is further enhanced by the variation of the Eddington tensor as f reaches small values significantly below 1/2. In particular, the flux should in reality satisfy f ≥ 1/2 in that region, since f should only increase from its LTE value of 1/2 as energy is transported away from the disk in optically thin regions.

To minimize these effects, we enforced the physicality condition (22) via Eq. (42b). This slightly enhances the outward flux (see Fig. 9), reducing the midplane temperature and improving the agreement with the MC distribution. For comparison, if Eq. (42b) is not applied, we obtained maximum temperature overestimations within 10 au of 44%, 17%, and 13% for N_g = 1, 3, and 22, respectively. Lastly, we also found that ignoring the constraint D_{x z} ≥ 0 (Eq. (20)) slightly increased the disk temperature, worsening the match with MC by a few percent.

Fig. 9 Half-moment (+) flux and energy density at r = 1 au for the highest-energy group in the disk tests with Ng = 3, given in erg−3. Solid and dotted lines correspond to solutions obtained with and without enforcing Eq. (42b), respectively.

7 Discussion

7.1 Disk modeling

The multigroup HM temperatures are significantly closer to the MC distributions than those obtained with the M1 closure, even when using only three frequency groups. Compared to dynamical disk models relying on moment-based or single-group methods, this approach allows for more accurate modeling of HD and MHD disk processes self-consistently with radiative transfer without requiring high-resolution angular discretization. In particular, hydrodynamical instabilities are sensitive to temperature-dependent cooling timescales, and therefore the regions where they can operate have a strong dependence on the disk temperature (Pfeil & Klahr 2019; Melon Fuksman et al. 2024b). An accurate and self-consistent thermal modeling is also required for the precise modeling of disk observations (Ziampras et al. 2023, 2025a; Muley et al. 2024b), planetary migration (Guilera et al. 2021; Ziampras et al. 2024; Ziampras et al. 2025b), and shadowing effects (Melon Fuksman & Klahr 2022; Muley et al. 2024a).

While the application of Eq. (42b) in our HM models minimizes the artificial reduction of the outgoing flux, it makes the method more diffusive in the $\hat{\mathit{d}}$ direction, as characterized in Section 5.2. In the disk case, this prevents the formation of shadows cast by reprocessed radiation in the radial direction. It is important to note that this should not play a role in the modeling of phenomena produced by shadows cast by starlight (e.g., Melon Fuksman & Klahr 2022; Muley et al. 2024a), as these are produced by the irradiation source term. Still, it remains to be explored if the interaction between radial and vertical fluxes can be reduced while also minimizing diffusivity.

In some cases, radiation-HD simulations with more accurate temperature distributions than in this work may be needed. In particular, not only the temperature but also its derivatives control the local growth rates of hydrodynamical instabilities, such as the vertical shear instability and the convective overstability (Klahr 2024; Klahr et al. 2023). Naturally, one way to achieve this would be by means of frequency-dependent discrete ordinate methods. On the other hand, computing the temperature via MC methods provides an interesting alternative in scenarios where the temperature can be assumed constant or locally relaxed to an initial value. However, this approach is subject to the case-dependent drawbacks outlined in the introduction whenever compression and expansion work cannot be neglected. Alternatively, one could seek an approximate maximum-entropy closure for four quadrants, thus splitting the radial and vertical fluxes and removing their interaction. Additionally, frequency-dependent scattering, absent in this work to ease the method comparison, should also be included in realistic models. Further improvements could include obtaining separate different maximum-entropy closures for each frequency group (Section 2.6), and including non-LTE emissivities (e.g., Colombo et al. 2019).

7.2 Computational efficiency

Our method’s computational efficiency is subject to its implementation. In particular, with an IMEX method, the overhead depends almost exclusively on the value of the reduced speed of light, which determines the time step for the radiation solver. Conversely, reducing the speed of light would not be necessary in a globally implicit implementation, at the cost of worse parallel scalability than IMEX methods.

For the disk test shown here with $\hat{c}/c=2\times {10}^{-5}$, a typical minimum value for disk simulations with our chosen parameters, though problem-dependent, we found that solving the single-group radiation-HD problem advecting all three components of vector fields and using Newton’s method in the implicit step is ∼7 times more expensive with M1 and ∼16 times more expensive with HM than an equivalent HD simulation. The cost depends on the convergence speed of the implicit step, so we measure it after reaching radiative equilibrium, representing a typical case for a radiation-HD simulation. When employing fixed-point solvers, these factors are slightly reduced to 6 and 15 for M1 and HM, respectively. The difference between the Newton and fixed-point methods increases with N_g as the cost of these methods for large N_g is $O\left(N_g^3\right)$ and O(N_g), respectively. For N_g = 22 with HM, we obtained t_HM/t_HD ∼795 and 282 with Newton and fixed point methods, respectively, while for M1 we obtained t_M1/t_HD ∼120 and 65.

To characterize the scaling of both Newton and fixed-point methods for large N_g, we replicated the disk test with N_g up to 100 at a resolution of 100 × 100 using a single CPU (Fig. 10). We computed the cost of the integration step fixing the number of iterations to 1 in the implicit solver to avoid comparing runs with a different number of iterations. For low N_g, the cost of all methods is mostly due to the transport terms and therefore increases linearly with N_g, proportionally to the increase in the number of advected fields. The cost increase is slightly shallower for N_g < 3 due to the overhead of operations involving the HD fields, despite the fact that we do not include HD advection. With Newton’s method, the dependence on N_g becomes steeper for larger N_g, eventually reaching $t\propto N_g^3$ when the overhead is dominated by the computational cost of the Gaussian elimination algorithm employed by the implicit solver. For the HM method, this occurs for lower N_g than for M1 due to the higher cost of the implicit method, which iterates twice as many fields as in the M1 case. Since the average number of iterations in simulations is always ≥ 1, Fig. 10 represents a best-case scaling scenario where the transition to the $O\left(N_g^3\right)$ scaling occurs as late as possible. On the contrary, the fixed point method remains linear for large N_g, leading to substantial speedup compared to Newton’s method. In contrast, globally implicit methods (e.g., González et al. 2015) scale as $N_g^2$.

Our particular implementation via IMEX schemes can be compared with discrete ordinate methods such as that in Jiang et al. (2014), which also relies on IMEX integration with a reduced speed of light. This method requires advecting as many fields as directions defined for the angle discretization. In that work, 8 directions are sufficient to produce a slightly more diffusive crossing beam test compared to that shown in Fig. 2, obtained advecting 6 fields ($E_{\pm },F_{\pm }^x,F_{\pm }^y$). While we have employed a larger resolution, it is likely that more directions are required in discrete ordinate methods to produce sharply transported light beams comparable to that figure, which would increase the overall cost compared to the HM method. In axisymmetric protoplanetary disk models with a single-frequency group using a similar globally implicit method by Jiang (2021), Zhang et al. (2024) found convergence with 16 discrete directions, which is larger than the 6 radiation fields advected in our single-group disk simulations. Generally, the required angular resolution depends on the problem at hand; in particular, Zhu et al. (2020) used 80 discrete directions in 3D models of FU Ori’s accretion disks.

In general, discrete ordinate methods naturally outperform the HM closure’s accuracy. For example, in the beam crossing test, they would allow for the addition of a third beam, which in the HM case would result in artificial interaction. Yet, our HM method can be preferable when desiring to prevent such interactions in a single chosen direction, or when requiring a high-resolution frequency discretization that renders a high-resolution angle discretization prohibitively expensive.

Fig. 10 Computational cost of the radiation integration step using a single implicit iteration as a function of Ng, normalized to the fixed-point M1 cost for Ng = 1.

8 Conclusions

We introduced an HM closure for radiative transfer that extends previous 1D HM formulations to multidimensional transport. In this strategy, a splitting direction is chosen, and half moments are defined by integrating the radiation specific intensity over each hemisphere defined by it. A closure for the half moments can then be obtained through maximum-entropy constraints. We derived our HM closure by finding an approximation to this maximum-entropy closure that agrees with it in the diffusion and free-streaming limits. The resulting closure is easy to implement, and the transport terms can be integrated using Riemann solvers with exactly computed eigenvalues.

We implemented this method as a modification of the M1 radiation-HD/MHD module in the open-source PLUTO code, relying on IMEX schemes to handle both the implicit integration of radiation-matter interaction terms and the explicit integration of transport terms. We also extended both the HM and M1 codes to support multiple frequency groups. We implemented Newton and fixed-point methods for the implicit step. For a large number of frequency groups N_g, the cost of these methods scales as $N_g^3$ and N_g, respectively.

In a series of numerical benchmarks, we demonstrated that the HM closure behaves as expected in both free-streaming and diffusive transport, including shadow-casting scenarios where radiation transitions between these regimes. In particular, the HM closure removes the artificial interaction between beams – unavoidable in the M1 closure – in selected spatial directions. This is achieved with minimal angular discretization, as integrating the specific intensity over two hemispheres suffices to produce sharp beams in arbitrary directions. Conversely, this artificial interaction cannot be avoided when the beams propagate in directions corresponding to the same half moment; that is, when their fluxes are simultaneously positive or negative along the splitting direction.

In models of protoplanetary disks around T-Tauri stars, the HM method outperforms the M1 closure by removing the artificial interaction between fluxes entering and leaving the disk. The significantly smaller temperature discrepancies compared to MC simulations represent an improvement over typical state-of-the-art moment-based and single-group dynamical simulations of circumstellar disks. In multigroup simulations with 22 frequency groups, the HM closure limits the maximum temperature overestimation to 6% compared to MC, whereas this overestimation increases to 21% with the M1 closure. We found that as few as three groups are enough to reduce the temperature overestimation to 8% with the HM closure. The remaining discrepancies with the MC distributions result from interactions between vertical and radial fluxes, which could be avoided using discrete ordinate methods or in future extensions to partial-moment closures over four quadrants.

Acknowledgements

We thank Anton Krieger and Shangjia Zhang for valuable comments and discussions that helped us improve this manuscript, as well as W. Bear for his support. M.F. acknowledges support from the European Research Council (ERC), under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 757957). H.K. acknowledges funding from the European Research Council (ERC) via the ERC Advanced Grant “TiPPi – Turbulence, Pebbles and Planetesimals: The Origin of Minor Bodies in the Solar System” (project ID 855130) and from the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data,” funded by the German Excellence Strategy. Numerical simulations were run on MPIA’s VERA cluster, hosted at the Max-Planck Computing and Data Facility in Garching, Germany. Data analysis and plotting were performed using the Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), and PyPLUTO (Mattia et al. 2025) Python libraries.

Appendix A Comparison of the approximate and exact half-moment closures

Here we assess the accuracy of the HM closure introduced in this work compared with the exact maximum-entropy HM closure resulting from a numerical inversion of the Lagrange multipliers (Section 2.4). To do this, we need to compare the Eddington tensor $\mathbb{D}(E,\mathit{F})$ obtained in both approaches. Here, as in the rest of this appendix, we omit the subscripts ± to simplify the notation.

This comparison can be simplified by noting that both the exact and approximate closures satisfy the same transformation law for $\mathbb{D}$ when F is rotated about $\hat{\mathit{d}}$ according to a rotation matrix, $R_{\hat{d}}$; namely, $$\mathbb{D}(E,R_{\hat{d}}\mathit{F})=R_{\hat{d}}\mathbb{D}(E,\mathit{F})R_{\hat{d}}^{\mathrm{\top }}.$$(A.1)

Thus, a comparison for F in any plane containing $\hat{\mathit{d}}$ suffices to compare these closures for arbitrary (E, F). One way to prove this is to start from the computation of the half moments from Eq. (17) to show that the transformation $\mathit{b}\to R_{\hat{d}}\mathit{b}$ results in the transformations $$E\to E,\mathit{F}\to R_{\hat{d}}\mathit{F},\mathbb{D}\to R_{\hat{d}}\mathbb{D}{R}_{\hat{d}}^{\mathrm{\top }}.$$(A.2)

This means that the transformation of $\mathbb{D}$ as a function of the Lagrange multipliers when $\mathit{b}\to R_{\hat{d}}\mathit{b}$ is equivalent to its transformation as a function of F when $\mathit{F}\to R_{\hat{d}}\mathit{F}$. An alternative way is to view $R_{\hat{d}}$ as a transformation of coordinates, noting that the rotation symmetry about $\hat{\mathit{d}}$ of the radiation fields is not broken by the choice of the splitting direction.

To perform the comparison, we define a basis $\{\hat{\mathit{x}},\hat{\mathit{y}},\hat{\mathit{z}}\}$ of ${\mathbb{R}}^3$ such that, without loss of generality, $\hat{\mathit{d}}=\hat{\mathit{z}}$ and F is in the xz plane. Therefore, b must also be contained in that plane³, and so we can parameterize it as b = b(sin β, 0, cos β), with 0 ≤ b < 1 and 0 ≤ β ≤ π. This parameter reduction largely simplifies the calculation of the $\mathbb{D}$ components in the exact maximum-entropy closure, which can be obtained as D^ij = P^ij/E and displayed as functions of f_∥ = F_z/E and f_⊥ = F_x/E, computing all of these quantities as functions of (b, β).

While the $\mathbb{D}$ components are always nonnegative in this basis, $D_{xz}^{\text{approx}}$ defined by Eqs. (13) and (18) is negative for f < 1/2, with minimum values close to − 0.1. To improve the match between the approximate and exact closures, we enforce $D_{xz}^{\text{approx}}\ge 0$ in the defined basis by zeroing (3ξ–1)/2 in the diagonal terms for f < 1/2.

The resulting Eddington tensor for f < 1/2 can be expressed in coordinate-independent form as $${\mathbb{D}}^{\text{approx}}=\frac{1-\xi }{2}\mathbb{I}+\frac{3\xi -1}{2}[(\hat{\mathit{n}}\cdot \hat{\mathit{x}}{)}^2\hat{\mathit{x}}{\hat{\mathit{x}}}^{\mathrm{\top }}+(\hat{\mathit{n}}\cdot \hat{\mathit{d}}{)}^2\hat{\mathit{d}}{\hat{\mathit{d}}}^{\mathrm{\top }}],$$(A.3)

where $\mathbb{I}$ is the identity tensor and we have already zeroed D_xz. Since $\hat{\mathit{x}},\hat{\mathit{n}}$, and $\hat{\mathit{d}}$ are coplanar, we have $\hat{\mathit{x}}=k_1\hat{\mathit{n}}+k_2\hat{\mathit{d}}$. Together with the orthonormality conditions $\hat{\mathit{x}}\cdot \hat{\mathit{d}}=0$ and $\hat{\mathit{x}}\cdot \hat{\mathit{x}}=1$, this results in $\hat{\mathit{x}}=\pm (\hat{\mathit{n}}-(\hat{\mathit{n}}\cdot \hat{\mathit{d}})\hat{\mathit{d}})/{(1-(\hat{\mathit{n}}\cdot \hat{\mathit{d}}{)}^2)}^{1/2}$. Lastly, replacing this expression in Eq. (A.3), we obtain the general form of ${\mathbb{D}}^{\text{approx}}$ for f < 1/2 (Eq. (20)).

In Fig. A.1, we show the nonzero components of $\mathbb{D}$ and ${\mathbb{D}}^{\text{approx}}$, along with their differences, as functions of (f_∥, f_⊥). We obtain that both closures have similar functional forms, coinciding near the free-streaming and diffusion limits. Away from these limits, differences are typically under 10% except in some regions near the limit given by Eq. (22). In particular, D_zz is overestimated up to a factor of ∼2 in a narrow band near that limit, and D_xz is underestimated by a similar factor away from the isotropic limit when f ∼1/2.

Often the radiation field is either near the diffusion or the free-streaming regime (Fig. 4), where the approximate and exact maximum-entropy closures coincide. However, in some cases, it can reach (f_∥, f_⊥) regions where differences are maximal, particularly in scenarios involving beam crossing (Section 5.1). In such pathological cases, inaccuracies stem from the HM approach itself, regardless of whether the exact or approximate closures are used. Therefore, both closures should fail in similar scenarios and perform similarly whenever such effects are absent.

Appendix B Radiation-matter interaction terms

For each frequency, monochromatic HM radiative transfer equations can be obtained by multiplying Eq. (1) by ($1,\hat{\mathit{n}}$) and integrating over ${\mathcal{S}}_{\pm }$. This leads to the following system of equations: $$\begin{array}{rl}\frac{1}{\hat{c}}{\mathrm{\partial }}_t{E}_{v,\pm }+\mathrm{\nabla }\cdot {\mathit{F}}_{v,\pm }& =\rho {\kappa }_v(2\pi B_v(T)-E_{v,\pm })\\ & +\rho {\sigma }_v(\frac{E_v}{2}-E_{v,\pm })\\ \frac{1}{\hat{c}}{\mathrm{\partial }}_t{\mathit{F}}_{v,\pm }+\mathrm{\nabla }\cdot {\mathbb{P}}_{v,\pm }& =-\rho {\chi }_v{\mathit{F}}_{v,\pm }\\ & \pm \frac{\rho }{4}({\kappa }_{v}4\pi B_v(T)+{\sigma }_v{E}_v)\hat{\mathit{d}}\end{array}$$(B.1)

where E_{v, ±}, F_{v, ±}, and ${\mathbb{P}}_{v,\pm }$ are monochromatic half moments and E_v = E_v,++E_v,−. Here we have replaced $c\to \hat{c}$ in the factors multiplying the time derivatives. Opacity coefficients are computed in the fluid’s comoving frame, and relativistic corrections are ignored.

The frequency integration of the first of these equations is straightforward if we replace κ_v and σ_v with their Planck-averaged opacities, with the criterion that this is a valid approximation when I_v ∼ B_v(T) (Mihalas & Mihalas 1984). This results in the $G_{\ell ,\pm }^0$ term in Eq. (27a).

On the other hand, to average the second equation in such a way that the flux is accurate in the diffusion regime, we need to use Rosseland-averaged opacities. This is simple in the pure-absorption and pure-scattering cases, where it suffices to replace κ_v → κ_{R ℓ} or σ_v → σ_{R ℓ}, respectively. This results in the source terms in Eqs. (30) and (31). For combined absorption and scattering, we can follow similar steps to the usual derivation of the Rosseland opacity (Mihalas & Mihalas 1984) by zeroing the time derivatives and assuming I_v ∼ B_v(T) to compute the pressure tensor in the diffusion regime, resulting in $\mathrm{\nabla }\cdot {\mathbb{P}}_{v,\pm }\approx \frac{2\pi }{3}\mathrm{\nabla }{B}_v(T)$. In this way, we obtain the frequency-integrated source term $$\begin{array}{rl}{\mathit{G}}_{\ell ,\pm }& =\rho {\chi }_{R\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{B_{\ell }(T)}{4}\hat{\mathit{d}})\\ & \pm \frac{\rho {\chi }_{R\ell }}{4c}{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\mathrm{d}v\frac{{\sigma }_v}{{\chi }_v}(4\pi B_v(T)-E_v)\end{array}$$(B.2) $$\begin{array}{rl}& =\rho {\chi }_{R\ell }({\mathit{F}}_{\ell ,\pm }\mp \frac{E_{\ell }}{4}\hat{\mathit{d}})\\ & \pm \frac{\rho {\chi }_{R\ell }}{4c}{\int }_{v_{\ell }^{min}}^{v_{\ell }^{max}}\mathrm{d}v\frac{{\kappa }_v}{{\chi }_v}(E_v-4\pi B_v(T))\end{array}$$(B.3)

Fig. A.1 Nonzero components of the Eddington tensor in the exact maximum-entropy HM closure (Dij) and the approximate closure introduced in this work ($D_{ij}^{\text{approx}}$). The bottom row shows the difference between the two. Solid lines highlight the bounds given by Eqs. (21) and (22), while dotted lines correspond to f = 1/2. The isotropic and free-streaming limits correspond to (f∥, f⊥) = (1/2, 0) and f = 1, respectively. Here (f∥, f⊥) is defined as (fz, fx). The components shown are computed assuming zero y-flux, which suffices to represent all other flux directions, as is shown in Appendix A.

Since E_v ∼4 π B_v(T) in the diffusion regime and χ_Rℓ is small away from it, the integrals on the right-hand side can likely be ignored, leaving only the first term. Under this criterion, it should be convenient to use Eq. (B.2) for absorption-dominated problems and Eq. (B.3) for scattering-dominated ones. Alternatively, one could either Planck-average σ_ν/χ_ν in Eq. (B.2) or κ_ν/χ_ν in Eq. (B.3). While this is not necessary in this work as we only use frequency-dependent opacities for absorption-only problems (Section 6), future studies should compare these approaches in problems involving both frequency-dependent absorption and scattering.

Appendix C Signal speeds

Signal speeds are computed as the eigenvalues of the Jacobian matrices of Eq. (36). This calculation can be significantly simplified by writing the system in Cartesian coordinates and considering transport along single directions. Taking x as the chosen direction without loss of generality, we consider the following system: $$\frac{1}{\hat{c}}{\mathrm{\partial }}_t\mathcal{U}+{\mathrm{\partial }}_x{\mathcal{F}}^x=0,$$(C.1)

which has the Jacobian $J_{ij}^x=\mathrm{\partial }{\mathcal{F}}_i^x/\mathrm{\partial }{\mathcal{U}}^j$, where ${\mathcal{F}}^x=(F^x,P^{xx},P^{xy},P^{xz})$. Here, and in the rest of this section, we omit the group and sign subindices (ℓ, s).

Here we show the general result of this calculation for transport in an arbitrary direction ${\hat{\mathit{e}}}_i$ aligned with a coordinate axis. As in the M1 closure, the resulting eigenvalues for f ≥ 1/2 only depend on f and $\mu ={\hat{\mathit{e}}}_i\cdot \hat{\mathit{n}}$. We obtain a set of three eigenvalues, {λ₁, λ₂, λ₃}, where λ₂ is twice degenerate. In units of $\hat{c}$, these are of the form $$\begin{array}{rl}& {\lambda }_1={\beta }_1\mu -\sqrt{{\beta }_0+{\beta }_2{\mu }^2}\\ & {\lambda }_2=\frac{6f-1-\sqrt{\alpha }}{4f}\mu \\ & {\lambda }_3={\beta }_1\mu +\sqrt{{\beta }_0+{\beta }_2{\mu }^2},\end{array}$$(C.2)

where $$\begin{array}{rl}& \alpha (f)=1+12f(1-f)\\ & {\beta }_0(f)=\frac{1-\xi }{2}+(f+\frac{1-3\xi }{2f})\frac{{\xi }^{\mathrm{\prime }}}{2}\\ & {\beta }_1(f)=\frac{{\xi }^{\mathrm{\prime }}}{2}\\ & {\beta }_2(f)=\frac{{\left({\xi }^{\mathrm{\prime }}\right)}^2}{4}-f{\xi }^{\mathrm{\prime }}+\xi -{\beta }_0\end{array}$$(C.3)

Fig. C.1 Eigenvalues of the HM closure for F±aligned parallel (μ = 1), perpendicular (μ = 0), and antiparallel (μ = −1) to the propagation direction.

Here ξ is defined in Eq. (18), while ξ^′ denotes its derivative, $${\xi }^{\mathrm{\prime }}(f)=\frac{d\xi }{df}=\frac{1-\sqrt{\alpha }+2f(6-6f+\sqrt{\alpha })}{\alpha }.$$(C.4)

When D_xz = 0 is enforced for f < 1/2, the rotation symmetry of $\mathbb{D}$ in arbitrary directions is lost, leading to different eigenvalues for ${\hat{\mathit{e}}}_i$ parallel and orthogonal to $\hat{\mathit{d}}$. In both cases, the eigenvalues are not necessarily bound by $\pm \hat{c}$, and are complex for some values of f and $\hat{\mathit{n}}$, meaning that the system is not always strictly hyperbolic. Instead of computing the eigenvalues exactly, we maintain the presented form of the signal speeds also for f < 1/2, with the advantage that the resulting speeds are continuous at f = 0.5. We also enforce $\left|{\lambda }_i\right|\le \hat{c}$ to prevent superluminal transport for some μ values when f → 2/7. Improvements to these estimates would require an exhaustive analysis of the system’s characteristic waves for f < 1/2, which we defer to future investigations.

In highly optically thick regions, signal speeds must be reduced to prevent excessive numerical diffusion. To do this, we adopt the prescription in Eq. 47 of Sądowski et al. (2013).

Appendix D Implicit integration of source terms

The implicit methods presented here are analogous to those used in the PLUTO radiation module (see Melon Fuksman & Klahr 2022, Appendix B). The energy evolution equations (Eq. (38)), together with the conservation of E_tot (Eq. (37)), are used to write one of the energies, either E_{ℓ, s} or $\mathcal{E}$, as a function of the remaining ones, yielding a system of 2N_g equations with 2N_g unknowns. In our implementation, we can either choose V₁ = {E_{ℓ, s}} as the iterated variables in Newton’s method, or V₂ = {ϵ, E_{ℓ, s ≠ ℓ0, s0}} for a chosen (ℓ₀, s₀), where $ϵ=(\mathcal{E}-\rho {\mathit{v}}^2/2)/\rho$ is the specific internal energy. To ensure convergence, it is convenient to exclude the largest energy density from the iterated variables. This helps prevent small relative errors in a large variable from becoming large relative errors in smaller variables, which would require additional iterations to reach convergence.

We discretize Eq. (38) implicitly by evaluating source terms at the next time step as $$E_{\ell ,s}=E_{\ell ,s}^{\mathrm{\prime }}-\mathrm{\Delta }t\hat{c}{G}_{\ell ,s}^0,$$(D.1)

where we use tildes to denote fields evaluated at the previous time step, while untilded terms are evaluated at the next time step. Expanding this expression, we obtain $$0=E_{\ell ,s}^{\mathrm{\prime }}+{\zeta }_{\ell }^K\frac{B_{\ell }(T)}{2}-(1+{\zeta }_{\ell }^K+\frac{{\zeta }_{\ell }^{\sigma }}{2})E_{\ell ,s}+\frac{{\zeta }_{\ell }^{\sigma }}{2}{E}_{\ell ,-s},$$(D.2)

where we have defined ${\zeta }_{\ell }^k=\mathrm{\Delta }t\hat{c}\rho {\kappa }_{\ell }$ and ${\zeta }_{\ell }^{\sigma }=\mathrm{\Delta }t\hat{c}\rho {\sigma }_{\ell }$. If the set of variables V₁ is iterated, T(ϵ) is obtained as a function of V₁ using energy conservation ($E_{\text{tot}}=E_{\text{tot}}^{\mathrm{\prime }}$ and Eqs. (37)). Conversely, if the variables V₂ are iterated, we use energy conservation to obtain E_{ℓ0, s0} as a function of V₂.

The resulting system of equations is solved via a standard Newton method. Variable opacities, depending for instance on the temperature, are updated at each iteration and treated as constants for the computation of the Jacobian coefficients, in such a way that these have a simple polynomial form. To improve the scaling with N_g, we also implemented an alternative fixed-point algorithm by using Eq. (D.2) to express each E_{ℓ, s} as a function of T. The algorithm iterates all E_{ℓ, s} using these expressions and updates T via energy conservation.

After each iteration, fluxes are updated as $${\mathit{F}}_{\ell ,\pm }=\frac{1}{1+{\zeta }_{\ell }^{\chi }}[{\mathit{F}}_{\ell ,\pm }^{\mathrm{\prime }}\pm {\zeta }_{\ell }^K\frac{B_{\ell }(T)}{4}\hat{\mathit{d}}\pm {\zeta }_{\ell }^{\sigma }\frac{E_{\ell }}{4}\hat{\mathit{d}}],$$(D.3)

with ${\zeta }_{\ell }^{\chi }=\mathrm{\Delta }t\hat{c}\rho {\chi }_{\ell }$, while momentum is updated in such a way that the modified total momentum defined in Eq. (37) is conserved. If opacities are computed via frequency-averaging within each frequency group, the ζ functions in this expression use the Rosseland-averaged values of κ, σ, and χ. In that case, Eq. (D.3) is modified according to the form of the flux source terms in the pure-absorption, pure-scattering, and combined absorption and scattering cases (Appendix B). Conversely, Eq. (D.2) uses the Planck-averaged values of κ and σ.

If external irradiation is included, the corresponding source term is simply integrated by adding −Δt∇ · F_Irr to E_tot, as done in Muley et al. (2023). Compared to adding this term explicitly during the HD/MHD integration, this choice ensures faster convergence, as the initial guess for the iterated fields is closer to their converged values.