# Discontinuous Galerkin finite element method for the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes

Jérémy Perdigon, Marianne Faurobert, Gilles Niccolini

► To cite this version:

Jérémy Perdigon, Marianne Faurobert, Gilles Niccolini. Discontinuous Galerkin finite element method for the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes. Astronomy & Astrophysics - A&A, 2024, 690, pp.A192. 10.1051/0004-6361/202244322. hal-04732365

HAL Id: hal-04732365

https://hal.science/hal-04732365v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discontinuous Galerkin finite element method for the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes

J. Perdigon*, M. Faurobert, and G. Niccolini

Université de la Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France

Received 22 June 2022 / Accepted 16 August 2024

## ABSTRACT

**Context.** The study of the continuum radiative transfer problem inside circumstellar envelopes is both a theoretical and numerical challenge, especially in the frequency-dependent and multi-dimensional case. While approximate methods are easier to handle numerically, they often fail to accurately describe the radiation field inside complex geometries. For these cases, it is necessary to directly solve the radiative transfer equation numerically.

**Aims.** We investigate the accuracy of the discontinuous Galerkin finite element method (DGFEM hereafter) applied to the frequency-dependent two-dimensional radiative transfer problem, and coupled with the radiative equilibrium equation. We next used this method in the context of axis-symmetric circumstellar envelopes.

**Methods.** The DGFEM is a variant of finite element methods. It employs discontinuous elements and flux integrals along their boundaries, ensuring local flux conservation. However, as opposed to the classical finite element methods, the solution is discontinuous across element edges. We implemented this approach in a code and tested its accuracy by comparing our results with the benchmarks from the literature.

**Results.** For all the tested cases, the temperatures profiles agree within one percent. Additionally, the emerging spectral energy distributions (SEDs) and images, obtained by ray-tracing techniques from the DGFEM emissivity, agree on average within 5% and 10%, respectively.

**Conclusions.** We show that the DGFEM can accurately describe the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes. Consecutively the emerging SEDs and images are also well reproduced. The DGFEM provides an alternative method (other than Monte-Carlo methods for instance) for solving the radiative transfer equation, and it could be used in cases that are more difficult to handle with the other methods.

**Key words.** radiative transfer – methods: numerical – circumstellar matter

## 1. Introduction

The study of the continuum radiative transfer problem is crucial for the characterisation of circumstellar environments. Radiative processes play a major role in the determination of physical observables such as, for example, the temperature, abundances and velocity fields. The description of the radiation field is both a theoretical and numerical challenge, especially in the frequency-dependent and multi-dimensional case.

Several directions have been followed to tackle this problem. One approach involves approximate methods. It is usually done by assuming a particular form for the radiation field, most of the time based on symmetry arguments. The radiative transfer equation is then recast into a presumably simpler equation. This is the case, for example, for the Eddington approximation, in which the radiation is assumed to be a correction to an isotropic field, yielding an accurate description in the optically thick regime. For this approximation, the radiation is described by a simple linear diffusion equation. More sophisticated approximations were later developed, for example the flux-limited diffusion (Levermore & Pomraning 1981), which is asymptotically correct to both optically thin and thick regimes. In the latter case, the radiation is described by a non-linear diffusion equation.

While approximate methods are easier to handle numerically, they often fail to describe the radiation field accurately inside complex geometries (Kuiper & Klessen 2013). For these cases, it is necessary to numerically solve the radiative transfer equation directly (see Steinacker et al. 2013, for a thorough review of the different methods). A first approach is to approximate the transport operator with finite difference (e.g Steinacker et al. 2003), yielding a system of linear equations. This technique has the disadvantage of introducing spurious numerical oscillations and possible negative values for the intensity, due to the strong spatial and angular variations of radiation field. Other techniques, such as long and short characteristics (e.g Woitke et al. 2009), rely on the integral form of the radiative transfer equation. They are generally harder to implement than finite difference methods, are numerically demanding, and can also exhibit pathological behaviours such as negative values (Kunasz & Auer 1988). However, both approaches allow explicit error control. Finally, Monte-Carlo methods (e.g Wolf 2003) are amongst the most popular ones because they are easy and fast to implement, can handle complex geometries, and do not suffer from the same flaws as the other methods. However, they introduce noises that are inherent to their statistical nature. This class of methods additionally suffers from a reliable way of estimating the error, especially in optically thick regions where radiation is trapped inside the medium.

<footnote>* Corresponding author; jeremy.perdigon@oca.eu</footnote>
Aside from these techniques, finite element methods have already been used to solve the radiative transfer equation. A variant of it, the discontinuous Galerkin finite element method (Reed & Hill 1973, DGFEM hereafter), makes use of discontinuous elements and flux integrals along their boundaries, ensuring local flux conservation. However, as opposed to the classical finite element methods, the reconstructed solution is discontinuous across element edges. One of the main strengths of the method is its ability to produce a high-order numerical scheme, meaning a coarse computational grid can be used to achieve a small error. This feature is particularly interesting in the context of radiative transfer, where we are often limited by computational resources. Furthermore, the solution presents little to no oscillations, even in the cases where the solution displays a non-smooth behaviour (Nair et al. 2011), which often happens in the case of radiative transfer.

The DGFEM was successfully applied to the one-dimensional (1D) spherical transport problem, in the context of neutron transport (Machorro 2007), or more recently in grey stellar atmospheres (Kitzmann et al. 2016). Extensions to two-dimensional (2D) radiative transfer have been developed, in cylindrical coordinates, in the context of a non-local thermodynamic equilibrium (non-LTE) multilevel atom. Dykema et al. (1996) (and references therein) used a DGFEM with spatial linear elements for the evaluation of the Eddington tensor. The Eddington tensor was then subsequently used for closing the radiation field tensor moment system of equations. In another study, Cui & Li (2005) developed a DGFEM with a linear basis function in space and a step function in angle, for the computation of the radiative transfer in participating media.

In this study, we present the DGFEM applied to the frequency-dependent continuum radiative transfer problem, with isotropic scattering and coupled with the radiative equilibrium equation. We then used this method in the context of axis-symmetric dusty circumstellar environments. We show that this method can successfully determine the correct temperature profile, and allows for accurate images and spectral energy distributions (SEDs) to be computed by subsequent ray-tracing techniques. The paper is organised as follows: in Sect. 2 we describe the radiative transfer problem for axis-symmetric configurations. In Sect. 3 we present the DGFEM, and some of its numerical implementation features in Sect. 4. In Sect. 5, we compare this method against spherically and axis-symmetric benchmarks from the literature. Finally, in Sect. 6, we conclude and present some perspectives.

## 2. Description of the problem

We want to describe the radiation field inside an axis-symmetric circumstellar envelope. A central star of radius $R_*$ lies inside an inner cavity, free of matter, with a radius of $R_{in}$. The envelope spans from $R_{in}$ up to the outer radius $R_{out}$. We assume that no matter is present after this radius. We consider the material to be exclusively made of dust. We choose to describe the problem with the spherical coordinate system, both for the spatial and angular variables (see Fig. 1). In addition to the axial symmetry around the polar axis, we also assume a planar symmetry with respect to the equatorial plane, at $\Theta = \pi/2$. Given the symmetries, the domain of definition of radiation is $D \subset \mathbb{R}^4$, with $D \ni x = [r, \Theta, \mu = \cos \theta, \varphi] \in [R_{in}, R_{out}] \times (0, \pi/2] \times [-1, 1] \times [0, \pi]^1$. The radiation field is described by the time-independent

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration showing the coordinate system. We note that $(\hat{r}, \hat{\Theta}, \hat{\Phi})$ is the standard spherical basis. Given the symmetry around the polar axis, the radiation field, at a given position $r$, in a given direction $\Omega$ is a function of two spatial $(r, \Theta)$ and two angular $(\mu = \cos \theta, \varphi)$ coordinates.</figcaption></figure>

radiation transfer equation (for a thorough derivation see e.g. Pomraning 1973, II-5),

$$ \begin{gathered} \Omega \cdot \nabla I_v + \kappa_v^{\text{ext}} I_v = \eta_v, \\ \text{with } \Omega \cdot \nabla I_v = \mu \partial_r I_v + \frac{\sqrt{1-\mu^2} \cos \varphi}{r} \partial_{\Theta} I_v + \frac{1-\mu^2}{r} \partial_{\mu} I_v - \frac{\cot \Theta \sqrt{1-\mu^2} \sin \varphi}{r} \partial_{\varphi} I_v. \end{gathered} \quad (1) $$

The subscript $v$ denotes the frequency dependence, $I_v$ is the specific intensity, and $\Omega \cdot \nabla$ is the transport operator, which corresponds to a spatial derivative in the direction $\Omega$. We note that $\kappa_v^{\text{ext}} = \kappa_v^{\text{abs}} + \kappa_v^{\text{sca}}$ is the extinction coefficient, with $\kappa_v^{\text{abs}}$ and $\kappa_v^{\text{sca}}$ being the absorption and scattering coefficients, respectively. For dust species, we generally express these coefficients as $\kappa_v^{\text{abs,sca}} = C_v^{\text{abs,sca}} n$, where $C_v^{\text{abs,sca}}$ are the absorption and scattering optical cross-sections and $n$ is the number density. Dust grains are usually described as homogeneous spheres and their optical cross-sections are commonly computed with the help of Mie theory (Bohren & Huffman 1998). The emissivity $\eta_v$ includes the thermal emission (Kirchoff-Planck law), proportional to the Planck function, $B_v$, at the local temperature $T$, and a scattering term. For this study, we assume the scattering to be isotropic, for the purpose of our numerical tests (although it can be generalised to anisotropic scattering with no difficulty). The emissivity is then $\eta_v = \kappa_v^{\text{abs}} B_v(T) + \kappa_v^{\text{sca}} J_v$, with $J_v$, the mean specific intensity, which is the zeroth-order angular moment of the specific intensity $I_v$,

$$ J_v = \frac{1}{4\pi} \int_{4\pi} I_v d\Omega = \frac{1}{2\pi} \int_{0}^{\pi} d\varphi \int_{-1}^{1} I_v d\mu. \quad (2) $$

Additionally, the radiation field and the dust temperature are coupled via the equation of radiative equilibrium,

$$ 4\pi \int_{0}^{\infty} \kappa_{v}^{\text{abs}} B_{v}(T) d\nu = 4\pi \int_{0}^{\infty} \kappa_{v}^{\text{abs}} J_{v} d\nu. \quad (3) $$

