Issue 
A&A
Volume 690, October 2024



Article Number  A192  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202244322  
Published online  09 October 2024 
Discontinuous Galerkin finite element method for the continuum radiative transfer problem inside axissymmetric circumstellar envelopes
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
^{★} Corresponding author; email: jeremy.perdigon@oca.eu
Received:
22
June
2022
Accepted:
16
August
2024
Context. The study of the continuum radiative transfer problem inside circumstellar envelopes is both a theoretical and numerical challenge, especially in the frequencydependent and multidimensional 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 frequencydependent twodimensional radiative transfer problem, and coupled with the radiative equilibrium equation. We next used this method in the context of axissymmetric 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 raytracing 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 axissymmetric circumstellar envelopes. Consecutively the emerging SEDs and images are also well reproduced. The DGFEM provides an alternative method (other than MonteCarlo 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
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 frequencydependent and multidimensional 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 fluxlimited 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 nonlinear 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, MonteCarlo 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.
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 highorder 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 nonsmooth behaviour (Nair et al. 2011), which often happens in the case of radiative transfer.
The DGFEM was successfully applied to the onedimensional (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 twodimensional (2D) radiative transfer have been developed, in cylindrical coordinates, in the context of a nonlocal thermodynamic equilibrium (nonLTE) 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 frequencydependent continuum radiative transfer problem, with isotropic scattering and coupled with the radiative equilibrium equation. We then used this method in the context of axissymmetric 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 raytracing techniques. The paper is organised as follows: in Sect. 2 we describe the radiative transfer problem for axissymmetric 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 axissymmetric 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 axissymmetric 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 Θ = π/2. Given the symmetries, the domain of definition of radiation is D ⊂ ℝ^{4}, with D ϶ x = [r, Θ, μ = cos θ, φ] ∈ [R_{in}, R_{out}] × (0, π/2] × [−1, 1] × [0, π]^{1}. The radiation field is described by the timeindependent radiation transfer equation (for a thorough derivation see e.g. Pomraning 1973, II5), $$\begin{array}{r}\mathbf{\Omega}.\mathrm{\nabla}{I}_{v}+{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}{I}_{v}={\eta}_{v},\\ \text{with}\mathbf{\Omega}.\mathrm{\nabla}{I}_{v}=\mu {\mathrm{\partial}}_{r}{I}_{v}+\frac{\sqrt{1{\mu}^{2}}\mathrm{cos}\phi}{r}{\mathrm{\partial}}_{\mathrm{\Theta}}{I}_{v}+\frac{1{\mu}^{2}}{r}{\mathrm{\partial}}_{\mu}{I}_{v}\\ \frac{\mathrm{cot}\mathrm{\Theta}\sqrt{1{\mu}^{2}}\mathrm{sin}\phi}{r}{\mathrm{\partial}}_{\phi}{I}_{v}.\end{array}$$(1)
The subscript v denotes the frequency dependence, I_{v} is the specific intensity, and Ω.∇ is the transport operator, which corresponds to a spatial derivative in the direction Ω. We note that ${\kappa}_{v}^{\text{ext}}={\kappa}_{v}^{\mathrm{abs}}+{\kappa}_{v}^{\text{sca}}$ is the extinction coefficient, with ${\kappa}_{v}^{\mathrm{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 crosssections and n is the number density. Dust grains are usually described as homogeneous spheres and their optical crosssections are commonly computed with the help of Mie theory (Bohren & Huffman 1998). The emissivity η_{v} includes the thermal emission (KirchoffPlanck 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 zerothorder angular moment of the specific intensity I_{v}, $${J}_{v}=\frac{1}{4\pi}\underset{4\pi}{\int}{I}_{v}\text{}d\mathbf{\Omega}=\frac{1}{2\pi}{\displaystyle \underset{0}{\overset{\pi}{\int}}}d\phi {\displaystyle \underset{1}{\overset{1}{\int}}}{I}_{v}d\text{}\mu .$$(2)
Additionally, the radiation field and the dust temperature are coupled via the equation of radiative equilibrium, $$4\pi {\displaystyle \underset{0}{\overset{\mathrm{\infty}}{\int}}}{\kappa}_{v}^{\mathrm{a}\mathrm{b}\mathrm{s}}\text{}{B}_{v}\text{}(T)\text{}dv=4\pi {\displaystyle \underset{0}{\overset{\mathrm{\infty}}{\int}}}{\kappa}_{v}^{\mathrm{a}\mathrm{b}\mathrm{s}}\text{}{J}_{v}\text{}dv.$$(3)
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}^{\star}+{I}_{v}^{\text{env}}$, where ${I}_{v}^{\star}$ 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{array}{l}\mathbf{\Omega}\cdot \mathrm{\nabla}{I}_{v}^{\star}+{\kappa}_{v}^{\text{ext}}{I}_{v}^{\star}=0\\ \mathbf{\Omega}.\mathrm{\nabla}{I}_{v}^{\text{env}}+{\kappa}_{v}^{\text{evt}}{I}_{v}^{\text{env}}={\eta}_{v}={\kappa}_{v}^{\text{abs}}{B}_{v}+{\kappa}_{v}^{\text{sca}}({J}_{v}^{\star}+{J}_{v}^{\text{env}}),\end{array}$$(4)
with ${J}_{v}^{\star}+{J}_{v}^{\text{env}}={J}_{v}$. We note that in Steinacker et al. (2003), the thermal emission term was put on the righthand 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}^{\star}$ 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}^{\star}=\{\begin{array}{ll}{\mathrm{\Gamma}}_{v}^{\star}\text{}({R}_{\star},\text{}\mu )\mathrm{exp}({\displaystyle \underset{0}{\overset{s\left({R}_{\text{in}}\right)}{\int}}}{\kappa}_{v}^{\mathrm{ext}}\left({s}^{\mathrm{\prime}}\right)d{s}^{\mathrm{\prime}}),& \text{if}{\mu}_{\star}\le \mu \le 1,\\ 0,& \text{otherwise.}\end{array}$$(5)
The incident stellar radiation field is ${\mathrm{\Gamma}}_{\nu}^{\star}({R}_{\star},\mu )$. Again, for the purpose of the numerical tests, we assume the star to radiate as a black body at the temperature T_{*}. Furthermore, μ_{*} = (1 − (R_{*}/r)^{2})^{1/2} is the cosine of the angle subtended by the star at radius r, and s(R_{in}) is the distance between a given point r and R_{in}, in the direction Ω. The argument in the exponential is the opposite of the optical depth integrated along the ray. For dusty media, we usually have r ≫ R_{*}, thus the star can be treated as a point source. In the point source approximation, μ* ≈ 1 − (R_{*}/r)^{2} /2, and the stellar mean intensity can be expressed analytically as, $${J}_{v}^{\star}\approx \frac{1}{4}\text{}{\left(\frac{{R}_{\star}}{r}\right)}^{2}\text{}{B}_{v}\left({T}_{\star}\right)\text{}\mathrm{exp}({\displaystyle \underset{{R}_{\text{in}}}{\overset{r}{\int}}}{\kappa}_{v}^{\mathrm{ext}}({r}^{\mathrm{\prime}},\mathrm{\Theta})d{r}^{\mathrm{\prime}}).$$(6)
In general, the integral in Eq. (6) can be carried out numerically, providing the stellar source term ${J}_{v}^{\star}$ for Eq. (4).
To complete the description of the problem, we specify the boundary conditions for ${I}_{v}^{\text{env}}$. We do so by prescribing the incident intensity ${\mathrm{\Gamma}}_{v}^{\text{env}}({\mathit{r}}_{s},\mathbf{\Omega})$ upon the surface of the domain D located at r_{s}, $${I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({\mathit{r}}_{s},\text{}\mathbf{\Omega})={\mathrm{\Gamma}}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({\mathit{r}}_{s},\text{}\mathbf{\Omega}),\mathrm{\forall}\hat{s}.\mathbf{\Omega}0,$$(7)
where ŝ is the unit vector normal to the surface of D, pointing outside of the domain. At the inner radius, $r={R}_{\mathrm{in}},\hat{s}=$ $\hat{r}$, and the incident radiation, in a given direction $\mathbf{\Omega}=\mu \hat{r}+$ $\sqrt{1{\mu}^{2}}\text{}\mathrm{cos}\text{}\phi \text{}\hat{\mathrm{\Theta}}+\sqrt{1{\mu}^{2}}\text{}\mathrm{sin}\text{}\phi \text{}\hat{\mathrm{\Phi}}$, comes directly from the opposite point of the cavity, $${I}_{v}^{\text{env}}({R}_{\mathrm{i}\mathrm{n}},\mathrm{\Theta},\text{}\mu ,\text{}\phi )={I}_{v}^{\text{env}}({R}_{\mathrm{i}\mathrm{n}},{\mathrm{\Theta}}^{\mathrm{\prime}},\text{}{\mu}^{\mathrm{\prime}},\text{}{\phi}^{\mathrm{\prime}}),\text{}\mathrm{\forall}0\mu {\mu}_{\star},$$(8)
with (R_{in}, Θ′, μ′, φ′) 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}},\mathrm{\Theta},\text{}\mu ,\text{}\phi )=0,\mathrm{\forall}\text{}\mu 0.$$(9)
At the equator, $\mathrm{\Theta}=\pi /2,\hat{s}=\hat{\mathrm{\Theta}}$, and the planar symmetry requires the radiation field to be, as shown in Fig. 2, $${I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}(r,\frac{\pi}{2},\text{}\mu ,\text{}\phi )={I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}(r,\frac{\pi}{2},\text{}\mu ,\text{}\pi \phi ),\mathrm{\forall}\text{}\phi \frac{\pi}{2}.$$(10)
Fig. 1 Illustration showing the coordinate system. We note that $(\hat{r},\hat{\mathrm{\Theta}},\hat{\mathrm{\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 Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates. 
Fig. 2 Representation of the direction vector Ω at the equator. The dotted arrows represent the symmetric of Ω with respect to the equatorial plane $(\hat{r},\hat{\mathrm{\Phi}})$. The left and right 2D views allow the corresponding direction angles to be seen. 
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, $${\mathrm{\nabla}}_{x}.\mathit{F}+{\kappa}^{\mathrm{ext}}\text{}\stackrel{~}{I}=\stackrel{~}{\eta},$$(11)
with $\stackrel{~}{\eta}={r}^{2}\mathrm{sin}\mathrm{\Theta}\eta $. We introduced the variable $\stackrel{~}{I}={r}^{2}\mathrm{sin}\mathrm{\Theta}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 (Θ = 0) where Eq. (1) is not defined when using the spherical coordinate system. The transport operator, ∇_{x}., corresponds to the Cartesian divergence operator with respect to x = (r, Θ, μ, φ), applied to the flux vector F of $\stackrel{~}{I}$, $$\mathit{F}=\mathit{a}\text{}\stackrel{~}{I}=\left(\begin{array}{l}{a}_{r}\\ {a}_{\mathrm{\Theta}}\\ {a}_{\mu}\\ {a}_{\phi}\end{array}\right)\stackrel{~}{I}=\left(\begin{array}{c}\mu \\ \sqrt{1{\mu}^{2}}\mathrm{cos}\phi /r\\ (1{\mu}^{2})/r\\ \mathrm{cot}\mathrm{\Theta}\sqrt{1{\mu}^{2}}\mathrm{sin}\phi /r\end{array}\right)\stackrel{~}{I}.$$(12)
We decomposed the domain D into N = N_{r} × N_{Θ} × N_{μ} × N_{φ} nonoverlapping rectangular elements D^{i,j,k,l}, with N_{r}, N_{Θ}, N_{μ} and N_{φ} being the number of elements along the r, Θ, μ and φ 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_{Θ} − 1, N_{μ} − 1 and, N_{φ} − 1. Inside each element D^{i,j,k,l}, we used the nodal representation and we approximated the local solution by a fourdimensional (4D) polynomial expansion, $${\stackrel{~}{I}}_{h}^{i,j,k,l}(\mathit{x})=\{\begin{array}{ll}{\displaystyle \sum _{a=0}^{{n}_{a}1}}{\displaystyle \sum _{b=0}^{{n}_{b}1}}{\displaystyle \sum _{c=0}^{{n}_{c}1}}{\displaystyle \sum _{d=0}^{{n}_{d}1}}{\stackrel{~}{I}}_{a,b,c,d}^{i,j,k}{h}_{a,b,c,d}(\mathit{x}),& \text{if}\mathit{x}\in {D}^{i,j,k,l},\\ 0,& \text{otherwise,}\end{array}$$(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, Θ, μ, and φ coordinate, respectively. We note that h_{a,b,c,d} is the 4D Lagrange polynomial, defined as, $$\begin{array}{rl}& {h}_{a,b,c,d}(\mathit{x})={h}_{a}\text{}(r)\text{}{h}_{b}\text{}(\mathrm{\Theta})\text{}{h}_{c}\text{}(\mu )\text{}{h}_{d}\text{}(\phi )\\ & =\prod _{{\scriptstyle \begin{array}{c}\alpha =0\\ \alpha \ne a\end{array}}}^{{n}_{a}1}\frac{r{r}_{\alpha}}{{r}_{a}{r}_{\alpha}}\prod _{{\scriptstyle \begin{array}{c}\beta =0\\ \beta \ne b\end{array}}}^{{n}_{b}1}\frac{\mathrm{\Theta}{\mathrm{\Theta}}_{\beta}}{{\mathrm{\Theta}}_{b}{\mathrm{\Theta}}_{\beta}}\prod _{{\scriptstyle \begin{array}{c}\gamma =0\\ \gamma \ne c\end{array}}}^{{n}_{c}1}\frac{\mu {\mu}_{\gamma}}{{\mu}_{c}{\mu}_{\gamma}}\prod _{{\scriptstyle \begin{array}{c}\delta =0\\ \delta \ne d\end{array}}}^{{n}_{d}1}\frac{\phi {\phi}_{\delta}}{{\phi}_{d}{\phi}_{\delta}}.\end{array}$$(14)
By definition, the coefficients ${\stackrel{~}{I}}_{a,b,c,d}^{i,j,k,l}={\stackrel{~}{I}}_{h}^{i,j,k,l}\left({\mathit{x}}_{a,b,c,d}\right)$ correspond to the value of ${\stackrel{~}{I}}_{h}^{i,j,k,l}$ at the nodes of coordinates x_{a,b,c,d} = (r_{a}, Θ_{b}, μ_{c}, φ_{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 piecewise continuous solutions inside each element, $$\stackrel{~}{I}(\mathit{x})\approx {\stackrel{~}{I}}_{h}(\mathit{x})={\displaystyle \sum _{i=0}^{{N}_{r}1}}{\displaystyle \sum _{j=0}^{{N}_{\mathrm{\Theta}}1}}{\displaystyle \sum _{k=0}^{{N}_{\mu}1}}{\displaystyle \sum _{l=0}^{{N}_{\phi}1}}{\stackrel{~}{I}}_{h}^{i,j,k,l}(\mathit{x}).$$(15)
Now we shall introduce ℛ_{h}, the residual of Eq. (11), $${\mathcal{R}}_{h}(\mathit{x})={\mathrm{\nabla}}_{x}.{\mathit{F}}_{h}+{\kappa}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}{\stackrel{~}{I}}_{h}\stackrel{~}{\eta},$$(16)
with ${\mathit{F}}_{h}=\mathit{a}{\stackrel{~}{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, $$\underset{{D}^{i,j,k,l}}{\int}}{\mathcal{R}}_{h}\text{}(\mathit{x})\text{}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}(\mathit{x})\text{}{d}^{4}\mathit{x}=0,\text{}\mathrm{\forall}\text{}{D}^{i,j,k,l}\text{, and}\mathrm{\forall}\text{}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{.$$(17)
The divergence term that appears in Eq. (17) can be recast with the help of the divergence theorem (GreenOstrogradsky), yielding the following socalled weak formulation of the radiative transfer equation Eq. (11): $$\begin{array}{rl}& {\displaystyle \underset{\mathrm{\partial}{D}^{i,j,k,l}}{\oint}}\text{}\hat{s}.{\mathit{F}}^{\ast}\text{}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{\mathrm{d}}^{3}\mathit{x}{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{\mathit{F}}_{h}.\mathrm{\nabla}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{\mathrm{d}}^{4}\mathit{x}\\ & +{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}({\kappa}^{\mathrm{e}\mathrm{x}\mathrm{t}}{\stackrel{~}{I}}_{h}\stackrel{~}{\eta}){h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{\mathrm{d}}^{4}\mathit{x}=0,\text{}\mathrm{\forall}\text{}{D}^{i,j,k,l}\text{, and}\mathrm{\forall}\text{}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{.}\end{array}$$(18)
The first term, in the lefthand side of Eq. (18), is a surface integral, along the boundaries of D^{i,j,k,l}. Furthermore, ŝ is the outward normal vector to the surface element, and 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 × 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 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 LaxFriedrichs numerical flux, defined as (e.g. on the radial right element edge where $\hat{s}=\hat{r}$), $$\hat{r}.{\mathit{F}}^{\ast}=\frac{1}{2}({F}_{r}\left({\stackrel{~}{I}}^{}\right)+{F}_{r}\left({\stackrel{~}{I}}^{+}\right)\left{a}_{r}\right({\stackrel{~}{I}}^{+}{\stackrel{~}{I}}^{})),$$(19)
where ${\stackrel{~}{I}}^{}$ and ${\stackrel{~}{I}}^{+}$ denote the left and right values of $\stackrel{~}{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 × n equations with n = n_{a} × n_{b} × n_{c} × n_{d}, relating the ${\stackrel{~}{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 problemdependent (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}{\stackrel{~}{\mathit{I}}}_{h}=\mathcal{B},$$(20)
with ${\stackrel{~}{\mathit{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′ × N′ coupling the elements of ${\stackrel{~}{\mathit{I}}}_{h}$, and ℬ a vector of size N′ containing the emissivity term $\stackrel{~}{\eta}$.
4 Solution strategy and numerical considerations
The solution of the problem presented in Sec. 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 Sec. 5, the star was treated in the point source approximation, and hence we directly used Eq. (6) to compute ${J}_{v}^{\star}$. The radiation field of the envelope ${I}_{v}^{\text{env}}$, described by the second Eq. (4), was solved with the DGFEM presented in Sec. 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, ${I}_{a,b,c,d}^{i,j,k,l}$, inside a given element, D^{i,j,kl}, only depends on the other points inside the same element and the neighbouring ones, inside D^{i±}^{1,j±1,k−1,l+1}. This property allows us to rewrite Eq. (20) as, $$\begin{array}{rl}& {\mathcal{A}}^{i,j,k,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l}+{\mathcal{A}}^{i+1,j,k,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i+1,j,k,l}+{\mathcal{A}}^{i1,j,k,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i1,j,k,l}\\ & +{\mathcal{A}}^{i,j+1,k,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i,j+1,k,l}+{\mathcal{A}}^{i,j1,k,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i,j1,k,l}+{\mathcal{A}}^{i,j,k1,l}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k1,l}\\ & +{\mathcal{A}}^{i,j,k,l+1}\text{}{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l+1}={\mathcal{B}}^{i,j,k,l},\text{}\mathrm{\forall}\text{}i,\text{}j,\text{}k,\text{and}l.\end{array}$$(21)
We note that $\mathcal{A}$^{i,j,k,l} are the diagonal blocks of $\mathcal{A}$, with a size of n × n while $\mathcal{A}$^{i±}^{1,j±1,k−1,l+1} (with a size of n × n) are the only nonzero, nondiagonal blocks. Furthermore, the elements D^{i}^{,}^{j}^{,}^{k}^{+1}^{,}^{1} and D^{i}^{,}^{j}^{,}^{k}^{,}^{l}^{−1} do not contribute because of the expression of the LaxFriedrichs numerical flux, Eq. (19), with a_{μ} ≥ 0 and a_{φ} ≤ 0, ∀x ∈ D. We note that ${\stackrel{~}{I}}_{h}^{i,j,k,l}$ and ℬ^{i,j,k,l} are the subvector of ${\stackrel{~}{\mathit{I}}}_{h}$ and ℬ, 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 α = i N_{Θ} N_{μ} + j N_{μ} N_{φ}, + k N_{φ} + 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 Λiteration in the literature, becomes very slow and does not converge for large optical depths (Mihalas & WeibelMihalas 1999, VI83). 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 Λiteration and yields satisfying results up to moderately thick envelopes. We however note that, as for the usual Λ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}^{\star}$ with Eq. (6) and set ${\left[{I}_{v}^{\text{env}}\right]}^{n=0}=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}^{\star}$. Then (ii), for each frequency, we computed ${\left[{I}_{V}^{\text{env}}\right]}^{n+1}$ with the help of Eq. (21), which we rewrote, performing a block GaussSeidel sweep (see e.g Karniadakis & Kirby II 2003, VII2), $$\begin{array}{rl}& {\mathcal{A}}^{i,j,k,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l}\right]}^{n+1}={\left[{\mathcal{B}}^{i,j,k,l}\right]}^{n}{\mathcal{A}}^{i+1,j,k,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i+1,j,k,l}\right]}^{n}\\ & \phantom{\rule{2em}{0ex}}{\mathcal{A}}^{i1,j,k,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i1,j,k,l}\right]}^{n+1}\text{}{\mathcal{A}}^{i,j+1,k,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i,j+1,k,l}\right]}^{n}\\ & \phantom{\rule{1em}{0ex}}{\mathcal{A}}^{i,j1,k,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i,j1,k,l}\right]}^{n+1}{\mathcal{A}}^{i,j,k1,l}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k1,l}\right]}^{n+1}\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{\mathcal{A}}^{i,j,k,l+1}{\left[{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l+1}\right]}^{n}.\end{array}$$(22)
We give in appendix B the expressions of $\mathcal{A}$^{i,j,k,l} and of the righthand 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, IX1). (iii) We computed ${J}_{v}^{\text{env}}$ and consequently updated the temperature via Eq. (3). The new temperature allowed us to update the righthand 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 Λ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).
Summary of the main results.
Fig. 3 Example of the sparse structure of $\mathcal{A}$ with N_{r} = N_{Θ} = N_{μ} = N_{φ} = 2. The squares represent the blocks of size n × n containing nonzero values. The matrix has bands corresponding to the coupling blocks ${\mathcal{A}}^{\alpha ,\alpha},{\mathcal{A}}^{\alpha ,\alpha +1},{\mathcal{A}}^{\alpha ,\alpha {N}_{\phi}},{\mathcal{A}}^{\alpha ,\alpha \pm {N}_{\mu}{N}_{\phi}}$, and ${\mathcal{A}}^{\alpha ,\alpha \pm {N}_{\mathrm{\Theta}}{N}_{\mu}{N}_{\phi}}$. 
5 Numerical tests
The lack of analytic solutions for the radiative transfer problem, especially in the multidimensional and frequencydependent 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 frequencydependent radiative transfer problem inside a spherically symmetric envelope, from Ivezic et al. (1997). The second test case, in Sec. 5.2, is about the frequencydependent radiative transfer inside an axissymmetric envelope (disc), from Pascucci et al. (2004). For the latter test, we note that five different codes were used (including Monte Carlo and gridbased 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 Sec. 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_{in} to the outer radius R_{out} = 1000 R_{in}. The inner radius was set so that its temperature would always be T_{in} = T(R_{in}) = 800 K. The number density profile, n(r), is assumed to be a power law, of the form n(r) = n_{0} (R_{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 though the envelope, at v_{0} = c/λ_{0}, with λ_{0} = 1 μm, $${\tau}_{{v}_{0}}={\displaystyle \underset{{R}_{\text{in}}}{\overset{{R}_{\text{out}}}{\int}}}{\kappa}_{{v}_{0}}^{\text{ext}}dr={C}_{{v}_{0}}^{\text{ext}}{\displaystyle \underset{{R}_{\text{in}}}{\overset{{R}_{\text{out}}}{\int}}}n\text{}dr={C}_{{v}_{0}}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}{n}_{0}\text{}{R}_{\text{in}}(1\frac{{R}_{\text{in}}}{{R}_{\text{out}}}).$$(23)
We note that ${C}_{{\nu}_{0}}^{\text{ext}}={C}_{{\nu}_{0}}^{\text{abs}}+{C}_{{\nu}_{0}}^{\text{sca}}$ is the extinction crosssection coefficient at v_{0}. In this test we consider ${\tau}_{{v}_{0}}=1$ and 10^{2} to correspond to a moderately thin and thick envelopes, respectively. The absorption and scattering crosssections, ${C}_{v}^{\text{abs}}$ and ${C}_{v}^{\text{sca}}$, feature a bilinear behaviour in loglog scaling, with a constant profile for v ≥ v_{0} and a powerlaw dependence ∝ (v/v_{0})^{α}, for v ≤ v_{0}, with α = 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^{2}. Although the envelope is spherically symmetric, we used a grid of N = N_{r} × N_{Θ} × N_{μ} × N_{φ} = 16^{4} elements in our DGFEM code, with each elements containing n_{a} ×n_{b} × n_{c} × n_{d} = 3 × 2 × 3 × 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 logaritmically spaced points, ranging from λ = 10^{−2} μm to λ = 3.6 × 10^{4} μm.
The temperature profiles, T, and the SEDs, λF_{λ}/F with $F={\int}_{0}^{\mathrm{\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 Θ. Additionally, the normalised SEDs were computed with the help of a raytracing 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 loglog 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.
Fig. 4 Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with ${\tau}_{{v}_{0}}=1$ (blue curve) and ${\tau}_{{v}_{0}}={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. 
5.2 2D axissymmetric envelope
A point source surrounded by an axissymmetric disc of dust, at radiative equilibrium, radiates as a black body at the temperature T_{*} = 5800 K. This disc extends from the inner radius R_{in} = 1 AU to the outer radius R_{out} = 1000 AU. The density profile n(r, Θ) is assumed to be $$n(r,\mathrm{\Theta})={n}_{0}\frac{{r}_{d}}{r\mathrm{sin}\mathrm{\Theta}}\mathrm{exp}(\frac{\pi}{4}{\left(\frac{r\mathrm{cos}\mathrm{\Theta}}{h}\right)}^{2}),$$(24)
with $h={z}_{d}{(r\mathrm{sin}\mathrm{\Theta}/{r}_{d})}^{\frac{9}{8}},{r}_{d}={R}_{\text{out}}/2$, and z_{d} = R_{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}_{{v}_{0}}$, the radial optical depth, integrated through the disc midplane (Θ = π/2), at frequency v_{0} = c/λ_{0} with λ_{0} = 0.55 μm, $${n}_{0}=\frac{{\tau}_{{v}_{0}}}{{C}_{{v}_{0}}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}{r}_{d}\text{}\mathrm{ln}\left(\frac{{R}_{\text{out}}}{{R}_{\text{in}}}\right)}.$$(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 g cm^{−3}. A table of precomputed 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}_{{\nu}_{0}}={10}^{1}$ and 10^{2}, respectively. The benchmarks were produced with RADMC3D^{3} (Dullemond et al. 2012), a MonteCarlo 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 RADMC3D, 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 highorder numerical scheme. The radiation field can exhibit discontinuities, because of boundary conditions 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 midplane is very well reproduced, highlighting that the method is able to correctly represent the shadowed regions of the disc (the cold outer midplane 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 raytracing of the emissivity η_{v}. They consist of a stellar and an envelope contribution (see Eq. (C.1)). The stellar component that peaks at around λ ≈ 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 gridbased and a MonteCarlo 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 λ = 2.3, 4.5, and 12.1 μm 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 RADMC3D 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 λ = 2.3, 4.5, and 12.1 μm. In general, the disc emitting inner regions (peaks in Fig. 7) are reproduced very well $(\u03f5\left({I}_{v}^{\text{env}}\right)\le 3\mathrm{\%})$. The biggest discrepancies occur in the wings of the peaks, where the gradient of intensity is the steepest (in logarithmic scaling).
Fig. 5 Temperature profiles for the axissymmetric envelope with ${\tau}_{{v}_{0}}={10}^{1}$ (blue curve) and ${\tau}_{{v}_{0}}={10}^{2}$ (orange curve), in the in the disc midplane (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 RADMC3D. The lower panels show the relative differences between the two codes. 
Fig. 6 SED profiles for the axissymmetric envelope with ${\tau}_{{v}_{0}}={10}^{1}$ (blue curve) and ${\tau}_{{v}_{0}}={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 RADMC3D. The lower panels show the relative differences between the two codes. 
Fig. 7 Images at λ = 2.3, 4.5, and 12.1 μm (left, middle, and right panels, respectively) of the 10 AU inner regions of the axissymmetric envelope (Sect. 5.2) and computed from the DGFEM solution with the raytracing 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}_{v}^{\text{env}}$, inside each pixel. The solid white lines show the isocontours of ${I}_{v}^{\text{env}}$. 
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 raytracing procedure to accurately compute SEDs. However, in practice, we already know the intensity I_{v}(r = R_{out}, Θ, μ, φ) 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 postprocessing treatment required. It is, however, not advised to proceed in this way, at least in multidimensional 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_{out} ΔΘ 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 postprocessing raytracing step is usually needed in order to capture all the features of the disc accurately. The postprocessing raytracing 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 blackbody emission B_{v}(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_{DG} and the RADMC3D code T_{ref}. It is defined as, $${L}^{2}={\left(\frac{{\displaystyle \underset{\mathcal{V}}{\int}}{\left(\frac{{T}_{\mathrm{D}\mathrm{G}}{T}_{\mathrm{r}\mathrm{e}\mathrm{f}}}{{T}_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)}^{2}{r}^{2}\mathrm{sin}\mathrm{\Theta}\text{}dr\text{}d\mathrm{\Theta}}{{\displaystyle \underset{\mathcal{V}}{\int}}{r}^{2}\mathrm{sin}\mathrm{\Theta}\text{}dr\text{}d\mathrm{\Theta}}\right)}^{\frac{1}{2}},$$(26)
where the integral runs across the volume $\mathcal{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} × N_{Θ} × N_{μ} × N_{φ} 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_{Θ} = N_{μ} = N_{φ}).
We can see that the L^{2} norm saturates when N ≈ 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), shortcharacteristics, (Dullemond & Turolla 2000) or MonteCarlo (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 postprocessing raytracing or directly by using I_{v}(r = R_{out}, Θ, μ, φ). We compare these SEDs with the benchmark curves we already presented in Sects. 5.1 and 5.2 (τ = 100). While for the spherically symmetric case, they agree quite well, I_{v}(r = R_{out}, Θ, μ, φ) 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 Θ and φ variables. We note that the increase in the grid size should make the SED from I_{v}(r = R_{out}, Θ, μ, φ) converge to the raytracing 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 posttreatment raytracing is thus needed if one wants to compute the emerging fluxes from the disc with precision.
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 RADMC3D 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. 
Fig. 9 L^{2} norm of the relative differences in temperature between this study and RADMC3D versus the number of elements in the computational grid (N = N_{r} × N_{Θ} × N_{μ} × N_{φ}. The L^{2} norm was computed for the benchmark problem presented in Sec. 5.2. 
6 Conclusions
For this study, we applied the DGFEM to the frequencydependent 2D radiative transfer equation, coupled with the radiative equilibrium equation, inside axissymmetric 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 raytracing 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 highorder 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 MonteCarlo and raytracing techniques, the DGFEM provides 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.
Fig. 10 Emerging SEDs of RADMC3D (solid black lines) and of the DGFEM code, computed from the postprocessing raytracing 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 (τ = 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. 
Acknowledgements
This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côted’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 RADMC3D: https://www.ita.uniheidelberg.de/~dullemond/software/radmc3d/ Part of the computations were carried out with the help of OPALMeso 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.
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}, Θ) on the inner cavity of radius R_{in} and in a direction $\mathbf{\Omega}(\mu ,\phi )(\mathbf{\Omega}.\hat{r}>0)$, the coordinates of the opposite point B (R_{in}, Θ′, μ′ = cos θ′, φ′). 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 axissymmetric around the polar axis and the radiation field does not depend on the azimutal angle Φ (see Fig. 1). Consequently, we can arbitrarily choose the Φ coordinate of the point A and set it to be Φ = π/2. This convention simplifies the computations.
Fig. A.1 Geometry of a ray in the inner cavity. Point B is the opposite point of A, in the direction Ω. 
First, we see that θ′ = π − θ because the OAB triangle is isosceles, we then have, $${\mu}^{\mathrm{\prime}}=\mu .$$(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 $${\mathit{r}}_{B}={\mathit{r}}_{A}{s}_{AB}\text{}\mathbf{\Omega},$$(A.2)
with s_{AB} the distance between the point A and B. This length is equal to s_{AB} = 2 μ R_{in}, because it corresponds to the base of the isosceles triangle OAB, with OA = OB = R_{in}. The direction vector, Ω, is, $$\begin{array}{rl}& \mathbf{\Omega}={\left[\begin{array}{cc}\mu & \\ \sqrt{1{\mu}^{2}}& \mathrm{cos}\phi \\ \sqrt{1{\mu}^{2}}& \mathrm{sin}\phi \end{array}\right]}_{(\hat{r},\hat{\mathrm{\Theta}},\hat{\mathrm{\Phi}})}\\ & ={\left[\begin{array}{c}\sqrt{1{\mu}^{2}}\mathrm{sin}\phi \\ \mathrm{sin}\mathrm{\Theta}\mu +\mathrm{cos}\mathrm{\Theta}\sqrt{1{\mu}^{2}}\mathrm{cos}\phi \\ \mathrm{cos}\mathrm{\Theta}\mu \mathrm{sin}\mathrm{\Theta}\sqrt{1{\mu}^{2}}\mathrm{cos}\phi \end{array}\right]}_{(\hat{x},\hat{y},\hat{z})}.\end{array}$$(A.3)
Eq. (A.2) is then $$\begin{array}{rl}{\mathit{r}}_{B}& ={R}_{\mathrm{i}\mathrm{n}}{\left[\begin{array}{c}\mathrm{sin}{\mathrm{\Theta}}^{\mathrm{\prime}}\mathrm{cos}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ \mathrm{sin}{\mathrm{\Theta}}^{\mathrm{\prime}}\mathrm{sin}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ \mathrm{cos}{\mathrm{\Theta}}^{\mathrm{\prime}}\end{array}\right]}_{(\hat{x},\hat{y},\hat{z})}\\ & ={R}_{\mathrm{i}\mathrm{n}}{\left[\begin{array}{c}2\mu \sqrt{1{\mu}^{2}}\mathrm{sin}\phi \\ \mathrm{sin}\mathrm{\Theta}(12{\mu}^{2})2\mathrm{cos}\mathrm{\Theta}\mu \sqrt{1{\mu}^{2}}\text{\hspace{0.17em}cos\hspace{0.17em}}\phi \\ \text{\hspace{0.17em}cos\hspace{0.17em}}\mathrm{\Theta}(12{\mu}^{2})+2\text{\hspace{0.17em}sin}\mathrm{\Theta}\mu \sqrt{1{\mu}^{2}}\text{\hspace{0.17em}cos\hspace{0.17em}}\phi \end{array}\right]}_{(\hat{x},\hat{y},\hat{z})}.\end{array}$$(A.4)
From the zcomponent of Eq. (A.4) we get, $$\mathrm{cos}{\mathrm{\Theta}}^{\mathrm{\prime}}=\mathrm{cos}\mathrm{\Theta}(12{\mu}^{2})+2\mathrm{sin}\mathrm{\Theta}\mu \sqrt{1{\mu}^{2}}\text{\hspace{0.17em}cos\hspace{0.17em}}\phi ,$$(A.5)
and consequently Θ′ = arccos  cos Θ′. The absolute value occurs if we consider the planar symmetry with respect to the equatorial plane.
The angle φ′ at the point B verifies, $$\mathrm{tan}{\phi}^{\mathrm{\prime}}=\frac{\mathbf{\Omega}\cdot {\hat{\mathrm{\Phi}}}^{\mathrm{\prime}}}{\mathbf{\Omega}\cdot {\hat{\mathrm{\Theta}}}^{\mathrm{\prime}}},$$(A.6)
with, $${\hat{\mathrm{\Phi}}}^{\mathrm{\prime}}={\left[\begin{array}{c}\mathrm{sin}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ \mathrm{cos}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ 0\end{array}\right]}_{(\hat{x},\hat{y},\hat{z})},{\hat{\mathrm{\Theta}}}^{\mathrm{\prime}}={\left[\begin{array}{c}\mathrm{cos}{\mathrm{\Theta}}^{\mathrm{\prime}}\mathrm{cos}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ \mathrm{cos}{\mathrm{\Theta}}^{\mathrm{\prime}}\mathrm{sin}{\mathrm{\Phi}}^{\mathrm{\prime}}\\ \mathrm{sin}{\mathrm{\Theta}}^{\mathrm{\prime}}\end{array}\right]}_{(\hat{x},\hat{y},\hat{z})},$$(A.7)
and Ω given by Eq. (A.3). Then, using Eq. (A.4) in Eq. (A.6) yields, after tedious calculations, $$\mathrm{tan}{\phi}^{\mathrm{\prime}}=\frac{\mathrm{sin}\phi \mathrm{sin}\mathrm{\Theta}}{(12\text{}{\mu}^{2})\mathrm{cos}\phi \mathrm{sin}\mathrm{\Theta}2\mu \text{}\sqrt{1{\mu}^{2}}\mathrm{cos}\mathrm{\Theta}}.$$(A.8)
Consequently φ′ = arctan u/υ with u and υ 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 φ′ angle to take its values in [0, π]. We note that the previous formula is not valid for the couples of coordinates (μ, φ) = (sin (Θ/2), 0) and (μ, φ) = (cos (Θ/2), π), because it corresponds to have the opposite point B on the poles, where ${\hat{\mathrm{\Theta}}}^{\prime}$ and ${\hat{\mathrm{\Phi}}}^{\prime}$ are not defined. These particular cases are in practice not a problem since the radiation field is independent of φ′ 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, Θ, μ, φ). For the r, μ, φ coordinates, we make use of the GaussLobatto 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 ${\mathit{F}}_{h}^{\ast}$ at the element edges. The GaussLobatto quadrature can exactly integrate polynomials up to degree 2 n − 3, with n the number of nodes. For the coordinate Θ, we cannot use the GaussLobatto quadrature, at least in the cells that are touching Θ = 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 Θ, we use a GaussLegendre quadrature for this coordinate. The GaussLegendre 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.
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, Θ) and (μ, φ) 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}^{\ast}$ at the interface, along the Θ coordinate. 
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{array}{r}{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}({\kappa}^{\mathrm{e}\mathrm{x}\mathrm{t}}{\stackrel{~}{I}}_{h}\stackrel{~}{\eta}){h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{d}^{4}\mathit{x}\\ =\frac{\mathrm{\Delta}{x}^{i,j,k,l}}{16}{\displaystyle \underset{1}{\overset{1}{\int}}}({\kappa}^{\mathrm{e}\mathrm{x}\mathrm{t}}{\stackrel{~}{I}}_{h}^{i,j,k,l}{\stackrel{~}{\eta}}^{i,j}){h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{d}^{4}\stackrel{~}{\mathit{x}}\\ =\frac{\mathrm{\Delta}{x}^{i,j,k,l}}{16}{\displaystyle \sum _{a,b,c,d}}{W}_{a,b,c,d}({\kappa}_{a,b}^{{\mathrm{e}\mathrm{x}\mathrm{t}}^{i,j}}{\stackrel{~}{I}}_{a,b,c,d}^{i,j,k,l}{\stackrel{~}{\eta}}_{a,b}^{i,j}){h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\left({\stackrel{~}{\mathit{x}}}_{a,b,c,d}\right).\end{array}$$(B.1)
We note that Δx^{i,j,k,l} = Δr^{i}, Δ Θ^{j} Δμ^{k} Δφ^{l} is the 4D volume of the element D^{i,j,k,l} For integration, we rather use the local coordinates $(\stackrel{~}{r},\stackrel{~}{\mathrm{\Theta}},\stackrel{~}{\mu},\stackrel{~}{\phi})$, defined as (e.g for the r coordinate), $$\stackrel{~}{r}=\frac{2}{\mathrm{\Delta}{r}^{i}}(r{r}^{i+1/2})\text{,}$$(B.2)
with Δ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}\left({\stackrel{~}{r}}_{a}\right){W}_{\mathrm{\Theta}}\left({\stackrel{~}{\mathrm{\Theta}}}_{b}\right){W}_{\mu}\left({\stackrel{~}{\mu}}_{c}\right){W}_{\phi}\left({\stackrel{~}{\phi}}_{d}\right)$ are the weights associated with the different quadrature in each direction. Finally, ${h}_{{a}^{\prime},{b}^{\prime},{c}^{\prime},{d}^{\prime}}\left({\stackrel{~}{\mathit{x}}}_{a,b,c,d}\right)$ is the 4D Lagrange polynomials, defined in Eq. (14), evaluated at the node ${\stackrel{~}{\mathit{x}}}_{a,b,c,d}=({\stackrel{~}{r}}_{a},{\stackrel{~}{\mathrm{\Theta}}}_{b},{\stackrel{~}{\mu}}_{c},{\stackrel{~}{\phi}}_{d})$. By definition of the Lagrange polynomials, we have, $${h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}({\stackrel{~}{r}}_{a},{\stackrel{~}{\mathrm{\Theta}}}_{b},\text{}{\stackrel{~}{\mu}}_{c},\text{}{\stackrel{~}{\phi}}_{d})={\delta}_{{a}^{\mathrm{\prime}},a}\text{}{\delta}_{{b}^{\mathrm{\prime}},b}\text{}{\delta}_{{c}^{\mathrm{\prime}},c}\text{}{\delta}_{{d}^{\mathrm{\prime}},d},$$(B.3)
with δ_{a}_{′,a} the usual delta Kronecker. The other volume term in Eq. (18) is expressed as, $$\begin{array}{rl}& \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{\mathit{F}}_{h}.{\mathrm{\nabla}}_{x}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{d}^{4}\mathit{x}\\ & ={\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{a}_{r}\text{}{\stackrel{~}{\mathit{I}}}_{h}\text{}{\mathrm{\partial}}_{r}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{d}^{4}\mathit{x}+{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{a}_{\mathrm{\Theta}}\text{}{\stackrel{~}{\mathit{I}}}_{h}\text{}{\mathrm{\partial}}_{\mathrm{\Theta}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{d}^{4}\mathit{x}\\ & +{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{a}_{\mu}\text{}{\stackrel{~}{\mathit{I}}}_{h}\text{}{\mathrm{\partial}}_{\mu}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{d}^{4}\mathit{x}+{\displaystyle \underset{{D}^{i,j,k,l}}{\int}}{a}_{\phi}{\stackrel{~}{\mathit{I}}}_{h}\text{}{\mathrm{\partial}}_{\phi}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{d}^{4}\mathit{x}\\ & =\frac{\mathrm{\Delta}{x}^{i,j,k,l}}{8}\sum _{a,b,c,d}{W}_{a,b,c,d}\text{}{\stackrel{~}{I}}_{a,b,c,d}^{i,j,k,l}(\frac{{a}_{{\stackrel{~}{r}}_{a,b,c,d}}^{i,j,k,l}\text{}{\mathrm{\partial}}_{\stackrel{~}{r}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{{l}^{\mathrm{\prime}}}}{\stackrel{~}{x}}_{a,b,c,d}}{\mathrm{\Delta}{r}^{i}}\\ & +\frac{{a}_{{\stackrel{~}{\mathrm{\Theta}}}_{a,b,c,d}}^{i,j,k,l}\text{}{\mathrm{\partial}}_{\stackrel{~}{\mathrm{\Theta}}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{\stackrel{~}{x}}_{a,b,c,d}}{\mathrm{\Delta}{\mathrm{\Theta}}^{j}}+\frac{{a}_{{\stackrel{~}{\mu}}_{a,b,c,d}}^{i,j,k,l}\text{}{\mathrm{\partial}}_{\stackrel{~}{\mu}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{\stackrel{~}{x}}_{a,b,c,d}}{\mathrm{\Delta}{\mu}^{k}}\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}+\frac{{a}_{\stackrel{~}{\phi}a,b,c,d}^{i,j,k,l}\text{}{\mathrm{\partial}}_{\stackrel{~}{\phi}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{\stackrel{~}{x}}_{a,b,c,d}}{\mathrm{\Delta}{\phi}^{l}}).\end{array}$$(B.4)
In Eq. (B.4), we used the definition of the flux Eq. (12). The quantity ${\mathrm{\partial}}_{\stackrel{~}{r}}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{x}_{a,b,c,d}\text{}$ is the partial derivative of the Lagrange polynomial, with respect to the coordinate r, evaluated at the node ${\stackrel{~}{\mathit{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 × 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 $(\stackrel{~}{r}=\pm 1)$, respectively, $$\begin{array}{r}{\displaystyle \underset{\mathrm{\partial}{D}^{i,j,k,l}}{\oint}}{\left[{F}_{r}^{\ast}{h}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\right]}_{\stackrel{~}{r}=1}^{\stackrel{~}{r}=1}d\mathrm{\Theta}\text{}d\mu \text{}d\phi \\ =\frac{\mathrm{\Delta}{x}^{i,j,k,l}}{8\text{}\mathrm{\Delta}{r}^{i}}{W}_{{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{[{F}_{\stackrel{~}{r}}^{\ast}(\stackrel{~}{r},{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}){h}_{{a}^{\mathrm{\prime}}}(\stackrel{~}{r})]}_{\stackrel{~}{r}=1}^{\stackrel{~}{r}=1}.\end{array}$$(B.5)
At the cell boundaries, we use the upwind numerical flux, for example, at the right edge, $$\begin{array}{r}{F}_{\stackrel{~}{r}}^{\ast}(\stackrel{~}{r}=1,{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}})=\\ max\{{a}_{\stackrel{~}{r}}^{i,j,k,l}(1,{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}),0\}{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l}(1,{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}})\\ +min\{{a}_{\stackrel{~}{r}}^{i,j,k,l}(1,{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}),0\}{\stackrel{~}{\mathit{I}}}_{h}^{i+1,j,k,l}(1,{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}).\end{array}$$(B.6)
The use of the GaussLobatto quadrature further simplify the computation because the quadrature nodes include the endpoints $\stackrel{~}{r}=\pm 1$. The numerical flux can then be expressed as $$\begin{array}{rl}& \phantom{\rule{2em}{0ex}}{[{F}_{\stackrel{~}{r}}^{\ast}(\stackrel{~}{r},{\stackrel{~}{\mathrm{\Theta}}}_{{b}^{\mathrm{\prime}}},{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}){h}_{{a}^{\mathrm{\prime}}}(\stackrel{~}{r})]}_{\stackrel{~}{r}=1}^{\stackrel{~}{r}=1}=\\ & {\delta}_{{a}^{\mathrm{\prime}},{n}_{a}1}max\{{a}_{{\stackrel{~}{r}}_{{n}_{a}}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l},0\}{\stackrel{~}{I}}_{{n}_{a}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l}\\ & \phantom{\rule{1em}{0ex}}+{\delta}_{{a}^{\mathrm{\prime}},{n}_{a}1}min\{{a}_{{\stackrel{~}{r}}_{{n}_{a}}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l},0\}{\stackrel{~}{I}}_{0,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i+1,j,k,l}\\ & \phantom{\rule{2em}{0ex}}{\delta}_{{a}^{\mathrm{\prime}},0}max\{{a}_{{\stackrel{~}{r}}_{0},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l},0\}{\hat{I}}_{{n}_{a}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i1,j,k,l}\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}{\delta}_{{a}^{\mathrm{\prime}},0}min\{{a}_{{\stackrel{~}{r}}_{0},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l},0\}{\hat{I}}_{0,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l}.\end{array}$$(B.7)
The same form holds for the surface integrals normal to the μ and φ coordinates. We however note that we have a_{μ} ≥ 0, a_{φ} ≤ 0 ∀ x ∈ D, so there are no terms proportional to ${\stackrel{~}{I}}^{i,j,k+1,l}$ and ${\stackrel{~}{I}}^{i,j,k,l1}$. For the flux computation normal to the Θ 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{array}{rl}& \phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{F}_{\stackrel{~}{\mathrm{\Theta}}}^{\ast}({\stackrel{~}{r}}_{{a}^{\mathrm{\prime}}},\stackrel{~}{\mathrm{\Theta}}=1,{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}})=\\ & \phantom{\rule{2em}{0ex}}max\{{a}_{\stackrel{~}{\mathrm{\Theta}}}^{i,j,k,l}({\stackrel{~}{r}}_{{a}^{\mathrm{\prime}}},1,{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}),0\}\sum _{b}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},b,{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j,k,l}\text{}{h}_{b}(1)\\ & +min\{{a}_{\stackrel{~}{\mathrm{\Theta}}}^{i,j,k,l}({\stackrel{~}{r}}_{{a}^{\mathrm{\prime}}},1,{\stackrel{~}{\mu}}_{{c}^{\mathrm{\prime}}},{\stackrel{~}{\phi}}_{{d}^{\mathrm{\prime}}}),0\}\sum _{b}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},b,{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j+1,k,l}\text{}{h}_{b}(1).\end{array}$$(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,}{\stackrel{~}{\mathit{I}}}_{h}^{i,j,k,l}={\mathit{b}}^{i,j,k,l}.$$(B.9)
In Eq. (B.9), we replaced the RHS by b^{i,j,k,l} because we do not need to formally write the nondiagonal matrices $\mathcal{A}$^{i}^{±1,j±1,k−1,l+1}. To assemble $\mathcal{A}$^{i,j,k,l}, we make use of the global index α = a n_{b}n_{c}n_{d} + b n_{c}n_{d} + c n_{d} + d. The elements of $\mathcal{A}$^{i,j,k,l} and b^{i,j,k,l} are then, $$\begin{array}{rl}& \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{\mathcal{A}}_{{\alpha}^{\mathrm{\prime}}\alpha}^{i,j,k,l}={W}_{a,b,c,d}\text{}{\kappa}_{a,b}^{{\mathrm{e}\mathrm{x}\mathrm{t}}^{i,j}}{\delta}_{{a}^{\mathrm{\prime}},a}\text{}{\delta}_{{b}^{\mathrm{\prime}},b}\text{}{\delta}_{{c}^{\mathrm{\prime}},c}\text{}{\delta}_{{d}^{\mathrm{\prime}},d}\\ & \text{}\phantom{\rule{1em}{0ex}}+\frac{2}{\mathrm{\Delta}{r}^{i}}{\delta}_{{b}^{\mathrm{\prime}},b}{\delta}_{{c}^{\mathrm{\prime}},c}{\delta}_{{d}^{\mathrm{\prime}},d}{W}_{b,c,d}(max\{{a}_{{\stackrel{~}{r}}_{a,b,c,d}}^{i,j,k,l},0\}{\delta}_{{a}^{\mathrm{\prime}},{n}_{a}1}{\delta}_{a,{n}_{a}1}\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}min\{{a}_{{\stackrel{~}{r}}_{a,b,c,d}^{i,j,k,l},\text{}0\}}{\delta}_{{a}^{\mathrm{\prime}},0}{\delta}_{a,0}{{W}_{a}{a}_{\stackrel{~}{r}a,b,c,d}^{i,j,k,l}{\mathrm{\partial}}_{\stackrel{~}{r}}{h}_{{a}^{\mathrm{\prime}}}}_{{\stackrel{~}{r}}_{a}})\\ & +\frac{2}{\mathrm{\Delta}{\mathrm{\Theta}}^{j}}{\delta}_{{a}^{\mathrm{\prime}},a}{\delta}_{{b}^{\mathrm{\prime}},b}{\delta}_{{d}^{\mathrm{\prime}},d}{W}_{a,c,d}({{max\{{{a}_{{\stackrel{~}{\mathrm{\Theta}}}_{a,c,d}}^{i,j,k,l}}_{\stackrel{~}{\mathrm{\Theta}}=1},0\}{h}_{{b}^{\mathrm{\prime}}}}_{\stackrel{~}{\mathrm{\Theta}}=1}{h}_{b}}_{\stackrel{~}{\mathrm{\Theta}}=1}\\ & {{min\{{{a}_{{\stackrel{~}{\mathrm{\Theta}}}_{a,c,d}}^{i,j,k,l}}_{\stackrel{~}{\mathrm{\Theta}}=1},0\}{h}_{{b}^{\mathrm{\prime}}}}_{\stackrel{~}{\mathrm{\Theta}}=1}{h}_{b}}_{\stackrel{~}{\mathrm{\Theta}}=1}{{W}_{b}\text{}{a}_{{\stackrel{~}{\mathrm{\Theta}}}_{a,b,c,d}}^{i,j,k,l}{\mathrm{\partial}}_{\stackrel{~}{\mathrm{\Theta}}}{h}_{{b}^{\mathrm{\prime}}}}_{{\stackrel{~}{\mathrm{\Theta}}}_{b}})\\ & +\frac{2}{\mathrm{\Delta}{\mu}^{k}}{\delta}_{{a}^{\mathrm{\prime}},a}{\delta}_{{b}^{\mathrm{\prime}},b}{\delta}_{{d}^{\mathrm{\prime}},d}{W}_{a,b,d}{a}_{{\stackrel{~}{\mu}}_{a,b,c,d}}^{i,j,k,l}({\delta}_{{c}^{\mathrm{\prime}},{n}_{c}1}{\delta}_{c,{n}_{c}1}{W}_{c}{\mathrm{\partial}}_{\stackrel{~}{\mu}}{h}_{{c}^{\mathrm{\prime}}}{\stackrel{~}{\mu}}_{c})\\ & \phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}\frac{2}{\mathrm{\Delta}{\phi}^{l}}{\delta}_{{a}^{\mathrm{\prime}},a}{\delta}_{{b}^{\mathrm{\prime}},b}{\delta}_{{c}^{\mathrm{\prime}},c}{W}_{a,b,c}{a}_{\stackrel{~}{\phi}a,b,c,d}^{i,j,k,l}({\delta}_{{d}^{\mathrm{\prime}},0}{\delta}_{d,0}+{W}_{d}{\mathrm{\partial}}_{\stackrel{~}{\phi}}{h}_{{d}^{\mathrm{\prime}}{\stackrel{~}{\phi}}_{{}_{d}}}).\end{array}$$(B.10) $$\begin{array}{r}{b}_{{\alpha}^{\mathrm{\prime}}}^{i,j,k,l}={W}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}\text{}{\stackrel{~}{\eta}}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}}}^{i,j}\\ +\frac{2}{\mathrm{\Delta}{r}^{i}}{W}_{{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}(max\{{a}_{{\stackrel{~}{r}}_{{n}_{a}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}}^{i1,j,k,l},0\}{\stackrel{~}{I}}_{{n}_{a}1,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i1,j,k,l}{\delta}_{{a}^{\mathrm{\prime}},0}\\ min\{{a}_{{\stackrel{~}{r}}_{0},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i+1,j,k,l},0\}{\stackrel{~}{I}}_{0,{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i+1,j,k,l}{\delta}_{{a}^{\mathrm{\prime}},{n}_{a}1})\\ +\frac{2}{\mathrm{\Delta}{\mathrm{\Theta}}^{j}}{W}_{{a}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}({{max\{{a}_{{\stackrel{~}{\mathrm{\Theta}}}_{a}{}^{\mathrm{\prime}},c{}^{\mathrm{\prime}},d{}^{\mathrm{\prime}}}^{i,j,1,k,l}{}_{\stackrel{~}{\mathrm{\Theta}}=1},0\}{h}_{{b}^{\mathrm{\prime}}}}_{\stackrel{~}{\mathrm{\Theta}}=1}\sum _{b}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},b,{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j1,k,l}{h}_{b}}_{\stackrel{~}{\mathrm{\Theta}}=1}\\ {{min\{{{a}_{\stackrel{~}{\mathrm{\Theta}}a{}^{\mathrm{\prime}}c{}^{\mathrm{\prime}}d{}^{\mathrm{\prime}}}^{i,j+1,k,l}}_{\stackrel{~}{\mathrm{\Theta}}=1},0\}{h}_{{b}^{\mathrm{\prime}}}}_{\stackrel{~}{\mathrm{\Theta}}=1}\sum _{b}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},b,{c}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}^{i,j+1,k,l}{h}_{b}}_{\stackrel{~}{\mathrm{\Theta}}=1})\\ +\frac{2}{\mathrm{\Delta}{\mu}^{k}}{W}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{d}^{\mathrm{\prime}}}{a}_{{\stackrel{~}{\mu}}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{n}_{c}1,{d}^{\mathrm{\prime}}}}^{i,j,k1,l}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{n}_{c}1,{d}^{\mathrm{\prime}}}^{i,j,k1,l}{\delta}_{{c}^{\mathrm{\prime}},0}\\ \frac{2}{\mathrm{\Delta}{\phi}^{l}}{W}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}}}{a}_{{\stackrel{~}{\phi}}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},0}}^{i,j,k,l+1,}{\stackrel{~}{I}}_{{a}^{\mathrm{\prime}},{b}^{\mathrm{\prime}},{c}^{\mathrm{\prime}},0}^{i,j,k,l+1}{\delta}_{{d}^{\mathrm{\prime}},{n}_{d}1}.\end{array}$$(B.11)
Appendix C Raytracing module
The SEDs and intensity maps from the DGFEM code, shown in Sect 5, were computed with the help of the raytracing 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 ≫ R_{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,}\star}+{f}_{v}^{\text{obs,env}}.$$(C.1)
If we assume the star to be an unresolved blackbody 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_{*}/d)^{2}, $$\begin{array}{rl}& {f}_{v}^{\text{obs,}\star}(i,d)=\pi {\left(\frac{{R}_{\star}}{d}\right)}^{2}{B}_{v}\left({T}_{\star}\right)\text{}\mathrm{exp}\text{}\{\tau (i)\}\\ & \text{with}\tau (i)={\int}_{{R}_{\text{in}}}^{{R}_{\text{out}}}{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}(r,i)\text{}dr\end{array}$$(C.2)
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 × 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). 
To compute the envelope flux and intensity maps, we define an image plane $(\hat{x},\hat{y})$, at a distance d ≫ R_{out} 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 coordinates system $(\hat{r},\hat{\mathrm{\Theta}},\hat{\mathrm{\Phi}})$. 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_{out}, the flux inside the image can formally be written, $$\begin{array}{rl}{f}_{v}^{\text{obs,env}}(i)& ={\displaystyle \underset{\frac{L}{2}}{\overset{\frac{L}{2}}{\int}}}{\displaystyle \underset{\frac{L}{2}}{\overset{\frac{L}{2}}{\int}}}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}(x,y,\hat{r})\frac{dx\text{}dy}{{d}^{2}},\\ \text{or}\phantom{\rule{1em}{0ex}}{f}_{v}^{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{e}\mathrm{n}\mathrm{v}}(i)& ={\displaystyle \underset{0}{\overset{2\pi}{\int}}}{\displaystyle \underset{0}{\overset{R}{\int}}}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}(r,\omega ,\hat{r})\frac{r\text{}d\omega \text{}dr}{{d}^{2}}.\end{array}$$(C.3)
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 × N (N_{r} × N_{ω}) pixels and the flux in the image can then be rewritten, $$\begin{array}{rl}{f}_{v}^{\text{obs,env}}(i)& \approx \sum _{i=0}^{N1}\sum _{j=0}^{N1}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({x}_{i},{y}_{j},\hat{r})\frac{\mathrm{\Delta}{x}_{i}\text{}\mathrm{\Delta}{y}_{i}}{{d}^{2}},\\ \text{or}\phantom{\rule{1em}{0ex}}{f}_{v}^{\text{obs,env}}(i)& \approx \sum _{i=0}^{{N}_{r}}\sum _{j=0}^{{N}_{\omega}}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({r}_{i},{\omega}_{j},\hat{r})\frac{{r}_{i}\text{}\mathrm{\Delta}{\omega}_{j}\text{}\mathrm{\Delta}{r}_{i}}{{d}^{2}},\end{array}$$(C.4)
with Δx_{i} Δy_{i} (r_{i} Δw_{j} Δr_{i};) the pixel size of the square (circular) image. We note that the circular image is particularly wellsuited 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 ${\mathit{r}}_{0}={x}_{i}\hat{x}+{y}_{j}\hat{y}+d\hat{r}({\mathit{r}}_{0}={r}_{i}\mathrm{sin}{\omega}_{j}\hat{x}+{r}_{i}\mathrm{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, $$\begin{array}{rl}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({\mathit{r}}_{0},\hat{r})& ={\displaystyle \underset{{s}_{min}}{\overset{{s}_{max}}{\int}}}{\eta}_{v}\mathrm{exp}\{{\tau}_{v}(s)\}\text{}d\text{}s,\\ \text{with}\phantom{\rule{1em}{0ex}}{\tau}_{v}(s)& ={\displaystyle \underset{{s}_{min}}{\overset{s}{\int}}}{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}d\text{}{s}^{\mathrm{\prime}}.\end{array}$$(C.5)
We define s as to be the distance from the pixel centre r_{0} to a given point along the ray. The quantities η_{v} and ${\kappa}_{v}^{\text{ext}}$ are the emissivity and the extinction coefficient, respectively, as defined in Eq. (1). The points s_{min} and s_{max} correspond to the two intersections of the ray with the sphere of radius R_{out}.
In practice, η_{v} and ${\kappa}_{v}^{\text{ext}}$ are defined on a discrete grid and the previous integral can be rewritten as, $${I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}}({\mathit{r}}_{0},\hat{r})={\displaystyle \sum _{i=0}^{n2}}{\displaystyle \underset{{s}_{i}}{\overset{{s}_{i+1}}{\int}}}{\eta}_{v}\text{}\mathrm{exp}\text{}\{{\tau}_{v}(s)\}\text{}d\text{}s.$$(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_{min} and s_{n}_{−1} = s_{max}. The coordinates of all intersections can be computed, in the Cartesian coordinate system. We can express $\mathrm{\Delta}{I}_{v}^{i}$, the contribution to the total intensity ${I}_{v}^{\mathrm{env}}(x,y,\hat{r})$ of each portion between two consecutive intersections as, $$\begin{array}{rl}{I}_{v}^{\mathrm{e}\mathrm{n}\mathrm{v}v}({\mathit{r}}_{0},\hat{r})& =\sum _{i=0}^{n2}\mathrm{exp}({\tau}_{v}^{i})\mathrm{\Delta}{I}_{v}^{i}\\ \text{with}\mathrm{\Delta}{I}_{v}^{i}& ={\displaystyle \underset{{s}_{i}}{\overset{{s}_{i+1}}{\int}}}{\eta}_{v}\mathrm{exp}({\displaystyle \underset{{s}_{i}}{\overset{s}{\int}}}{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}d{s}^{\mathrm{\prime}})\text{}d\text{}s,\\ \text{and}{\tau}_{v}^{i}& ={\displaystyle \underset{{s}_{min}}{\overset{{s}_{i}}{\int}}}{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}d\text{}{s}^{\mathrm{\prime}}.\end{array}$$(C.7)
Following Olson et al. (1986), we assume that η_{v} and ${\kappa}_{v}^{\text{ext}}$ are linear functions between two consecutive intersections. Each contribution $\mathrm{\Delta}{I}_{v}^{i}$ is given by, $$\begin{array}{rl}\mathrm{\Delta}{I}_{v}^{i}& \approx (1\mathrm{exp}\{\mathrm{\Delta}{\tau}_{v}^{i}\}\beta )S{\text{}}_{v}\left({s}_{i}\right)+\beta \text{}S{\text{}}_{v}\left({s}_{i+1}\right),\\ \text{with}\beta & =\frac{\mathrm{\Delta}{\tau}_{v}^{i}1+\mathrm{exp}\{\mathrm{\Delta}{\tau}_{v}^{i}\}}{\mathrm{\Delta}{\tau}_{v}^{i}},\\ \mathrm{\Delta}{\tau}_{v}^{i}& =\frac{{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left({s}_{i}\right)+{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left({s}_{i+1}\right)}{2}({s}_{i+1}{s}_{i}),\end{array}$$(C.8)
and ${S}_{v}={\eta}_{v}/{\kappa}_{v}^{\text{ext}}$. The values of S_{v} and ${\kappa}_{v}^{\text{ext}}$ at the intersections s_{i} and s_{i}_{+1} are estimated by linear interpolation from the gridadjacent values (red dots in Fig. C.1). The optical depth, ${\tau}_{v}^{i}$, can be computed recursively, $${\tau}_{v}^{i+1}={\tau}_{v}^{i}+{\displaystyle \underset{{s}_{i}}{\overset{{s}_{i+1}}{\int}}}{\kappa}_{v}^{\mathrm{e}\mathrm{x}\mathrm{t}}\text{}d\text{}{s}^{\mathrm{\prime}}={\tau}_{v}^{i}+\mathrm{\Delta}{\tau}_{v}^{i},$$(C.9)
with ${\tau}_{v}^{0}=0$ and $\mathrm{\Delta}{\tau}_{v}^{i}$ defined in Eq. (C.8).
References
 Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering by a Sphere (John Wiley & Sons, Ltd), 82 [Google Scholar]
 Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
 Cockburn, B. 2003, Z. Angew. Math. Mech., 83, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, X., & Li, B. Q. 2005, JQSRT, 96, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC3D: a multipurpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
 Dykema, P. G., Klein, R. I., & Castor, J. I. 1996, ApJ, 457, 892 [NASA ADS] [CrossRef] [Google Scholar]
 Frisken, S., & Perry, R. 2002, J. Graph. Tools, 7, 1 [CrossRef] [Google Scholar]
 Gravity Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hesthaven, J. S., & Warburton, T. 2007, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (Springer Science & Business Media), 19 [Google Scholar]
 Ivezic, Z., Groenewegen, M. A. T., Men’shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] [CrossRef] [Google Scholar]
 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) [Google Scholar]
 Kitzmann, D., Bolte, J., & Patzer, A. B. C. 2016, A&A, 595, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kopriva, D. A., & Gassner, G. 2010, J. Sci. Comput., 44, 136 [CrossRef] [Google Scholar]
 Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Machorro, E. 2007, J. Computat. Phys., 223, 67 [CrossRef] [Google Scholar]
 Mihalas, D., & WeibelMihalas, B. 1999, Foundations of Radiation Hydrodynamics, Dover Books on Physics (Dover Publications) [Google Scholar]
 Nair, R. D., Levy, M. N., & Lauritzen, P. H. 2011, Emerging Numerical Methods for Atmospheric Modeling (Berlin, Heidelberg: Springer Berlin Heidelberg), 251 [Google Scholar]
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perdigon, J., Niccolini, G., & Faurobert, M. 2021, A&A, 653, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
 Reed, W. H., & Hill, T. R. 1973, Triangular mesh methods for the neutron transport equation, Tech. rep. (North Mexico, USA: Los Alamos Scientific Laboratory) [Google Scholar]
 Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63 [CrossRef] [Google Scholar]
 Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolf, S. 2003, Comput. Phys. Commun., 150, 99 [CrossRef] [Google Scholar]
Available at http://faculty.washington.edu/ivezic/dusty_web/
Available at https://www.ita.uniheidelberg.de/~dullemond/software/radmc3d/
All Tables
All Figures
Fig. 1 Illustration showing the coordinate system. We note that $(\hat{r},\hat{\mathrm{\Theta}},\hat{\mathrm{\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 Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates. 

In the text 
Fig. 2 Representation of the direction vector Ω at the equator. The dotted arrows represent the symmetric of Ω with respect to the equatorial plane $(\hat{r},\hat{\mathrm{\Phi}})$. The left and right 2D views allow the corresponding direction angles to be seen. 

In the text 
Fig. 3 Example of the sparse structure of $\mathcal{A}$ with N_{r} = N_{Θ} = N_{μ} = N_{φ} = 2. The squares represent the blocks of size n × n containing nonzero values. The matrix has bands corresponding to the coupling blocks ${\mathcal{A}}^{\alpha ,\alpha},{\mathcal{A}}^{\alpha ,\alpha +1},{\mathcal{A}}^{\alpha ,\alpha {N}_{\phi}},{\mathcal{A}}^{\alpha ,\alpha \pm {N}_{\mu}{N}_{\phi}}$, and ${\mathcal{A}}^{\alpha ,\alpha \pm {N}_{\mathrm{\Theta}}{N}_{\mu}{N}_{\phi}}$. 

In the text 
Fig. 4 Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with ${\tau}_{{v}_{0}}=1$ (blue curve) and ${\tau}_{{v}_{0}}={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. 

In the text 
Fig. 5 Temperature profiles for the axissymmetric envelope with ${\tau}_{{v}_{0}}={10}^{1}$ (blue curve) and ${\tau}_{{v}_{0}}={10}^{2}$ (orange curve), in the in the disc midplane (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 RADMC3D. The lower panels show the relative differences between the two codes. 

In the text 
Fig. 6 SED profiles for the axissymmetric envelope with ${\tau}_{{v}_{0}}={10}^{1}$ (blue curve) and ${\tau}_{{v}_{0}}={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 RADMC3D. The lower panels show the relative differences between the two codes. 

In the text 
Fig. 7 Images at λ = 2.3, 4.5, and 12.1 μm (left, middle, and right panels, respectively) of the 10 AU inner regions of the axissymmetric envelope (Sect. 5.2) and computed from the DGFEM solution with the raytracing 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}_{v}^{\text{env}}$, inside each pixel. The solid white lines show the isocontours of ${I}_{v}^{\text{env}}$. 

In the text 
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 RADMC3D 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. 

In the text 
Fig. 9 L^{2} norm of the relative differences in temperature between this study and RADMC3D versus the number of elements in the computational grid (N = N_{r} × N_{Θ} × N_{μ} × N_{φ}. The L^{2} norm was computed for the benchmark problem presented in Sec. 5.2. 

In the text 
Fig. 10 Emerging SEDs of RADMC3D (solid black lines) and of the DGFEM code, computed from the postprocessing raytracing 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 (τ = 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. 

In the text 
Fig. A.1 Geometry of a ray in the inner cavity. Point B is the opposite point of A, in the direction Ω. 

In the text 
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, Θ) and (μ, φ) 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}^{\ast}$ at the interface, along the Θ coordinate. 

In the text 
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 × 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). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.