<footnote n="1">¹ Note that $(\hat{r}, \hat{\Theta}, \hat{\Phi})$ is a plane of symmetry when axis-symmetry is assumed.</footnote>
The circumstellar matter is illuminated by a central star and it is customary to decompose the radiation field into two contributions (see e.g Steinacker et al. 2003), $I_v = I_v^* + I_v^{\text{env}}$, where $I_v^*$ is the direct unprocessed stellar radiation field attenuated by the circumstellar extinction, and $I_v^{\text{env}}$ is the radiation emitted by the envelope (either via thermal emission or scattering). Eq. (1) can then be recast into the following system of equations,

$$ \begin{cases} \Omega \cdot \nabla I_v^* + \kappa_v^{\text{ext}} I_v^* = 0, \\ \Omega \cdot \nabla I_v^{\text{env}} + \kappa_v^{\text{ext}} I_v^{\text{env}} = \eta_v = \kappa_v^{\text{abs}} B_v + \kappa_v^{\text{ sca}} (J_v^* + J_v^{\text{env}}), \end{cases} \quad (4) $$

with $J_v^* + J_v^{\text{env}} = J_v$. We note that in Steinacker et al. (2003), the thermal emission term was put on the right-hand side of the first Eq. (4). However, having it in the second equation is advantageous if $\kappa_v^{\text{abs}}$, $\kappa_v^{\text{ sca}}$, and $\kappa_v^{\text{ext}}$ are independent of temperature. In such cases, the equations decouple and $I_v^*$ can be solved and computed definitively. The first Eq. (4) can be integrated along a given ray, yielding the formal solution for the stellar contribution,

$$ I_v^* = \begin{cases} \Gamma_v^*(R_*, \mu) \exp\left(-\int_0^{s(R_{\text{in}})} \kappa_v^{\text{ext}}(s') ds'\right), & \text{if } \mu_* \le \mu \le 1, \\ 0, & \text{otherwise.} \end{cases} \quad (5) $$

The incident stellar radiation field is $\Gamma_v^*(R_*, \mu)$. Again, for the purpose of the numerical tests, we assume the star to radiate as a black body at the temperature $T_*$. Furthermore, $\mu_* = (1 - (R_*/r)^2)^{1/2}$ is the cosine of the angle subtended by the star at radius $r$, and $s(R_{\text{in}})$ is the distance between a given point $r$ and $R_{\text{in}}$, in the direction $\Omega$. The argument in the exponential is the opposite of the optical depth integrated along the ray. For dusty media, we usually have $r \gg R_*$, thus the star can be treated as a point source. In the point source approximation, $\mu_* \approx 1 - (R_*/r)^2/2$, and the stellar mean intensity can be expressed analytically as,

$$ J_v^* \approx \frac{1}{4} \left(\frac{R_*}{r}\right)^2 B_v(T_*) \exp\left(-\int_{R_{\text{in}}}^r \kappa_v^{\text{ext}}(r', \Theta) dr'\right). \quad (6) $$

In general, the integral in Eq. (6) can be carried out numerically, providing the stellar source term $J_v^*$ for Eq. (4).

To complete the description of the problem, we specify the boundary conditions for $I_v^{\text{env}}(r_s, \Omega)$ upon the surface of the domain $D$ located at $r_s$,

$$ I_v^{\text{env}}(r_s, \Omega) = \Gamma_v^{\text{env}}(r_s, \Omega), \quad \forall \hat{s} \cdot \Omega < 0, \quad (7) $$

where $\hat{s}$ is the unit vector normal to the surface of $D$, pointing outside of the domain. At the inner radius, $r = R_{\text{in}}$, $\hat{s} = -\hat{r}$, and the incident radiation, in a given direction $\Omega = \mu \hat{r} + \sqrt{1 - \mu^2} \cos \varphi \hat{\Theta} + \sqrt{1 - \mu^2} \sin \varphi \hat{\Phi}$, comes directly from the opposite point of the cavity,

$$ I_v^{\text{env}}(R_{\text{in}}, \Theta, \mu, \varphi) = I_v^{\text{env}}(R_{\text{in}}, \Theta', \mu', \varphi'), \quad \forall 0 < \mu < \mu_*, \quad (8) $$

with $(R_{\text{in}}, \Theta', \mu', \varphi')$ being the coordinates of the opposite point. Their derivation is given in Appendix A. On the outer edge, $r = R_{\text{out}}$, $\hat{s} = \hat{r}$, and we assume that there is no incident radiation upon the surface,

$$ I_v^{\text{env}}(R_{\text{out}}, \Theta, \mu, \varphi) = 0, \quad \forall \mu < 0. \quad (9) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Representation of the direction vector $\Omega$ at the equator. The dotted arrows represent the symmetric of $\Omega$ with respect to the equatorial plane $(\hat{r}, \hat{\Theta})$. The left and right 2D views allow the corresponding direction angles to be seen.</figcaption></figure>

At the equator, $\Theta = \pi/2$, $\hat{s} = \hat{\Theta}$, and the planar symmetry requires the radiation field to be, as shown in Fig. 2,

$$ J_v^{\text{env}}(r, \frac{\pi}{2}, \mu, \varphi) = I_v^{\text{env}}(r, \frac{\pi}{2}, \mu, \pi - \varphi), \quad \forall \varphi > \frac{\pi}{2}. \quad (10) $$

### 3. The radiative transfer equation with the discontinuous Galerkin finite element method

We present the DGFEM applied to the radiative transfer Eq. (1). We used some elements of notation from Kitzmann et al. (2016) and Hesthaven & Warburton (2007). For simplicity, we omitted the frequency subscript. The conservative form of Eq. (1) is,

$$ \nabla_x \mathbf{F} + \kappa^{\text{ext}} \tilde{\mathbf{I}} = \tilde{\eta}, \quad (11) $$

with $\tilde{\eta} = r^2 \sin \Theta \eta$. We introduced the variable $\tilde{\mathbf{I}} = r^2 \sin \Theta \mathbf{I}$, which is the quantity that is being conserved in Eq. (11). Using this quantity is important because it improves the stability of the numerical scheme, especially near the polar axis ($\Theta = 0$) where Eq. (1) is not defined when using the spherical coordinate system. The transport operator, $\nabla_x$, corresponds to the Cartesian divergence operator with respect to $\mathbf{x} = (r, \Theta, \mu, \varphi)$, applied to the flux vector $\mathbf{F}$ of $\tilde{\mathbf{I}}$,

$$ \mathbf{F} = a \tilde{\mathbf{I}} = \begin{pmatrix} a_r \\ a_\Theta \\ a_\mu \\ a_\varphi \end{pmatrix}, \quad \tilde{\mathbf{I}} = \begin{pmatrix} \frac{\mu}{\sqrt{1 - \mu^2} \cos \varphi / r} \\ \frac{(1 - \mu^2)/r}{\sqrt{1 - \mu^2} \sin \varphi / r} \\ -\cot \Theta \end{pmatrix}, \quad \tilde{\mathbf{I}}. \quad (12) $$

We decomposed the domain $D$ into $N = N_r \times N_\Theta \times N_\mu \times N_\varphi$ non-overlapping rectangular elements $D^{i,j,k,l}$, with $N_r, N_\Theta, N_\mu$ and $N_\varphi$ being the number of elements along the $r, \Theta, \mu$ and $\varphi$ coordinates, respectively. Each element is denoted with the help of four indexes $i, j, k$, and $l$, ranging from 0 to $N_r - 1$, $N_\Theta - 1$, $N_\mu - 1$ and, $N_\varphi - 1$. Inside each element $D^{i,j,k,l}$, we used the nodal representation and we approximated the local solution by a four-dimensional (4D) polynomial expansion,

$$ \tilde{I}_h^{i,j,k,l}(\mathbf{x}) = \begin{cases} \sum_{a=0}^{n_a-1} n_b-1 \sum_{b=0}^{n_b-1} n_c-1 \sum_{c=0}^{n_c-1} n_d-1 \tilde{I}_{a,b,c,d}^{i,j,k,l} h_{a,b,c,d}(\mathbf{x}), & \text{if } \mathbf{x} \in D^{i,j,k,l}, \\ 0, & \text{otherwise,} \end{cases} \quad (13) $$
with $n_a, n_b, n_c$, and $n_d$, being the number of nodes inside each element $D^{i,j,k,l}$, along the $r, \Theta, \mu$, and $\varphi$ coordinate, respectively. We note that $h_{a,b,c,d}$ is the 4D Lagrange polynomial, defined as,

$$h_{a,b,c,d}(\mathbf{x}) = h_a(r) h_b(\Theta) h_c(\mu) h_d(\varphi) \\ = \prod_{\substack{n_a=0 \\ \alpha \neq a}}^{n_a-1} \frac{r-r_\alpha}{r_a-r_\alpha} \prod_{\substack{\beta=0 \\ \beta \neq b}}^{n_b-1} \frac{\Theta-\Theta_\beta}{\Theta_b-\Theta_\beta} \prod_{\substack{\gamma=0 \\ \gamma \neq c}}^{n_c-1} \frac{\mu-\mu_\gamma}{\mu_c-\mu_\gamma} \prod_{\substack{\delta=0 \\ \delta \neq d}}^{n_d-1} \frac{\varphi-\varphi_\delta}{\varphi_d-\varphi_\delta}. \quad (14)$$

By definition, the coefficients $\tilde{I}_{a,b,c,d}^{i,j,k,l} = \tilde{I}_h^{i,j,k,l}(\mathbf{x}_{a,b,c,d})$ correspond to the value of $\tilde{I}_h^{i,j,k,l}$ at the nodes of coordinates $\mathbf{x}_{a,b,c,d} = (r_a, \Theta_b, \mu_c, \varphi_d)$. An example of an element $D^{i,j,k,l}$ with the associated nodes is shown in Fig. B.1. The global numerical approximation of the solution $I_h$ across the domain $D$ is formed by the sum of the $N$ piece-wise continuous solutions inside each element,

$$\tilde{I}(\mathbf{x}) \approx \tilde{I}_h(\mathbf{x}) = \sum_{i=0}^{N_r-1} \sum_{j=0}^{N_\Theta-1} \sum_{k=0}^{N_\mu-1} \sum_{l=0}^{N_\varphi-1} \tilde{I}_h^{i,j,k,l}(\mathbf{x}). \quad (15)$$

Now we shall introduce $\mathcal{R}_h$, the residual of Eq. (11),

$$\mathcal{R}_h(\mathbf{x}) = \nabla_x \cdot \mathbf{F}_h + \kappa^{\text{ext}} \tilde{I}_h - \tilde{\eta}, \quad (16)$$

with $\mathbf{F}_h = a \tilde{I}_h$. We form the classical Galerkin formulation (see e.g. Eq. (2.3) of Hesthaven & Warburton 2007) by requiring the residual to be orthogonal, inside each element, to the same set of functions used for the solution representation,

$$\int_{D^{i,j,k,l}} \mathcal{R}_h(\mathbf{x}) h_{a',b',c',d'}(\mathbf{x}) d^4\mathbf{x} = 0, \quad \forall D^{i,j,k,l}, \text{ and } \forall h_{a',b',c',d'}. \quad (17)$$

The divergence term that appears in Eq. (17) can be recast with the help of the divergence theorem (Green-Ostrogradsky), yielding the following so-called weak formulation of the radiative transfer equation Eq. (11):

$$\oint_{\partial D^{i,j,k,l}} \hat{s} \cdot \mathbf{F}^* h_{a',b',c',d'} d^3\mathbf{x} - \int_{D^{i,j,k,l}} \mathbf{F}_h \cdot \nabla h_{a',b',c',d'} d^4\mathbf{x} \\ + \int_{D^{i,j,k,l}} (\kappa^{\text{ext}} \tilde{I}_h - \tilde{\eta}) h_{a',b',c',d'} d^4\mathbf{x} = 0, \quad \forall D^{i,j,k,l}, \text{ and } \forall h_{a',b',c',d'}. \quad (18)$$

The first term, in the left-hand side of Eq. (18), is a surface integral, along the boundaries of $D^{i,j,k,l}$. Furthermore, $\hat{s}$ is the outward normal vector to the surface element, and $\mathbf{F}^*$ is an estimate of the flux at the cell interface, referred to as the numerical flux. It arises because the solution is not uniquely defined at the edge of the element, due to the discontinuous nature of the solution. We would like to emphasise that $D^{i,j,k,l}$ is a 4D rectangular element, and hence that each element has $2 \times 4$ boundary surfaces. The second and third terms in Eq. (18) are volume integrals and are purely local terms. Consequently, they only depend on the solution inside the element considered. For $\mathbf{F}^*$, several choices are possible, depending on the nature of the problem (Cockburn 2003). In our case, and as is commonly employed for transport problems, we used the Lax-Friedrichs numerical flux, defined as (e.g. on the radial right element edge where $\hat{s} = \hat{r}$),

$$\hat{r} \cdot \mathbf{F}^* = \frac{1}{2} (F_r(\tilde{I}^-) + F_r(\tilde{I}^+) - |a_r| (\tilde{I}^+ - \tilde{I}^-)), \quad (19)$$

where $\tilde{I}^-$ and $\tilde{I}^+$ denote the left and right values of $\tilde{I}$ at the element edge, respectively, and $a_r$ and $F_r$ are as defined in Eq. (12). The same form of expression holds for the other surfaces of the element.

Eq. (18) can be assembled into a system of $N' = N \times n$ equations with $n = n_a \times n_b \times n_c \times n_d$, relating the $\tilde{I}_{a,b,c,d}^{i,j,k,l}$ coefficients. The integrals are numerically estimated with the help of a quadrature formula. The choice of the quadrature, with the associated roots (nodes), is not unique and is usually problem-dependent (for an extensive review see e.g. Kopriva & Gassner 2010). In general, we can put the system of Eq. (18) in the form of the following linear system:

$$\mathcal{A} \tilde{I}_h = \mathcal{B}, \quad (20)$$

with $\tilde{I}_h$ being the vector of size $N'$, which contains the solution points for the full domain $D$, $\mathcal{A}$, a sparse matrix with a size of $N' \times N'$ coupling the elements of $\tilde{I}_h$, and $\mathcal{B}$ a vector of size $N'$ containing the emissivity term $\tilde{\eta}$.

## 4. Solution strategy and numerical considerations

The solution of the problem presented in Sect. 2 involves determining the radiation field, given by Eq. (4), coupled with the equation of radiative equilibrium, Eq. (3). For the test cases we consider in Sect. 5, the star was treated in the point source approximation, and hence we directly used Eq. (6) to compute $J_v^*$. The radiation field of the envelope $I_v^{\text{env}}$, described by the second Eq. (4), was solved with the DGFEM presented in Sect. 3. For simplicity, we omitted the envelope superscript but implicitly refer to this contribution below.

Because the radiative transfer problem is intrinsically highly dimensional, the size of the matrix $\mathcal{A}$ is often huge and solving Eq. (20) becomes numerically tedious. However, $\mathcal{A}$ is sparse, because each solution point, $\tilde{I}_{a,b,c,d}^{i,j,k,l}$, inside a given element, $D^{i,j,k,l}$, only depends on the other points inside the same element and the neighbouring ones, inside $D^{i \pm 1, j \pm 1, k-1, l+1}$. This property allows us to rewrite Eq. (20) as,

$$\begin{aligned} & \mathcal{A}^{i,j,k,l} \tilde{I}_h^{i,j,k,l} + \mathcal{A}^{i+1,j,k,l} \tilde{I}_h^{i+1,j,k,l} + \mathcal{A}^{i-1,j,k,l} \tilde{I}_h^{i-1,j,k,l} \\ & + \mathcal{A}^{i,j+1,k,l} \tilde{I}_h^{i,j+1,k,l} + \mathcal{A}^{i,j-1,k,l} \tilde{I}_h^{i,j-1,k,l} + \mathcal{A}^{i,j,k-1,l} \tilde{I}_h^{i,j,k-1,l} \\ & + \mathcal{A}^{i,j,k,l+1} \tilde{I}_h^{i,j,k,l+1} = \mathcal{B}^{i,j,k,l}, \quad \forall i, j, k, \text{ and } l. \end{aligned} \quad (21)$$

We note that $\mathcal{A}^{i,j,k,l}$ are the diagonal blocks of $\mathcal{A}$, with a size of $n \times n$ while $\mathcal{A}^{i \pm 1, j \pm 1, k-1, l+1}$ (with a size of $n \times n$) are the only non-zero, non-diagonal blocks. Furthermore, the elements $D^{i,j,k+1,l}$ and $D^{i,j,k,l-1}$ do not contribute because of the expression of the Lax-Friedrichs numerical flux, Eq. (19), with $a_\mu \ge 0$ and $a_\varphi \le 0$, $\forall x \in D$. We note that $\tilde{I}_h^{i,j,k,l}$ and $\mathcal{B}^{i,j,k,l}$ are the sub-vector of $\tilde{I}_h$ and $\mathcal{B}$, respectively, with a length of $n$, containing the local points in $D^{i,j,k,l}$.

To put $\mathcal{A}$ in a matrix form, we used a global index $\alpha = i N_\Theta N_\mu N_\varphi + j N_\mu N_\varphi + k N_\varphi + l$. An example of the structure of the matrix $\mathcal{A}$ with the associated blocks is displayed in Fig. 3. The formulation Eq. (21) avoids the storage and computation of the entire matrix $\mathcal{A}$, which reduces the computational effort.

In general, the simplest approach to obtain a solution from the problem is to solve Eq. (21) and to update the temperature with Eq. (3). Iterating between these two steps until convergence yields the solution to the problem. This procedure, commonly referred to as the $\Lambda$-iteration in the literature, becomes very slow and does not converge for large optical depths
<table><caption>**Table 1.** Summary of the main results.</caption><thead><tr><th></th><th>Ivezic et al. (1997)</th><th>Sect. 5.1</th><th>Pascucci et al. (2004)</th><th>Sect. 5.2</th></tr></thead><tbody><tr><td>Optical depth &tau;<sub>&nu;0</sub></td><td>1</td><td>10<sup>2</sup></td><td>10<sup>-1</sup></td><td>10<sup>2</sup></td></tr><tr><td>&epsilon;(T) (%)</td><td>0/0/0</td><td>0/0/1</td><td>0/0/0</td><td>1/0/2</td></tr><tr><td>&epsilon; (SED)<sub>i=12.5deg</sub> (%)</td><td>-/-/-</td><td>-/-/-</td><td>2/2/6</td><td>2/3/10</td></tr><tr><td>&epsilon; (SED)<sub>i=77.5deg</sub> (%)</td><td>2/2/8</td><td>1/2/2</td><td>2/2/5</td><td>3/4/24</td></tr></tbody></table>

**Notes.** Relative differences for the temperature profiles, $\epsilon(T)$, and for the SEDs, $\epsilon$ (SED), of the two test problems. The results are presented in the form mean(|$\epsilon$|) / std(|$\epsilon$|) / max(|$\epsilon$|) and rounded to the closest percent.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Example of the sparse structure of $\mathcal{A}$ with $N_r = N_\Theta = N_\mu = N_\phi = 2$. The squares represent the blocks of size $n \times n$ containing non-zero values. The matrix has bands corresponding to the coupling blocks $\mathcal{A}^{\alpha,\alpha}$, $\mathcal{A}^{\alpha,\alpha+1}$, $\mathcal{A}^{\alpha,\alpha+N_\Theta N_\mu N_\phi}$, $\mathcal{A}^{\alpha,\alpha+N_\mu N_\phi}$, and $\mathcal{A}^{\alpha,\alpha+N_\Theta N_\mu N_\phi}$.</figcaption></figure>

(Mihalas & Weibel-Mihalas 1999, VI-83). Our solution strategy is directly inspired from Perdigon et al. (2021). The key point of the method is to solve simultaneously instead of repetitively, Eq. (21) and Eq. (3). This strategy can be assimilated to an acceleration procedure to the usual $\Lambda$-iteration and yields satisfying results up to moderately thick envelopes. We however note that, as for the usual $\Lambda$-iteration, it converges very slowly for optically thick envelopes.

We proceeded with the following solution strategy: if we denote the iteration of the method by the superscript index $n$, we first (i) computed the stellar mean radiation field $J_v^*$ with Eq. (6) and set $[J_v^{\text{env}}]^n = 0$, as an initial condition. This allowed us to compute the initial temperature profile $T^0$, with the help of Eq. (3), where $J_v^0 = J_v^*$. Then (ii), for each frequency, we computed $[J_v^{\text{env}}]^{n+1}$ with the help of Eq. (21), which we rewrote, performing a block Gauss-Seidel sweep (see e.g Karniadakis & Kirby II 2003, VII-2),

$$ \begin{aligned} \mathcal{A}^{i,j,k,l} [\tilde{\mathbf{I}}_h^{i,j,k,l}]^{n+1} = & [\mathcal{B}^{i,j,k,l}]^n - \mathcal{A}^{i+1,j,k,l} [\tilde{\mathbf{I}}_h^{i+1,j,k,l}]^n \\ & - \mathcal{A}^{i-1,j,k,l} [\tilde{\mathbf{I}}_h^{i-1,j,k,l}]^{n+1} - \mathcal{A}^{i,j+1,k,l} [\tilde{\mathbf{I}}_h^{i,j+1,k,l}]^n \\ & - \mathcal{A}^{i,j-1,k,l} [\tilde{\mathbf{I}}_h^{i,j-1,k,l}]^{n+1} - \mathcal{A}^{i,j,k-1,l} [\tilde{\mathbf{I}}_h^{i,j,k-1,l}]^{n+1} \\ & - \mathcal{A}^{i,j,k,l+1} [\tilde{\mathbf{I}}_h^{i,j,k,l+1}]^n. \end{aligned} \quad (22) $$

We give in appendix B the expressions of $\mathcal{A}^{i,j,k,l}$ and of the right-hand side of Eq. (22). We note that Eq. (22) represents $N$ linear systems to solve. We solved the linear system Eq. (22) by direct inversion, using the Gauss elimination algorithm (see e.g Karniadakis & Kirby II 2003, IX-1). (iii) We computed $J_v^{\text{env}}$ and consequently updated the temperature via Eq. (3). The new temperature allowed us to update the right-hand side of Eq. (22) and to repeat steps (ii) and (iii), until convergence. We would like to point out that this scheme is not strictly identical to the usual $\Lambda$-iteration as we updated the temperature together with $I_v$ after each iteration index $n$. Finally, concerning the implementation of the boundary conditions, we directly followed the prescription from Kitzmann et al. (2016, Appendix A.6).

## 5. Numerical tests

The lack of analytic solutions for the radiative transfer problem, especially in the multi-dimensional and frequency-dependent case, limited our options for testing the validity of the method. In general, one has to compare numerical solutions with other ones from the literature. For this purpose, benchmark problems have been proposed. The first test case we consider, in Sect. 5.1, is the frequency-dependent radiative transfer problem inside a spherically symmetric envelope, from Ivezic et al. (1997). The second test case, in Sect. 5.2, is about the frequency-dependent radiative transfer inside an axis-symmetric envelope (disc), from Pascucci et al. (2004). For the latter test, we note that five different codes were used (including Monte Carlo and grid-based methods) to produce the benchmark and check the consistency of the results. We note that both of these tests are compatible with the boundary conditions presented in Sect. 2. A summary of the main results is presented in Table 1.

### 5.1. 1D spherically symmetric envelope

A point source, surrounded by a spherically symmetric envelope of dust, at radiative equilibrium, radiates as a black body at the temperature $T_* = 2500$ K. This envelope extends from the inner radius $R_{\text{in}}$ to the outer radius $R_{\text{out}} = 1000$ $R_{\text{in}}$. The inner radius was set so that its temperature would always be $T_{\text{in}} = T(R_{\text{in}}) = 800$ K. The number density profile, $n(r)$, is assumed to be a power law, of the form $n(r) = n_0 (R_{\text{in}}/r)^2$. The number density at the base of the disc, $n_0$, is determined by setting the value of $\tau_{v_0}$, which corresponds to the radial optical depth integrated through the envelope, at $v_0 = c/\lambda_0$, with $\lambda_0 = 1$ µm,

$$ \tau_{v_0} = \int_{R_{\text{in}}}^{R_{\text{out}}} \kappa_{v_0}^{\text{ext}} dr = C_{v_0}^{\text{ext}} \int_{R_{\text{in}}}^{R_{\text{out}}} n dr = C_{v_0}^{\text{ext}} n_0 R_{\text{in}} \left( 1 - \frac{R_{\text{in}}}{R_{\text{out}}} \right). \quad (23) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with $\tau_{v0} = 1$ (blue curve) and $\tau_{v0} = 10^2$ (orange curve). The cross marks represent the solution from this study and the dashed curves are from DUSTY (Ivezic et al. 1997). The lower panels show the relative differences between the two codes.</figcaption></figure>

We note that $C_{v0}^{\text{ext}} = C_{v0}^{\text{abs}} + C_{v0}^{\text{sca}}$ is the extinction cross-section coefficient at $v_0$. In this test we consider $\tau_{v0} = 1$ and $10^2$ to correspond to a moderately thin and thick envelopes, respectively. The absorption and scattering cross-sections, $C_v^{\text{abs}}$ and $C_v^{\text{sca}}$, feature a bi-linear behaviour in log-log scaling, with a constant profile for $v \ge v_0$ and a power-law dependence $\propto (v/v_0)^\alpha$, for $v \le v_0$, with $\alpha = 1$ and 4 for absorption and scattering, respectively. This dependence aims to mimic the behaviour of spherical astronomical silicate grains. We set the value of the thermal coupling parameter $C_v^{\text{abs}}/C_v^{\text{ext}}$ to be 1/2. The benchmarks from Ivezic et al. (1997) were reproduced with version 2 of DUSTY<fnref n="2" />. Although the envelope is spherically symmetric, we used a grid of $N = N_r \times N_\Theta \times N_\mu \times N_\varphi = 16^4$ elements in our DGFEM code, with each elements containing $n_a \times n_b \times n_c \times n_d = 3 \times 2 \times 3 \times 3$ nodes, as pictured in Fig. B.1. For the radial coordinate, the cell edges were logarithmically spaced, to account for the important dynamic of the solution with respect to the radius. The frequency grid consists of 60 logarithmically spaced points, ranging from $\lambda = 10^{-2}$ µm to $\lambda = 3.6 \times 10^4$ µm.

The temperature profiles, $T$, and the SEDs, $\lambda F_\lambda/F$ with $F = \int_0^\infty F_\lambda d\lambda$, of the envelope are shown in Fig. 4. For the DGFEM code, the temperature that is shown corresponds to the mean radial profile across all angular points $\Theta$. Additionally, the normalised SEDs were computed with the help of a ray-tracing module we present in Appendix C. We subsequently explain the reasons for such a choice in Sect. 5.3. The spatial and frequency grids differ between both codes and we performed a linear interpolation (in log-log scaling) of the DUSTY profiles at our grid points, in order to do the comparison. We observed a good agreement between the two codes. On average, the absolute relative differences stay below 1% for the temperatures and 2% for the SEDs.

## 5.2. 2D axis-symmetric envelope

A point source surrounded by an axis-symmetric disc of dust, at radiative equilibrium, radiates as a black body at the temperature

$T_\star = 5800$ K. This disc extends from the inner radius $R_{\text{in}} = 1$ AU to the outer radius $R_{\text{out}} = 1000$ AU. The density profile $n(r, \Theta)$ is assumed to be

$$n(r, \Theta) = n_0 \frac{r_d}{r \sin\Theta} \exp\left(-\frac{\pi}{4} \left(\frac{r \cos\Theta}{h}\right)^2\right), \quad (24)$$

with $h = z_d (r \sin\Theta/r_d)^{9/8}$, $r_d = R_{\text{out}}/2$, and $z_d = R_{\text{out}}/8$. This density law is characteristic of a Keplerian disc hydrostatically supported in the vertical direction, assuming a constant temperature along this direction (Chiang & Goldreich 1997). Again, $n_0$ was determined to set $\tau_{v0}$, the radial optical depth, integrated through the disc mid-plane ($\Theta = \pi/2$), at frequency $v_0 = c/\lambda_0$ with $\lambda_0 = 0.55$ µm,

$$n_0 = \frac{\tau_{v0}}{C_{v0}^{\text{ext}} r_d \ln\left(\frac{R_{\text{out}}}{R_{\text{in}}}\right)}. \quad (25)$$

The opacities were taken from Draine & Lee (1984). They are associated with a unique dust species composed of spherical astronomical silicate grains, with a radius of 0.12 µm and a density of $3.6 \text{ g cm}^{-3}$. A table of pre-computed values for $C_v^{\text{ext}}$ and $C_v^{\text{sca}}$ is available in Pascucci et al. (2004). We considered the optically thin and thick cases $\tau_{v0} = 10^{-1}$ and $10^2$, respectively. The benchmarks were produced with RADMC-3D<fnref n="3" /> (Dullemond et al. 2012), a Monte-Carlo radiative transfer code and whose previous version is part of the original benchmarks from Pascucci et al. (2004). For our code, we used the same spatial and frequency grids as in the previous test. For RADMC-3D, we used 128 points both in radial and angular direction, with a logarithmic sampling in the radial direction.

The temperature of the disc is displayed in Fig. 5. The DGFEM code successfully reproduces it. The regions of the disc with steep gradients are always well reproduced, even with a fairly reasonable number of nodes. This result is a direct consequence of having a high-order numerical scheme. The radiation field can exhibit discontinuities, because of boundary conditions

<footnote n="2">² Available at http://faculty.washington.edu/ivezic/dusty_web/</footnote>

<footnote n="3">³ Available at https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Temperature profiles for the axis-symmetric envelope with $\tau_{v0} = 10^{-1}$ (blue curve) and $\tau_{v0} = 10^{2}$ (orange curve), in the in the disc mid-plane (left panel) and at $r = 2$ AU (right panel). The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** SED profiles for the axis-symmetric envelope with $\tau_{v0} = 10^{-1}$ (blue curve) and $\tau_{v0} = 10^{2}$ (orange curve). The left and right panels correspond to $i = 12.5$ and 77.5 deg, respectively. The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes.</figcaption></figure>

or very strong density gradients. Our numerical tests revealed little to no oscillations and very few negative values for the specific intensity, which did not pollute the computation of the mean radiation field. On average, the absolute relative differences stay below one percent for both test cases. The temperature in the disc mid-plane is very well reproduced, highlighting that the method is able to correctly represent the shadowed regions of the disc (the cold outer mid-plane regions shadowed by the inner dense regions).

In Fig. 6, we show the corresponding emerging SEDs, for two inclinations with respect to the polar axis $i = 12.5$ and 77.5 deg. Again our SEDs were computed from the ray-tracing of the emissivity $\eta_v$. They consist of a stellar and an envelope contribution (see Eq. (C.1)). The stellar component that peaks at around $\lambda \approx 0.6$ µm was computed from of Eq. (C.2) and the discrepancies are always < 1% where this part dominates.

Concerning the contribution from the envelope, both the scattering ($<$ 1 µm) and emission ($>$10 µm) parts agree well. On average, the absolute relative differences stay below 5%, which are in the typical discrepancy levels of the original benchmark. We note a peak in the discrepancies for the optically thick case, between 8.7 and 10.2 µm (according to the resolution of our frequency grid). The same behaviour was previously observed in Pascucci et al. (2004), between a grid-based and a Monte-Carlo code. The authors suggested that, at these wavelengths, the flux mainly comes from the inner disc regions (between 1 and 2 AU) and the numerical simulations are particularly sensitive to the resolution of the inner parts. However, we tried to increase the resolution in these regions, which did not result in any improvement.

In Fig. 7, we show a set of images from the DGFEM code of the 10 AU disc inner regions, at $\lambda = 2.3, 4.5, \text{and } 12.1$ µm
<figure><img src="image_8.png" /><figcaption>Fig. 7. Images at $\lambda = 2.3, 4.5$, and $12.1 \mu m$ (left, middle, and right panels, respectively) of the 10 AU inner regions of the axis-symmetric envelope (Sect. 5.2) and computed from the DGFEM solution with the ray-tracing module (Appendix C). The top, middle and bottom panels correspond to the inclinations $i = 12.5, 77.5$, and $90$ deg, respectively. The colour code shows the specific intensity value (in W m$^{-2}$ Hz$^{-1}$ sr$^{-1}$) of the envelope $I_{\nu}^{\text{env}}$, inside each pixel. The solid white lines show the iso-contours of $I_{\nu}^{\text{env}}$.</figcaption></figure>

and for several inclinations $i = 12.5, 77.5$, and $90$ deg. These wavelengths are characteristic of the operating spectral bands of the interferometric instruments, such as GRAVITY (Gravity Collaboration 2017) and MATISSE (Lopez et al. 2022). On average, the agreement between the images from the DGFEM code and the images from RADMC-3D is around 10%, for all frequencies and inclinations. We show in Fig. 8 the comparison of an image slice at $x = 0$ AU, for two inclinations $i = 12.5$ and $77.5$ deg and for $\lambda = 2.3, 4.5$, and $12.1 \mu m$. In general, the disc emitting inner regions (peaks in Fig. 7) are reproduced very well ($\epsilon(I_{\nu}^{\text{env}}) \le 3\%$). The biggest discrepancies occur in the wings of the peaks, where the gradient of intensity is the steepest (in logarithmic scaling).

### 5.3. Effect of the grid resolution on temperatures and SEDs

As we previously mentioned in Sects. 5.1 and 5.2, we used a ray-tracing procedure to accurately compute SEDs. However, in practice, we already know the intensity $I_{\nu}(r = R_{\text{out}}, \Theta, \mu, \varphi)$ from the DGFEM solution at the outer edge of the envelope. This intensity can theoretically be used to compute the SEDs directly with no further post-processing treatment required. It is, however, not advised to proceed in this way, at least in multi-dimensional cases where the radiation field has a complex

geometry, because of the poor accuracy of the local solution on the outer edge due to the poor resolution of the grid ($R_{\text{out}} \Delta\Theta$ is large). This effect becomes even more important if we are interested in studying the inner parts of the disc and hence probing a very small part of the outer shell of the envelope. Consequently a post-processing ray-tracing step is usually needed in order to capture all the features of the disc accurately. The post-processing ray-tracing step requires a good estimate for the emissivity (or source function, see Eq. (11) in the whole envelope to compute accurate SEDs (see Appendix C). The emissivity is implicitly a function of the temperature of the medium (through the black-body emission $B_{\nu}(T)$ and the radiative equilibrium Eq. (3)), and hence an accurate temperature is needed to compute SEDs with precision.

In order to study the effect of the resolution of the computational grid on the temperature agreement, we computed the $L^2$ norm of the relative residual of the temperature, between our code $T_{\text{DG}}$ and the RADMC-3D code $T_{\text{ref}}$. It is defined as,

$$L^2 = \left( \frac{\int_{\nu} \left( \frac{T_{\text{DG}} - T_{\text{ref}}}{T_{\text{ref}}} \right)^2 r^2 \sin \Theta dr d\Theta}{\int_{\nu} r^2 \sin \Theta dr d\Theta} \right)^{\frac{1}{2}}, \quad (26)$$
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Image slices at *x* = 0 AU, for two inclinations *i* = 12.5 and 77.5 deg (left and right panels, respectively). The blue, orange, and green curves are the RADMC-3D intensities at *λ* = 2.3, 4.5, and 12.1 µm, respectively. The cross marks represent the solution from this study. The lower panels show the relative differences between the two codes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** $L^2$ norm of the relative differences in temperature between this study and RADMC-3D versus the number of elements in the computational grid ($N = N_r \times N_\Theta \times N_\mu \times N_\phi$). The $L^2$ norm was computed for the benchmark problem presented in Sect. 5.2.</figcaption></figure>

where the integral runs across the volume $V$ of the envelope. In Fig. 9, we display the $L^2$ norm versus the total number of elements of our grid (we would like to emphasise that the total number of elements is $N = N_r \times N_\Theta \times N_\mu \times N_\phi$ and that each element has 3 × 2 × 3 × 3 nodes). We varied $N$ in such a way as to keep the number of elements along each dimension identical ($N_r = N_\Theta = N_\mu = N_\phi$).

We can see that the $L^2$ norm saturates when $N \approx 10^4$, which corresponds to a number of elements along each dimension of ten. After that point, no substantial improvement in the temperature agreement is obtained. The computational grid at this point is still quite coarse, compared to the commonly used methods in the literature such as finite differences (Steinacker et al. 2003),

short-characteristics, (Dullemond & Turolla 2000) or Monte-Carlo (Wolf 2003) methods. This is one of the advantages of the DGFEM that can achieve a given error threshold with less resolution, thanks to the form of the intensity inside each element (see Eq. (13)).

In Fig. 10, we compare the SEDs obtained either from the post-processing ray-tracing or directly by using $I_v(r = R_{\text{out}}, \Theta, \mu, \phi)$. We compare these SEDs with the benchmark curves we already presented in Sects. 5.1 and 5.2 ($\tau = 100$). While for the spherically symmetric case, they agree quite well, $I_v(r = R_{\text{out}}, \Theta, \mu, \phi)$ fails to give an accurate SED for the 2D disc configuration. This is expected since for this case, the radiation field has a dependence on the $\Theta$ and $\phi$ variables. We note that the increase in the grid size should make the SED from $I_v(r = R_{\text{out}}, \Theta, \mu, \phi)$ converge to the ray-tracing one, but we could not investigate this hypothesis due to the limits of our current computational resources. In the absence of a finer grid, the post-treatment ray-tracing is thus needed if one wants to compute the emerging fluxes from the disc with precision.

## 6. Conclusions

For this study, we applied the DGFEM to the frequency-dependent 2D radiative transfer equation, coupled with the radiative equilibrium equation, inside axis-symmetric circumstellar envelopes. We have shown that it can accurately compute the temperature field and allow for the correct determination of images and SED profiles, via ray-tracing techniques. A desirable feature of the method is the ability to control the order of the numerical scheme via the number of nodes inside each element, meaning that a high-order numerical scheme can simply be achieved by increasing the number of nodes inside each element.

The DGFEM formulation, Eq. (18), is particularly adapted to parallelisation and adaptive mesh refinement (AMR) grids (e.g. Frisken & Perry 2002). In this regard, direct access to the residual, Eq. (16), might provide a robust estimate for the error that could be used as a criterion for the grid refinement (see also Richling et al. 2001). We also note that a further 3D generalisation is also possible and straightforward. Together with Monte-Carlo and ray-tracing techniques, the DGFEM provides
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Emerging SEDs of RADMC-3D (solid black lines) and of the DGFEM code, computed from the post-processing ray-tracing procedure (orange dots), and from the DGFEM solution at the outer edge (blue crosses). The left panel corresponds to the optically thick spherically symmetric case ($\tau = 100$, see Sect. 5.1) while the right panel corresponds to the optically thick disc benchmark from Sect. 5.2. The lower panels show the relative difference between the curves of this study and the benchmarks.</figcaption></figure>

an additional method for solving the radiative transfer problem, and it could be used in cases where the other methods are expected to be less efficient.

*Acknowledgements.* This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3-year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côte-d'Azur (UCA). The 1D benchmark profiles were computed with the radiative transfer code DUSTY, available at http://faculty.washington.edu/ivezic/dusty_web/. The 2D benchmark were computed with RADMC-3D: https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/ Part of the computations were carried out with the help of OPAL-Meso computing facilities. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d'Azur (CRIMSON) for providing resources and support. We are grateful to the referee and the editor whose insights helped improve this work. We also want to express our gratitude to the language editor who helped improve the clarity of this paper.

## References

* Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering by a Sphere (John Wiley & Sons, Ltd), 82

* Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

* Cockburn, B. 2003, Z. Angew. Math. Mech., 83, 731

* Cui, X., & Li, B. Q. 2005, JQSRT, 96, 383

* Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

* Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187

* Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: a multi-purpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015]

* Dykema, P. G., Klein, R. I., & Castor, J. I. 1996, ApJ, 457, 892

* Frisken, S., & Perry, R. 2002, J. Graph. Tools, 7, 1

* Gravity Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

* Hesthaven, J. S., & Warburton, T. 2007, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (Springer Science & Business Media), 19

* Ivezic, Z., Groenewegen, M. A. T., Men'shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121

* Karniadakis, G. E., & Kirby II, R. M. 2003, Parallel Scientific Computing in C and MPI: A Seamless Approach to Parallel Algorithms and their Implementation (Cambridge University Press)

* Kitzmann, D., Bolte, J., & Patzer, A. B. C. 2016, A&A, 595, A90

* Kopriva, D. A., & Gassner, G. 2010, J. Sci. Comput., 44, 136

* Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7

* Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67

* Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

* Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192

* Machorro, E. 2007, J. Computat. Phys., 223, 67

* Mihalas, D., & Weibel-Mihalas, B. 1999, Foundations of Radiation Hydrodynamics, Dover Books on Physics (Dover Publications)

* Nair, R. D., Levy, M. N., & Lauritzen, P. H. 2011, Emerging Numerical Methods for Atmospheric Modeling (Berlin, Heidelberg: Springer Berlin Heidelberg), 251

* Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431

* Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793

* Perdigon, J., Niccolini, G., & Faurobert, M. 2021, A&A, 653, A139

* Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

* Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press)

* Reed, W. H., & Hill, T. R. 1973, Triangular mesh methods for the neutron transport equation, Tech. rep. (North Mexico, USA:Los Alamos Scientific Laboratory)

* Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776

* Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405

* Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63

* Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383

* Wolf, S. 2003, Comput. Phys. Commun., 150, 99
<figure><img src="image_12.png" /><figcaption>Fig. A.1. Geometry of a ray in the inner cavity. Point B is the opposite point of A, in the direction $\Omega$.</figcaption></figure>

## Appendix A: Boundary conditions on a spherical enclosed cavity

The boundary condition in Eq. (7) requires to know, for a given point A ($R_{in}$, $\Theta$) on the inner cavity of radius $R_{in}$ and in a direction $\Omega$ ($\mu$, $\varphi$) ($\Omega, \hat{r} > 0$), the coordinates of the opposite point B ($R_{in}$, $\Theta'$, $\mu' = \cos \theta'$, $\varphi'$). In the following, the coordinates with a prime superscript are associated with the opposite point B. The problem is illustrated in Fig. A.1, where we show the plane containing the point A, B and the star. The problem is axisymmetric around the polar axis and the radiation field does not depend on the azimuthal angle $\Phi$ (see Fig. 1). Consequently, we can arbitrarily choose the $\Phi$ coordinate of the point A and set it to be $\Phi = \pi/2$. This convention simplifies the computations.

First, we see that $\theta' = \pi - \theta$ because the OAB triangle is isosceles, we then have,

$$ \mu' = -\mu. \qquad (A.1) $$

The position $r_B$ of the opposite point B, is linked to the position $r_A$ of the point A through the relation

$$ r_B = r_A - s_{AB} \Omega, \qquad (A.2) $$

with $s_{AB}$ the distance between the point A and B. This length is equal to $s_{AB} = 2\mu R_{in}$, because it corresponds to the base of the isosceles triangle OAB, with $OA = OB = R_{in}$. The direction vector, $\Omega$, is,

$$ \begin{aligned} \Omega &= \left[ \begin{matrix} \mu \\ \frac{\sqrt{1-\mu^2}}{\sqrt{1-\mu^2}} \cos \varphi \\ \frac{\sqrt{1-\mu^2}}{\sqrt{1-\mu^2}} \sin \varphi \end{matrix} \right]_{(\hat{x}, \hat{\theta}, \hat{\phi})} \\ &= \left[ \begin{matrix} -\sqrt{1-\mu^2} \sin \varphi \\ \sin \Theta \mu + \cos \Theta \sqrt{1-\mu^2} \cos \varphi \\ \cos \Theta \mu - \sin \Theta \sqrt{1-\mu^2} \cos \varphi \end{matrix} \right]_{(\hat{x}, \hat{y}, \hat{z})}. \end{aligned} \qquad (A.3) $$

Eq. (A.2) is then

$$ \begin{aligned} r_B &= R_{in} \begin{bmatrix} \sin \Theta' \cos \Phi' \\ \sin \Theta' \sin \Phi' \\ \cos \Theta' \end{bmatrix}_{(\hat{x}, \hat{y}, \hat{z})} \\ &= R_{in} \begin{bmatrix} 2\mu \sqrt{1-\mu^2} \sin \varphi \\ \sin \Theta (1-2\mu^2) - 2 \cos \Theta \mu \sqrt{1-\mu^2} \cos \varphi \\ \cos \Theta (1-2\mu^2) + 2 \sin \Theta \mu \sqrt{1-\mu^2} \cos \varphi \end{bmatrix}_{(\hat{x}, \hat{y}, \hat{z})}. \end{aligned} \qquad (A.4) $$

From the z-component of Eq. (A.4) we get,

$$ \cos \Theta' = \cos \Theta (1 - 2\mu^2) + 2 \sin \Theta \mu \sqrt{1 - \mu^2} \cos \varphi, \qquad (A.5) $$

<figure><img src="image_13.png" /><figcaption>Fig. B.1. Example of an element $D^{i,j,k,l}$ (dashed lines) using $n_a = n_c = n_d = 3$ and $n_b = 2$ points. The left and right panels correspond to 2D slices in the $(r, \Theta)$ and $(\mu, \varphi)$ planes, respectively (we display here the local coordinate system given by Eq. B.2). The black dots correspond to the nodes where the solution is computed while the crosses correspond to the interpolated value of I used to compute the numerical flux $F_h^*$ at the interface, along the $\Theta$ coordinate.</figcaption></figure>

and consequently $\Theta' = \arccos|\cos \Theta'|$. The absolute value occurs if we consider the planar symmetry with respect to the equatorial plane.

The angle $\varphi'$ at the point B verifies,

$$ \tan \varphi' = \frac{\Omega \cdot \hat{\Phi}'}{\Omega \cdot \hat{\Theta}'}, \qquad (A.6) $$

with,

$$ \hat{\Phi}' = \begin{bmatrix} -\sin \Phi' \\ \cos \Phi' \\ 0 \end{bmatrix}_{(\hat{x}, \hat{y}, \hat{z})}, \quad \hat{\Theta}' = \begin{bmatrix} \cos \Theta' \cos \Phi' \\ \cos \Theta' \sin \Phi' \\ -\sin \Theta' \end{bmatrix}_{(\hat{x}, \hat{y}, \hat{z})}, \qquad (A.7) $$

and $\Omega$ given by Eq. (A.3). Then, using Eq. (A.4) in Eq. (A.6) yields, after tedious calculations,

$$ \tan \varphi' = \frac{\sin \varphi \sin \Theta}{(1 - 2\mu^2) \cos \varphi \sin \Theta - 2\mu \sqrt{1 - \mu^2} \cos \Theta}. \qquad (A.8) $$

Consequently $\varphi' = \arctan|u|/v$ with $u$ and $v$ being the numerator and denominator in Eq. (A.8), respectively. This time, the absolute value is present because the symmetry around the polar axis causes the $\varphi'$ angle to take its values in $[0, \pi]$. We note that the previous formula is not valid for the couples of coordinates $(\mu, \varphi) = (\sin(\Theta/2), 0)$ and $(\mu, \varphi) = (\cos(\Theta/2), \pi)$, because it corresponds to have the opposite point B on the poles, where $\hat{\Theta}'$ and $\hat{\Phi}'$ are not defined. These particular cases are in practice not a problem since the radiation field is independent of $\varphi'$ at these points.

## Appendix B: DGFEM calculations

The computation of the terms in Eq. (18) requires the choice of a particular quadrature in the cell $D^{i,j,k,l}$, for each coordinate $(r, \Theta, \mu, \varphi)$. For the $r, \mu, \varphi$ coordinates, we make use of the Gauss-Lobatto quadrature, which has the advantage of having roots at the end points of the cell. This avoids the use of an interpolation formula when computing the numerical flux $F_h^*$ at the element edges. The Gauss-Lobatto quadrature can exactly integrate polynomials up to degree $2n-3$, with $n$ the number of nodes. For the coordinate $\Theta$, we cannot use the Gauss-Lobatto quadrature, at least in the cells that are touching $\Theta = 0$, because the radiative transfer equation in the spherical coordinates system is not defined at the pole. For this reason and to keep
an homogeneous method along $\Theta$, we use a Gauss-Legendre quadrature for this coordinate. The Gauss-Legendre quadrature is more precise and can exactly integrate polynomials up to degree $2n - 1$, but at the expense of an interpolation method required to compute the flux at the cell edges. An example of an element $D^{i,j,k,l}$ with the chosen nodes is shown in Fig. B.1.

In the following, all the superscript indexes refer to the element identification while the subscripts denote each node in the considered element. We start with the volume term in Eq. (18),

$$ 
\begin{aligned}
& \int_{D^{i,j,k,l}} (\kappa^{\text{ext}} \tilde{I}_h - \tilde{\eta}) h_{a',b',c',d'} d^4 \mathbf{x} \\
&= \frac{\Delta x^{i,j,k,l}}{16} \int_{-1}^{1} (\kappa^{\text{ext}} \tilde{I}_h^{i,j,k,l} - \tilde{\eta}^{i,j}) h_{a',b',c',d'} d^4 \tilde{\mathbf{x}} \\
&= \frac{\Delta x^{i,j,k,l}}{16} \sum_{a,b,c,d} W_{a,b,c,d} (\kappa_{a,b}^{\text{ext},i,j} \tilde{I}_{a,b,c,d}^{i,j,k,l} - \tilde{\eta}_{a,b}^{i,j}) h_{a',b',c',d'} (\tilde{\mathbf{x}}_{a,b,c,d}).
\end{aligned}
\quad \text{(B.1)}
$$

We note that $\Delta x^{i,j,k,l} = \Delta r^i \Delta \Theta^j \Delta \mu^k \Delta \varphi^l$ is the 4D volume of the element $D^{i,j,k,l}$. For integration, we rather use the local coordinates $(\tilde{r}, \tilde{\Theta}, \tilde{\mu}, \tilde{\varphi})$, defined as (e.g for the r coordinate),

$$ \tilde{r} = \frac{2}{\Delta r^i} (r - r^{i+1/2}), \quad \text{(B.2)} $$

with $\Delta r^i$, the element width along the coordinate r and $r^{i+1/2}$, the radial coordinate of the centre of the element. The same expression holds for the other coordinates. The quantities $W_{a,b,c,d} = W_r(\tilde{r}_a) W_\Theta(\tilde{\Theta}_b) W_\mu(\tilde{\mu}_c) W_\varphi(\tilde{\varphi}_d)$ are the weights associated with the different quadrature in each direction. Finally, $h_{a',b',c',d'}(\tilde{\mathbf{x}}_{a,b,c,d})$ is the 4D Lagrange polynomials, defined in Eq. (14), evaluated at the node $\tilde{\mathbf{x}}_{a,b,c,d} = (\tilde{r}_a, \tilde{\Theta}_b, \tilde{\mu}_c, \tilde{\varphi}_d)$. By definition of the Lagrange polynomials, we have,

$$ h_{a',b',c',d'}(\tilde{r}_a, \tilde{\Theta}_b, \tilde{\mu}_c, \tilde{\varphi}_d) = \delta_{a',a} \delta_{b',b} \delta_{c',c} \delta_{d',d}, \quad \text{(B.3)} $$

with $\delta_{a',a}$ the usual delta Kronecker. The other volume term in Eq. (18) is expressed as,

$$ 
\begin{aligned}
& \int_{D^{i,j,k,l}} F_h \cdot \nabla_x h_{a',b',c',d'} d^4 \mathbf{x} \\
&= \int_{D^{i,j,k,l}} a_r \tilde{I}_h \partial_r h_{a',b',c',d'} d^4 \mathbf{x} + \int_{D^{i,j,k,l}} a_\Theta \tilde{I}_h \partial_\Theta h_{a',b',c',d'} d^4 \mathbf{x} \\
&\quad + \int_{D^{i,j,k,l}} a_\mu \tilde{I}_h \partial_\mu h_{a',b',c',d'} d^4 \mathbf{x} + \int_{D^{i,j,k,l}} a_\varphi \tilde{I}_h \partial_\varphi h_{a',b',c',d'} d^4 \mathbf{x} \\
&= \frac{\Delta x^{i,j,k,l}}{8} \sum_{a,b,c,d} W_{a,b,c,d} \tilde{I}_{a,b,c,d}^{i,j,k,l} \left( \frac{a_{\tilde{r},a,b,c,d}^{i,j,k,l} \partial_{\tilde{r}} h_{a',b',c',d'} |\tilde{\mathbf{x}}_{a,b,c,d}}{\Delta r^i} + \frac{a_{\tilde{\Theta},a,b,c,d}^{i,j,k,l} \partial_{\tilde{\Theta}} h_{a',b',c',d'} |\tilde{\mathbf{x}}_{a,b,c,d}}{\Delta \Theta^j} + \frac{a_{\tilde{\mu},a,b,c,d}^{i,j,k,l} \partial_{\tilde{\mu}} h_{a',b',c',d'} |\tilde{\mathbf{x}}_{a,b,c,d}}{\Delta \mu^k} + \frac{a_{\tilde{\varphi},a,b,c,d}^{i,j,k,l} \partial_{\tilde{\varphi}} h_{a',b',c',d'} |\tilde{\mathbf{x}}_{a,b,c,d}}{\Delta \varphi^l} \right).
\end{aligned}
\quad \text{(B.4)}
$$

In Eq. (B.4), we used the definition of the flux Eq. (12). The quantity $\partial_{\tilde{r}} h_{a',b',c',d'} |_{x_{a,b,c,d}}$ is the partial derivative of the Lagrange

polynomial, with respect to the coordinate r, evaluated at the node $\tilde{\mathbf{x}}_{a,b,c,d}$ (with similar definitions for the other coordinates).

The last term to evaluate is the surface integral (first term in Eq. 18). The 4D element is delimited by $2 \times 4 = 8$ surfaces. We give here the derivation for the surface integrals normal to the coordinate r and the other terms will follow by substitution of the indices. We have $\hat{s} = \pm \hat{r}$ for the radial right and left surfaces ($\tilde{r} = \pm 1$), respectively,

$$ 
\begin{aligned}
& \oint_{\partial D^{i,j,k,l}} [F_r^* h_{a',b',c',d'}]_{\tilde{r}=-1}^{\tilde{r}=1} d\Theta d\mu d\varphi \\
&= \frac{\Delta x^{i,j,k,l}}{8 \Delta r^i} W_{b',c',d'} [F_{\tilde{r}}^* (\tilde{r}, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}) h_{a'}(\tilde{r})]_{\tilde{r}=-1}^{\tilde{r}=1}.
\end{aligned}
\quad \text{(B.5)}
$$

At the cell boundaries, we use the upwind numerical flux, for example, at the right edge,

$$ 
\begin{aligned}
F_{\tilde{r}}^*(\tilde{r} = 1, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}) = & \\
& \max \left\{ a_{\tilde{r}}^{i,j,k,l} (1, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}), 0 \right\} \tilde{I}_{h}^{i,j,k,l} (1, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}) \\
& + \min \left\{ a_{\tilde{r}}^{i,j,k,l} (1, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}), 0 \right\} \tilde{I}_{h}^{i+1,j,k,l} (-1, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}).
\end{aligned}
\quad \text{(B.6)}
$$

The use of the Gauss-Lobatto quadrature further simplify the computation because the quadrature nodes include the end-points $\tilde{r} = \pm 1$. The numerical flux can then be expressed as

$$ 
\begin{aligned}
& [F_{\tilde{r}}^*(\tilde{r}, \tilde{\Theta}_{b'}, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}) h_{a'}(\tilde{r})]_{\tilde{r}=-1}^{\tilde{r}=1} = \\
& \delta_{a',n_a-1} \max \left\{ a_{\tilde{r}_{n_a-1,1,b',c',d'}}^{i,j,k,l}, 0 \right\} \tilde{I}_{n_a-1,1,b',c',d'}^{i,j,k,l} \\
& + \delta_{a',n_a-1} \min \left\{ a_{\tilde{r}_{n_a-1,1,b',c',d'}}^{i,j,k,l}, 0 \right\} \tilde{I}_{0,b',c',d'}^{i+1,j,k,l} \\
& - \delta_{a',0} \max \left\{ a_{\tilde{r}_{0,b',c',d'}}^{i,j,k,l}, 0 \right\} \tilde{I}_{n_a-1,1,b',c',d'}^{i-1,j,k,l} \\
& - \delta_{a',0} \min \left\{ a_{\tilde{r}_{0,b',c',d'}}^{i,j,k,l}, 0 \right\} \tilde{I}_{0,b',c',d'}^{i,j,k,l}.
\end{aligned}
\quad \text{(B.7)}
$$

The same form holds for the surface integrals normal to the $\mu$ and $\varphi$ coordinates. We however note that we have $a_\mu \ge 0$, $a_\varphi \le 0$ $\forall \mathbf{x} \in D$, so there are no terms proportional to $\tilde{I}_{a',b,c',d'}^{i,j,k+1,l}$ and $\tilde{I}_{a',b,c',d'}^{i,j,k,l-1}$. For the flux computation normal to the $\Theta$ coordinate, we do not directly have the value of the solution at the element interface (see Fig. B.1), we need to interpolate the solution with the help of Eq. (13),

$$ 
\begin{aligned}
F_{\tilde{\Theta}}^*(\tilde{r}_{a'}, \tilde{\Theta} = 1, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}) = & \\
& \max \left\{ a_{\tilde{\Theta}}^{i,j,k,l} (\tilde{r}_{a'}, 1, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}), 0 \right\} \sum_b \tilde{I}_{a',b,c',d'}^{i,j,k,l} h_b (1) \\
& + \min \left\{ a_{\tilde{\Theta}}^{i,j,k,l} (\tilde{r}_{a'}, 1, \tilde{\mu}_{c'}, \tilde{\varphi}_{d'}), 0 \right\} \sum_b \tilde{I}_{a',b,c',d'}^{i,j,k+1,l} h_b (-1).
\end{aligned}
\quad \text{(B.8)}
$$

All the terms in this section can be put in the form of the system of equations, given by Eq. (22),

$$ \mathcal{A}^{i,j,k,l} \tilde{I}_h^{i,j,k,l} = \mathbf{b}^{i,j,k,l}. \quad \text{(B.9)} $$

In Eq. (B.9), we replaced the RHS by $\mathbf{b}^{i,j,k,l}$ because we do not need to formally write the non-diagonal matrices $\mathcal{A}^{i\pm 1, j\pm 1, k-1, l+1}$. To assemble $\mathcal{A}^{i,j,k,l}$, we make use of the global index $\alpha =$
$anbnncn_d + bnncn_d + cn_d + d$. The elements of $\mathcal{A}^{i,j,k,l}$ and $\mathbf{b}^{i,j,k,l}$ are then,

$$
\begin{equation}
\begin{split}
\mathcal{A}_{\alpha'\alpha}^{i,j,k,l} ={}& W_{a,b,c,d} \kappa_{a,b}^{\text{ext},i,j} \delta_{a',a} \delta_{b',b} \delta_{c',c} \delta_{d',d} \\
& + \frac{2}{\Delta r^l} \delta_{b',b} \delta_{c',c} \delta_{d',d} W_{b,c,d} \left( \max \left\{ a_{\bar{r}_{a,b,c,d},l}^{i,j,k,l}, 0 \right\} \delta_{a',n_a-1} \delta_{a,n_a-1} \right. \\
& \qquad \left. - \min \left\{ a_{\bar{r}_{a,b,c,d},l}^{i,j,k,l}, 0 \right\} \delta_{a',0} \delta_{a,0} - W_a a_{\bar{r}_{a,b,c,d},l}^{i,j,k,l} \partial_{\bar{r}} h_{a'} |_{\bar{r}_a} \right) \\
& + \frac{2}{\Delta \Theta^j} \delta_{a',a} \delta_{c',c} \delta_{d',d} W_{a,c,d} \left( \max \left\{ a_{\bar{\Theta}_{a,c,d},l}^{i,j,k,l} |_{\bar{\Theta}=1}, 0 \right\} h_{b'} |_{\bar{\Theta}=1} h_b |_{\bar{\Theta}=1} \right. \\
& \qquad \left. - \min \left\{ a_{\bar{\Theta}_{a,c,d},l}^{i,j,k,l} |_{\bar{\Theta}=-1}, 0 \right\} h_{b'} |_{\bar{\Theta}=-1} h_b |_{\bar{\Theta}=-1} - W_b a_{\bar{\Theta}_{a,b,c,d},l}^{i,j,k,l} \partial_{\bar{\Theta}} h_{b'} |_{\bar{\Theta}_b} \right) \\
& + \frac{2}{\Delta \mu^k} \delta_{a',a} \delta_{b',b} \delta_{d',d} W_{a,b,d} a_{\bar{\mu}_{a,b,c,d},l}^{i,j,k,l} \left( \delta_{c',n_c-1} \delta_{c,n_c-1} - W_c \partial_{\bar{\mu}} h_{c'} |_{\bar{\mu}_c} \right) \\
& \qquad \left. - \frac{2}{\Delta \varphi^l} \delta_{a',a} \delta_{b',b} \delta_{c',c} W_{a,b,c} a_{\bar{\varphi}_{a,b,c,d},l}^{i,j,k,l} \left( \delta_{d',0} \delta_{d,0} + W_d \partial_{\bar{\varphi}} h_{d'} |_{\bar{\varphi}_d} \right).
\end{split}
\tag{B.10}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
b_{\alpha'}^{i,j,k,l} ={}& W_{a',b',c',d'} \tilde{\eta}_{a',b'}^{i,j} \\
& + \frac{2}{\Delta r^i} W_{b',c',d'} \left( \max \left\{ a_{\bar{r}_{n_a-1,b',c',d'},0}^{i-1,j,k,l} \right\} \tilde{I}_{n_a-1,b',c',d'}^{i-1,j,k,l} \delta_{a',0} \right. \\
& \qquad \left. - \min \left\{ a_{\bar{r}_{0,b',c',d'},0}^{i+1,j,k,l} \right\} \tilde{I}_{0,b',c',d'}^{i+1,j,k,l} \delta_{a',n_a-1} \right) \\
& + \frac{2}{\Delta \Theta^j} W_{a',c',d'} \left( \max \left\{ a_{\bar{\Theta}_{a',c',d'},0}^{i,j-1,k,l} \right\} h_{b'} |_{\bar{\Theta}=1} h_b |_{\bar{\Theta}=1} \right. \\
& \qquad \left. - \min \left\{ a_{\bar{\Theta}_{a',c',d'},0}^{i,j+1,k,l} \right\} h_{b'} |_{\bar{\Theta}=1} h_b |_{\bar{\Theta}=-1} \right) \\
& + \frac{2}{\Delta \mu^k} W_{a',b',d'} a_{\bar{\mu}_{a',b',n_c-1,d'},0}^{i,j,k-1,l} \tilde{I}_{a',b',n_c-1,d'}^{i,j,k-1,l} \delta_{c',0} \\
& \qquad \left. - \frac{2}{\Delta \varphi^l} W_{a',b',c'} a_{\bar{\varphi}_{a',b',c',0}}^{i,j,k,l+1} \tilde{I}_{a',b',c',0}^{i,j,k,l+1} \delta_{d',n_d-1} \right).
\end{split}
\tag{B.11}
\end{equation}
$$

Appendix C: Ray-tracing module

The SEDs and intensity maps from the DGFEM code, shown in
Sect 5, were computed with the help of the ray-tracing routine
we present here. This procedure is quite generally used in the
literature, for all type of codes, and is not a limitation of the
DGFEM method itself (see e.g Pinte et al. 2009, Sect 2.2.3)

For the SED, we need to estimate the total flux $f_v^{\text{obs}}$ that an
observer receive from the object, situated at a distance $d \gg R_{\text{out}}$
and doing an angle $i$ with the polar axis (see Fig. C.1). Because
we decoupled the stellar from the envelope radiation (see Eq. 4),
the total flux is made of the stellar and envelope total flux,

$$
f_v^{\text{obs}}(i, d) = f_v^{\text{obs},*} + f_v^{\text{obs env}}. \tag{C.1}
$$

If we assume the star to be an unresolved black-body point
source, the stellar flux at distance *d* is the flux of the star at the
stellar surface attenuated by the circumstellar matter present in
the direction of the line of sight (black dotted line in Fig. C.1),
and with the dilution factor (*R*<sub>⋆</sub>/*d*<sup>2</sup>,

$$
\begin{align}
f_v^{\text{obs},*}(i, d) &= \pi \left(\frac{R_*}{d}\right)^2 B_v(T_*) \exp\{-\tau(i)\} \nonumber \\
&\quad \text{with } \tau(i) = \int_{R_{\text{in}}}^{R_{\text{out}}} \kappa_v^{\text{ext}}(r, i) \, dr \tag{C.2}
\end{align}
$$

<figure><img src="image_14.png" /><figcaption>**Fig. C.1.** Example of a ray (red line) normal to the image plane, crossing the spherical grid. In this example, we display a square image of size $L \times L$ and a circular image of radius $R$. The red crosses represent the intersections between the ray and the grid. The values of $\kappa_v^{\text{ext}}$ and $S_v$ at these intersections are linearly interpolated from the grid adjacent values (red dots).</figcaption></figure>

To compute the envelope flux and intensity maps, we define
an image plane (x̂, ŷ), at a distance d ≫ R<sub>out</sub> from the star
and tilted with an angle i with respect to the polar axis. The x
and y axes are oriented with the help of the spherical coordi-
nates system (r̂, Θ̂, Φ̂). In this plane, we can construct images
of any geometry but we only consider the special cases of a
square image of width L and a circular image of radius R. Since
d ≫ R<sub>out</sub>, the flux inside the image can formally be written,

$$
\begin{align}
f_v^{\text{obs,env}}(i) &= \int_{-\frac{L}{2}}^{\frac{L}{2}} \int_{-\frac{L}{2}}^{\frac{L}{2}} I_v^{\text{env}}(x, y, \hat{r}) \frac{dx dy}{d^2}, \tag{C.3} \\
\text{or } f_v^{\text{obs,env}}(i) &= \int_0^{2\pi} \int_0^R I_v^{\text{env}}(r, \omega, \hat{r}) \frac{r \, d\omega dr}{d^2}. \nonumber
\end{align}
$$

In practice, images are made of a collection of pixels in which
we evaluate the emerging specific intensity at the pixel centre. A
square (circular) image, is divided into $N \times N$ ($N_r \times N_\omega$) pixels
and the flux in the image can then be rewritten,

$$
\begin{align}
f_v^{\text{obs,env}}(i) &\approx \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} I_v^{\text{env}}(x_i, y_j, \hat{r}) \frac{\Delta x_i \Delta y_j}{d^2}, \tag{C.4} \\
\text{or } f_v^{\text{obs,env}}(i) &\approx \sum_{i=0}^{N_r} \sum_{j=0}^{N_\omega} I_v^{\text{env}}(r_i, \omega_j, \hat{r}) \frac{r_i \Delta \omega_j \Delta r_i}{d^2}, \nonumber
\end{align}
$$

with $\Delta x_i \Delta y_j (r_i \Delta \omega_j \Delta r_i)$ the pixel size of the square (circular) image. We note that the circular image is particularly well-suited for the computation of $f_v^{\text{obs,env}}(i)$ since we can easily increase the number of pixels in the centre of the image in order to resolve the disc inner parts.
The emerging specific intensity crossing each pixel centre $\mathbf{r}_0 = x_i \hat{x} + y_j \hat{y} + d \hat{r}$ ($\mathbf{r}_0 = r_i \sin \omega_j \hat{x} + r_i \cos \omega_j \hat{y} + d \hat{r}$ for a circular image), along the ray normal to the image plane (red line in Fig. C.1) is computed by integration of the emissivity along the ray,

$$I_{\nu}^{\text{env}}(\mathbf{r}_0, \hat{r}) = \int_{s_{\text{min}}}^{s_{\text{max}}} \eta_{\nu} \exp \{-\tau_{\nu}(s)\} ds, \qquad (C.5)$$

$$\text{with } \tau_{\nu}(s) = \int_{s_{\text{min}}}^{s} \kappa_{\nu}^{\text{ext}} ds'$$

We define $s$ as to be the distance from the pixel centre $\mathbf{r}_0$ to a given point along the ray. The quantities $\eta_{\nu}$ and $\kappa_{\nu}^{\text{ext}}$ are the emissivity and the extinction coefficient, respectively, as defined in Eq. (1). The points $s_{\text{min}}$ and $s_{\text{max}}$ correspond to the two intersections of the ray with the sphere of radius $R_{\text{out}}$.

In practice, $\eta_{\nu}$ and $\kappa_{\nu}^{\text{ext}}$ are defined on a discrete grid and the previous integral can be rewritten as,

$$I_{\nu}^{\text{env}}(\mathbf{r}_0, \hat{r}) = \sum_{i=0}^{n-2} \int_{s_i}^{s_{i+1}} \eta_{\nu} \exp \{-\tau_{\nu}(s)\} ds. \qquad (C.6)$$

The $\{s_i\}$ are the $n$ intersections between the ray and the grid (red crosses in Fig. C.1), with $s_0 = s_{\text{min}}$ and $s_{n-1} = s_{\text{max}}$. The coordinates of all intersections can be computed, in the Cartesian coordinate system. We can express $\Delta I_{\nu}^i$, the contribution to the total intensity $I_{\nu}^{\text{env}}(x, y, \hat{r})$ of each portion between two consecutive intersections as,

$$I_{\nu}^{\text{env}}(\mathbf{r}_0, \hat{r}) = \sum_{i=0}^{n-2} \exp(-\tau_{\nu}^i) \Delta I_{\nu}^i, \\ \text{with } \Delta I_{\nu}^i = \int_{s_i}^{s_{i+1}} \eta_{\nu} \exp\left(-\int_{s_i}^{s} \kappa_{\nu}^{\text{ext}} ds'\right) ds, \qquad (C.7) \\ \text{and } \tau_{\nu}^i = \int_{s_{\text{min}}}^{s_i} \kappa_{\nu}^{\text{ext}} ds'$$

Following Olson et al. (1986), we assume that $\eta_{\nu}$ and $\kappa_{\nu}^{\text{ext}}$ are linear functions between two consecutive intersections. Each contribution $\Delta I_{\nu}^i$ is given by,

$$\Delta I_{\nu}^{i} \approx (1 - \exp\{-\Delta \tau_{\nu}^{i}\} - \beta) S_{\nu}(s_{i}) + \beta S_{\nu}(s_{i+1}), \\ \text{with } \beta = \frac{\Delta \tau_{\nu}^{i} - 1 + \exp\{-\Delta \tau_{\nu}^{i}\}}{\Delta \tau_{\nu}^{i}}, \qquad (C.8)$$

$$\Delta \tau_{\nu}^{i} = \frac{\kappa_{\nu}^{\text{ext}}(s_{i}) + \kappa_{\nu}^{\text{ext}}(s_{i+1})}{2} (s_{i+1} - s_{i}),$$

and $S_{\nu} = \eta_{\nu}/\kappa_{\nu}^{\text{ext}}$. The values of $S_{\nu}$ and $\kappa_{\nu}^{\text{ext}}$ at the intersections $s_i$ and $s_{i+1}$ are estimated by linear interpolation from the grid-adjacent values (red dots in Fig. C.1). The optical depth, $\tau_{\nu}^i$, can be computed recursively,

$$\tau_{\nu}^{i+1} = \tau_{\nu}^{i} + \int_{s_i}^{s_{i+1}} \kappa_{\nu}^{\text{ext}} ds' = \tau_{\nu}^{i} + \Delta \tau_{\nu}^{i}, \qquad (C.9)$$

with $\tau_{\nu}^0 = 0$ and $\Delta \tau_{\nu}^i$ defined in Eq. (C.8).