Open Access
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/0004-6361/202244322
Published online 09 October 2024

© The Authors 2024

Licence Creative CommonsOpen 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 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.

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 Rin. The envelope spans from Rin up to the outer radius Rout. 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 θ, φ] ∈ [Rin, Rout] × (0, π/2] × [−1, 1] × [0, π]1. The radiation field is described by the time-independent radiation transfer equation (for a thorough derivation see e.g. Pomraning 1973, II-5), Ω.Iv+κvextIv=ηv, with Ω.Iv=μrIv+1μ2cosφrΘIv+1μ2rμIvcotΘ1μ2sinφrφIv.$\[\begin{array}{r}\boldsymbol{\Omega}. \nabla I_v+\kappa_v^{\mathrm{ext}} I_v=\eta_v, \\\text { with } \boldsymbol{\Omega}. \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{array}\]$(1)

The subscript v denotes the frequency dependence, Iv is the specific intensity, and Ω.∇ is the transport operator, which corresponds to a spatial derivative in the direction Ω. We note that κvext=κvabs+κvsca$\[\kappa_{v}^{\text {ext }}=\kappa_{v}^{\mathrm{abs}}+\kappa_{v}^{\text {sca }}\]$ is the extinction coefficient, with κvabs$\[\kappa_{v}^{\mathrm{abs}}\]$ and κvsca$\[\kappa_{v}^{\text {sca }}\]$ being the absorption and scattering coefficients, respectively. For dust species, we generally express these coefficients as κvabs,sca=Cvabs,scan$\[\kappa_{v}^{\text {abs,sca }}=C_{v}^{\text {abs,sca }} n\]$, where Cvabs,sca$\[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 ηv includes the thermal emission (Kirchoff-Planck law), proportional to the Planck function, Bv, 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 ηv=κvabsBv(T)+κvscaJv$\[\eta_{v}=\kappa_{v}^{\text {abs }} B_{v}(T)+\kappa_{v}^{\text {sca }} J_{v}\]$, with Jv, the mean specific intensity, which is the zeroth-order angular moment of the specific intensity Iv, Jv=14π4πIv dΩ=12π0πdφ11Ivd μ.$\[J_v=\frac{1}{4 \pi} \int_{4 \pi} I_v ~d \boldsymbol{\Omega}=\frac{1}{2 \pi} \int_0^\pi d \varphi \int_{-1}^1 I_v d ~\mu.\]$(2)

Additionally, the radiation field and the dust temperature are coupled via the equation of radiative equilibrium, 4π0κvabs Bv (T) dv=4π0κvabs Jv dv.$\[4 \pi \int_0^{\infty} \kappa_v^{\mathrm{abs}} ~B_v~(T) ~d v=4 \pi \int_0^{\infty} \kappa_v^{\mathrm{abs}} ~J_v ~d v.\]$(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), Iv=Iv+Ivenv$\[I_{v}=I_{v}^{\star}+I_{v}^{\text {env }}\]$, where Iv$\[I_{v}^{\star}\]$ is the direct unprocessed stellar radiation field attenuated by the circumstellar extinction, and Ivenv$\[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, {ΩIv+κvext Iv=0Ω.Ivenv +κvevt Ivenv =ηv=κvabs Bv+κvsca (Jv+Jvenv ),$\[\left\{\begin{array}{l}\boldsymbol{\Omega} \cdot \nabla I_v^{\star}+\kappa_v^{\text {ext }} I_v^{\star}=0 \\\boldsymbol{\Omega}. \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 }}\left(J_v^{\star}+J_v^{\text {env }}\right),\end{array}\right.\]$(4)

with Jv+Jvenv=Jv$\[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 right-hand side of the first Eq. (4). However, having it in the second equation is advantageous if κvabs,κvsca$\[\kappa_{v}^{\text {abs }}, \kappa_{v}^{\text {sca }}\]$, and κvext$\[\kappa_{v}^{\text {ext }}\]$ are independent of temperature. In such cases, the equations decouple and Iv$\[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, Iv={Γv (R, μ)exp(0s(Rin )κvext(s)ds), if μμ1,0, otherwise. $\[I_v^{\star}= \begin{cases}\Gamma_v^{\star}~\left(R_{\star}, ~\mu\right) \exp \left(-\int_0^{s\left(R_{\text {in }}\right)} \kappa_v^{\operatorname{ext}}\left(s^{\prime}\right) d s^{\prime}\right), & \text { if } \mu_{\star} \leq \mu \leq 1, \\ 0, & \text { otherwise. }\end{cases}\]$(5)

The incident stellar radiation field is Γν(R,μ)$\[\Gamma_{\nu}^{\star}\left(R_{\star}, \mu\right)\]$. 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(Rin) is the distance between a given point r and Rin, 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 rR*, 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, Jv14 (Rr)2 Bv(T) exp(Rin rκvext(r,Θ)dr).$\[J_v^{\star} \approx \frac{1}{4}~\left(\frac{R_{\star}}{r}\right)^2 ~B_v\left(T_{\star}\right) ~\exp \left(-\int_{R_{\text {in }}}^r \kappa_v^{\operatorname{ext}}\left(r^{\prime}, \Theta\right) d r^{\prime}\right).\]$(6)

In general, the integral in Eq. (6) can be carried out numerically, providing the stellar source term Jv$\[J_{v}^{\star}\]$ for Eq. (4).

To complete the description of the problem, we specify the boundary conditions for Ivenv$\[I_{v}^{\text {env }}\]$. We do so by prescribing the incident intensity Γvenv(rs,Ω)$\[\Gamma_{v}^{\text {env }}\left(\boldsymbol{r}_{s}, \boldsymbol{\Omega}\right)\]$ upon the surface of the domain D located at rs, Ivenv(rs, Ω)=Γvenv(rs, Ω),s^.Ω<0,$\[I_v^{\mathrm{env}}\left(\boldsymbol{r}_s, ~\boldsymbol{\Omega}\right)=\Gamma_v^{\mathrm{env}}\left(\boldsymbol{r}_s, ~\boldsymbol{\Omega}\right), \forall \hat{s}. \boldsymbol{\Omega}<0,\]$(7)

where ŝ is the unit vector normal to the surface of D, pointing outside of the domain. At the inner radius, r=Rin,s^=$\[r=R_{\mathrm{in}}, \hat{s}=\]$ r^$\[-\hat{r}\]$, and the incident radiation, in a given direction Ω=μr^+$\[\boldsymbol{\Omega}=\mu \hat{r}+\]$ 1μ2 cos φ Θ^+1μ2 sin φ Φ^$\[\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, Ivenv (Rin,Θ, μ, φ)=Ivenv (Rin,Θ, μ, φ), 0<μ<μ,$\[I_v^{\text {env }}\left(R_{\mathrm{in}}, \Theta, ~\mu, ~\varphi\right)=I_v^{\text {env }}\left(R_{\mathrm{in}}, \Theta^{\prime}, ~\mu^{\prime}, ~\varphi^{\prime}\right), ~\forall 0<\mu<\mu_{\star},\]$(8)

with (Rin, Θ′, μ′, φ′) being the coordinates of the opposite point. Their derivation is given in Appendix A. On the outer edge, r=$\[r=\]$ Rout,s^=r^$\[R_{\text {out }}, \hat{s}=\hat{r}\]$, and we assume that there is no incident radiation upon the surface, Ivenv (Rout ,Θ, μ, φ)=0, μ<0.$\[I_v^{\text {env }}\left(R_{\text {out }}, \Theta, ~\mu, ~\varphi\right)=0, \forall ~\mu<0.\]$(9)

At the equator, Θ=π/2,s^=Θ^$\[\Theta=\pi / 2, \hat{s}=\hat{\Theta}\]$, and the planar symmetry requires the radiation field to be, as shown in Fig. 2, Ivenv(r,π2, μ, φ)=Ivenv(r,π2, μ, πφ), φ>π2.$\[I_v^{\mathrm{env}}\left(r, \frac{\pi}{2}, ~\mu, ~\varphi\right)=I_v^{\mathrm{env}}\left(r, \frac{\pi}{2}, ~\mu, ~\pi-\varphi\right), \forall ~\varphi>\frac{\pi}{2}.\]$(10)

thumbnail Fig. 1

Illustration showing the coordinate system. We note that (r^,Θ^,Φ^)$\[(\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 Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates.

thumbnail Fig. 2

Representation of the direction vector Ω at the equator. The dotted arrows represent the symmetric of Ω with respect to the equatorial plane (r^,Φ^)$\[(\hat{r}, \hat{\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, x.F+κext I~=η~,$\[\nabla_x. \boldsymbol{F}+\kappa^{\operatorname{ext}} ~\tilde{I}=\tilde{\eta},\]$(11)

with η~=r2sinΘη$\[\tilde{\eta}=r^{2} \sin \Theta \eta\]$. We introduced the variable I~=r2sinΘI$\[\tilde{I}=r^{2} \sin \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 I~$\[\tilde{I}\]$, F=a I~=(araΘaμaφ)I~=(μ1μ2cosφ/r(1μ2)/rcotΘ1μ2sinφ/r)I~.$\[\boldsymbol{F}=\boldsymbol{a} ~\tilde{I}=\left(\begin{array}{l}a_r \\a_{\Theta} \\a_\mu \\a_{\varphi}\end{array}\right) \tilde{I}=\left(\begin{array}{c}\mu \\\sqrt{1-\mu^2} \cos \varphi / r \\\left(1-\mu^2\right) / r \\-\cot \Theta \sqrt{1-\mu^2} \sin \varphi / r\end{array}\right) \tilde{I}.\]$(12)

We decomposed the domain D into N = Nr × NΘ × Nμ × Nφ non-overlapping rectangular elements Di,j,k,l, with Nr, 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 Nr − 1, NΘ − 1, Nμ − 1 and, Nφ − 1. Inside each element Di,j,k,l, we used the nodal representation and we approximated the local solution by a four-dimensional (4D) polynomial expansion, I~hi,j,k,l(x)={a=0na1b=0nb1c=0nc1d=0nd1I~a,b,c,di,j,kha,b,c,d(x), if xDi,j,k,l,0, otherwise, $\[\tilde{I}_h^{i, j, k, l}(\boldsymbol{x})= \begin{cases}\sum_{a=0}^{n_a-1} \sum_{b=0}^{n_b-1} \sum_{c=0}^{n_c-1} \sum_{d=0}^{n_d-1} \tilde{I}_{a, b, c, d}^{i, j, k} h_{a, b, c, d}(\boldsymbol{x}), & \text { if } \boldsymbol{x} \in D^{i, j, k, l}, \\ 0, & \text { otherwise, }\end{cases}\]$(13)

with na, nb, nc, and nd, being the number of nodes inside each element Di,j,k,l, along the r, Θ, μ, and φ coordinate, respectively. We note that ha,b,c,d is the 4D Lagrange polynomial, defined as, ha,b,c,d(x)=ha (r) hb (Θ) hc (μ) hd (φ)=α=0αana1rrαrarαβ=0βbnb1ΘΘβΘbΘβγ=0γcnc1μμγμcμγδ=0δdnd1φφδφdφδ.$\[\begin{aligned}& h_{a, b, c, d}(\boldsymbol{x})=h_a~(r) ~h_b~(\Theta) ~h_c~(\mu) ~h_d~(\varphi) \\& =\prod_{\substack{\alpha=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}.\end{aligned}\]$(14)

By definition, the coefficients I~a,b,c,di,j,k,l=I~hi,j,k,l(xa,b,c,d)$\[\tilde{I}_{a, b, c, d}^{i, j, k, l}=\tilde{I}_{h}^{i, j, k, l}\left(\boldsymbol{x}_{a, b, c, d}\right)\]$ correspond to the value of I~hi,j,k,l$\[\tilde{I}_{h}^{i, j, k, l}\]$ at the nodes of coordinates xa,b,c,d = (ra, Θb, μc, φd). An example of an element Di,j,k,l with the associated nodes is shown in Fig. B.1. The global numerical approximation of the solution Ih across the domain D is formed by the sum of the N piece-wise continuous solutions inside each element, I~(x)I~h(x)=i=0Nr1j=0NΘ1k=0Nμ1l=0Nφ1I~hi,j,k,l(x).$\[\tilde{I}(\boldsymbol{x}) \approx \tilde{I}_h(\boldsymbol{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}(\boldsymbol{x}).\]$(15)

Now we shall introduce h, the residual of Eq. (11), Rh(x)=x.Fh+κext I~hη~,$\[\mathcal{R}_h(\boldsymbol{x})=\nabla_x. \boldsymbol{F}_h+\kappa^{\mathrm{ext}} ~\tilde{I}_h-\tilde{\eta},\]$(16)

with Fh=aI~h$\[\boldsymbol{F}_{h}=\boldsymbol{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, Di,j,k,lRh (x) ha,b,c,d (x) d4x=0,  Di,j,k,l, and   ha,b,c,d$\[\int_{D^{i, j, k, l}} \mathcal{R}_h~(\boldsymbol{x}) ~h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}~(\boldsymbol{x}) ~d^4 \boldsymbol{x}=0, ~\forall ~D^{i, j, k, l} \text {, and } ~\forall ~h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} \text {. }\]$(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): Di,j,k,l s^.F ha,b,c,d d3xDi,j,k,lFh.ha,b,c,d d4x+Di,j,k,l(κextI~hη~)ha,b,c,d d4x=0,  Di,j,k,l, and  ha,b,c,d$\[\begin{aligned}& \oint_{\partial D^{i, j, k, l}} ~\hat{s}. \boldsymbol{F}^* ~h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~\mathrm{d}^3 \boldsymbol{x}-\int_{D^{i, j, k, l}} \boldsymbol{F}_h. \nabla h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~\mathrm{d}^4 \boldsymbol{x} \\& +\int_{D^{i, j, k, l}}\left(\kappa^{\mathrm{ext}} \tilde{I}_h-\tilde{\eta}\right) h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~\mathrm{d}^4 \boldsymbol{x}=0, ~\forall ~D^{i, j, k, l} \text {, and } \forall ~h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} \text {. }\end{aligned}\]$(18)

The first term, in the left-hand side of Eq. (18), is a surface integral, along the boundaries of Di,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 Di,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 Lax-Friedrichs numerical flux, defined as (e.g. on the radial right element edge where s^=r^$\[\hat{s}=\hat{r}\]$), r^.F=12(Fr(I~)+Fr(I~+)|ar|(I~+I~)),$\[\hat{r}. \boldsymbol{F}^*=\frac{1}{2}\left(F_r\left(\tilde{I}^{-}\right)+F_r\left(\tilde{I}^{+}\right)-\left|a_r\right|\left(\tilde{I}^{+}-\tilde{I}^{-}\right)\right),\]$(19)

where I~$\[\tilde{I}^{-}\]$ and I~+$\[\tilde{I}^{+}\]$ denote the left and right values of I~$\[\tilde{I}\]$ at the element edge, respectively, and ar and Fr 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 = na × nb × nc × nd, relating the I~a,b,c,di,j,k,l$\[\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: AI~h=B,$\[\mathcal{A} \tilde{\boldsymbol{I}}_h=\mathcal{B},\]$(20)

with I~h$\[\tilde{\boldsymbol{I}}_{h}\]$ being the vector of size N′, which contains the solution points for the full domain D, A$\[\mathcal{A}\]$, a sparse matrix with a size of N′ × N′ coupling the elements of I~h$\[\tilde{\boldsymbol{I}}_{h}\]$, and a vector of size N′ containing the emissivity term η~$\[\tilde{\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 Jv$\[J_{v}^{\star}\]$. The radiation field of the envelope Ivenv$\[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 A$\[\mathcal{A}\]$ is often huge and solving Eq. (20) becomes numerically tedious. However, A$\[\mathcal{A}\]$ is sparse, because each solution point, Ia,b,c,di,j,k,l$\[I_{a, b, c, d}^{i, j, k, l}\]$, inside a given element, Di,j,kl, only depends on the other points inside the same element and the neighbouring ones, inside D1,j±1,k−1,l+1. This property allows us to rewrite Eq. (20) as, Ai,j,k,l I~hi,j,k,l+Ai+1,j,k,l I~hi+1,j,k,l+Ai1,j,k,l I~hi1,j,k,l+Ai,j+1,k,l I~hi,j+1,k,l+Ai,j1,k,l I~hi,j1,k,l+Ai,j,k1,l I~hi,j,k1,l+Ai,j,k,l+1 I~hi,j,k,l+1=Bi,j,k,l,  i, j, k, and l.$\[\begin{aligned}& \mathcal{A}^{i, j, k, l} ~\tilde{\boldsymbol{I}}_h^{i, j, k, l}+\mathcal{A}^{i+1, j, k, l} ~\tilde{\boldsymbol{I}}_h^{i+1, j, k, l}+\mathcal{A}^{i-1, j, k, l} ~\tilde{\boldsymbol{I}}_h^{i-1, j, k, l} \\& +\mathcal{A}^{i, j+1, k, l} ~\tilde{\boldsymbol{I}}_h^{i, j+1, k, l}+\mathcal{A}^{i, j-1, k, l} ~\tilde{\boldsymbol{I}}_h^{i, j-1, k, l}+\mathcal{A}^{i, j, k-1, l} ~\tilde{\boldsymbol{I}}_h^{i, j, k-1, l} \\& +\mathcal{A}^{i, j, k, l+1} ~\tilde{\boldsymbol{I}}_h^{i, j, k, l+1}=\mathcal{B}^{i, j, k, l}, ~\forall ~i, ~j, ~k, \text { and } l.\end{aligned}\]$(21)

We note that A$\[\mathcal{A}\]$i,j,k,l are the diagonal blocks of A$\[\mathcal{A}\]$, with a size of n × n while A$\[\mathcal{A}\]$1,j±1,k−1,l+1 (with a size of n × n) are the only nonzero, non-diagonal blocks. Furthermore, the elements Di,j,k+1,1 and Di,j,k,l−1 do not contribute because of the expression of the Lax-Friedrichs numerical flux, Eq. (19), with aμ ≥ 0 and aφ ≤ 0, ∀xD. We note that I~hi,j,k,l$\[\tilde{I}_{h}^{i, j, k, l}\]$ and ℬi,j,k,l are the sub-vector of I~h$\[\tilde{\boldsymbol{I}}_{h}\]$ and ℬ, respectively, with a length of n, containing the local points in Di,j,k,l.

To put A$\[\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 A$\[\mathcal{A}\]$ with the associated blocks is displayed in Fig. 3. The formulation Eq. (21) avoids the storage and computation of the entire matrix A$\[\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 & 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 Λ-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 Jv$\[J_{v}^{\star}\]$ with Eq. (6) and set [Ivenv]n=0=0$\[\left[I_{v}^{\text {env }}\right]^{n=0}=0\]$, as an initial condition. This allowed us to compute the initial temperature profile T0, with the help of Eq. (3), where Jv0=Jv$\[J_{v}^{0}=J_{v}^{\star}\]$. Then (ii), for each frequency, we computed [IVenv]n+1$\[\left[I_{V}^{\text {env }}\right]^{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), Ai,j,k,l[I~hi,j,k,l]n+1=[Bi,j,k,l]nAi+1,j,k,l[I~hi+1,j,k,l]nAi1,j,k,l[I~hi1,j,k,l]n+1 Ai,j+1,k,l[I~hi,j+1,k,l]nAi,j1,k,l[I~hi,j1,k,l]n+1Ai,j,k1,l[I~hi,j,k1,l]n+1Ai,j,k,l+1[I~hi,j,k,l+1]n.$\[\begin{aligned}& \mathcal{A}^{i, j, k, l}\left[\tilde{\boldsymbol{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[\tilde{\boldsymbol{I}}_h^{i+1, j, k, l}\right]^n \\&\qquad -\mathcal{A}^{i-1, j, k, l}\left[\tilde{\boldsymbol{I}}_h^{i-1, j, k, l}\right]^{n+1}~-\mathcal{A}^{i, j+1, k, l}\left[\tilde{\boldsymbol{I}}_h^{i, j+1, k, l}\right]^n \\&\quad -\mathcal{A}^{i, j-1, k, l}\left[\tilde{\boldsymbol{I}}_h^{i, j-1, k, l}\right]^{n+1}-\mathcal{A}^{i, j, k-1, l}\left[\tilde{\boldsymbol{I}}_h^{i, j, k-1, l}\right]^{n+1} \\&\qquad\qquad\qquad\qquad\qquad\qquad -\mathcal{A}^{i, j, k, l+1}\left[\tilde{\boldsymbol{I}}_h^{i, j, k, l+1}\right]^n.\end{aligned}\]$(22)

We give in appendix B the expressions of A$\[\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 Jvenv$\[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 Λ-iteration as we updated the temperature together with Iv 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).

Table 1

Summary of the main results.

thumbnail Fig. 3

Example of the sparse structure of A$\[\mathcal{A}\]$ with Nr = 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 Aα,α,Aα,α+1,Aα,αNφ,Aα,α±NμNφ$\[\mathcal{A}^{\alpha, \alpha}, \mathcal{A}^{\alpha, \alpha+1}, \mathcal{A}^{\alpha, \alpha-N_{\varphi}}, \mathcal{A}^{\alpha, \alpha \pm N_{\mu} N_{\varphi}}\]$, and Aα,α±NΘNμNφ$\[\mathcal{A}^{\alpha, \alpha \pm N_{\Theta} N_{\mu} N_{\varphi}}\]$.

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 Sec. 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 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 Rin to the outer radius Rout = 1000 Rin. The inner radius was set so that its temperature would always be Tin = T(Rin) = 800 K. The number density profile, n(r), is assumed to be a power law, of the form n(r) = n0 (Rin/r)2. The number density at the base of the disc, n0, is determined by setting the value of τv0$\[\tau_{v_{0}}\]$, which corresponds to the radial optical depth integrated though the envelope, at v0 = c/λ0, with λ0 = 1 μm, τv0=Rin Rout κv0ext dr=Cv0ext Rin Rout n dr=Cv0ext n0 Rin (1Rin Rout ).$\[\tau_{v_0}=\int_{R_{\text {in }}}^{R_{\text {out }}} \kappa_{v_0}^{\text {ext }} d r=C_{v_0}^{\text {ext }} \int_{R_{\text {in }}}^{R_{\text {out }}} n ~d r=C_{v_0}^{\mathrm{ext}} ~n_0 ~R_{\text {in }}\left(1-\frac{R_{\text {in }}}{R_{\text {out }}}\right).\]$(23)

We note that Cν0ext=Cν0abs+Cν0sca$\[C_{\nu_{0}}^{\text {ext }}=C_{\nu_{0}}^{\text {abs }}+C_{\nu_{0}}^{\text {sca }}\]$ is the extinction cross-section coefficient at v0. In this test we consider τv0=1$\[\tau_{v_{0}}=1\]$ and 102 to correspond to a moderately thin and thick envelopes, respectively. The absorption and scattering cross-sections, Cvabs$\[C_{v}^{\text {abs }}\]$ and Cvsca$\[C_{v}^{\text {sca }}\]$, feature a bi-linear behaviour in log-log scaling, with a constant profile for vv0 and a power-law dependence ∝ (v/v0)α, for v ≤ v0, 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 Cvabs/Cvext$\[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 DUSTY2. Although the envelope is spherically symmetric, we used a grid of N = Nr × NΘ × Nμ × Nφ = 164 elements in our DGFEM code, with each elements containing na ×nb × nc × nd = 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 × 104 μm.

The temperature profiles, T, and the SEDs, λFλ/F with F=0Fλdλ$\[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 Θ. 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.

thumbnail Fig. 4

Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with τv0=1$\[\tau_{v_{0}}=1\]$ (blue curve) and τv0=102$\[\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 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* = 5800 K. This disc extends from the inner radius Rin = 1 AU to the outer radius Rout = 1000 AU. The density profile n(r, Θ) is assumed to be n(r,Θ)=n0rdrsinΘexp(π4(rcosΘh)2),$\[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),\]$(24)

with h=zd(rsinΘ/rd)98,rd=Rout/2$\[h=z_{d}\left(r \sin \Theta / r_{d}\right)^{\frac{9}{8}}, r_{d}=R_{\text {out }} / 2\]$, and zd = Rout/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, n0 was determined to set τv0$\[\tau_{v_{0}}\]$, the radial optical depth, integrated through the disc mid-plane (Θ = π/2), at frequency v0 = c/λ0 with λ0 = 0.55 μm, n0=τv0Cv0ext rd ln(Rout Rin ).$\[n_0=\frac{\tau_{v_0}}{C_{v_0}^{\mathrm{ext}} ~r_d ~\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 pre-computed values for Cvext$\[C_{v}^{\text {ext }}\]$ and Cvsca$\[C_{v}^{\text {sca }}\]$ is available in Pascucci et al. (2004). We considered the optically thin and thick cases τν0=101$\[\tau_{\nu_{0}}=10^{-1}\]$ and 102, respectively. The benchmarks were produced with RADMC-3D3 (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 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 η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 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 λ = 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 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 λ = 2.3, 4.5, and 12.1 μm. In general, the disc emitting inner regions (peaks in Fig. 7) are reproduced very well (ϵ(Ivenv)3%)$\[\left(\epsilon\left(I_{v}^{\text {env }}\right) \leq 3 \%\right)\]$. The biggest discrepancies occur in the wings of the peaks, where the gradient of intensity is the steepest (in logarithmic scaling).

thumbnail Fig. 5

Temperature profiles for the axis-symmetric envelope with τv0=101$\[\tau_{v_{0}}=10^{-1}\]$ (blue curve) and τv0=102$\[\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 RADMC-3D. The lower panels show the relative differences between the two codes.

thumbnail Fig. 6

SED profiles for the axis-symmetric envelope with τv0=101$\[\tau_{v_{0}}=10^{-1}\]$ (blue curve) and τv0=102$\[\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 RADMC-3D. The lower panels show the relative differences between the two codes.

thumbnail 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 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 Ivenv$\[I_{v}^{\text {env }}\]$, inside each pixel. The solid white lines show the iso-contours of Ivenv$\[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 ray-tracing procedure to accurately compute SEDs. However, in practice, we already know the intensity Iv(r = Rout, Θ, μ, φ) 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 (Rout ΔΘ 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 Bv(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 L2 norm of the relative residual of the temperature, between our code TDG and the RADMC-3D code Tref. It is defined as, L2=(V(TDGTrefTref)2r2sinΘ dr dΘVr2sinΘ dr dΘ)12,$\[L^2=\left(\frac{\int_{\mathcal{V}}\left(\frac{T_{\mathrm{DG}}-T_{\mathrm{ref}}}{T_{\mathrm{ref}}}\right)^2 r^2 \sin \Theta ~d r ~d \Theta}{\int_{\mathcal{V}} r^2 \sin \Theta ~d r ~d \Theta}\right)^{\frac{1}{2}},\]$(26)

where the integral runs across the volume V$\[\mathcal{V}\]$ of the envelope. In Fig. 9, we display the L2 norm versus the total number of elements of our grid (we would like to emphasise that the total number of elements is N = Nr × 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 (Nr = NΘ = Nμ = Nφ).

We can see that the L2 norm saturates when N ≈ 104, 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 Iv(r = Rout, Θ, μ, φ). 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, Iv(r = Rout, Θ, μ, φ) 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 Iv(r = Rout, Θ, μ, φ) 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.

thumbnail 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.

thumbnail Fig. 9

L2 norm of the relative differences in temperature between this study and RADMC-3D versus the number of elements in the computational grid (N = Nr × NΘ × Nμ × Nφ. The L2 norm was computed for the benchmark problem presented in Sec. 5.2.

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

thumbnail 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 (τ = 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 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.

Appendix A Boundary conditions on a spherical enclosed cavity

The boundary condition in Eq. (7) requires to know, for a given point A (Rin, Θ) on the inner cavity of radius Rin and in a direction Ω(μ,φ)(Ω.r^>0)$\[\boldsymbol{\Omega}(\mu, \varphi)(\boldsymbol{\Omega}. \hat{r}&gt;0)\]$, the coordinates of the opposite point B (Rin, Θ′, μ′ = 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 axis-symmetric 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.

thumbnail 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^{\prime}=-\mu.\]$(A.1)

The position rB of the opposite point B, is linked to the position rA of the point A through the relation rB=rAsAB Ω,$\[\boldsymbol{r}_B=\boldsymbol{r}_A-s_{A B} ~\boldsymbol{\Omega},\]$(A.2)

with sAB the distance between the point A and B. This length is equal to sAB = 2 μ Rin, because it corresponds to the base of the isosceles triangle OAB, with OA = OB = Rin. The direction vector, Ω, is, Ω=[μ1μ2cosφ1μ2sinφ](r^,Θ^,Φ^)=[1μ2sinφsinΘμ+cosΘ1μ2cosφcosΘμsinΘ1μ2cosφ](x^,y^,z^).$\[\begin{aligned}& \boldsymbol{\Omega}=\left[\begin{array}{cc}\mu & \\\sqrt{1-\mu^2} & \cos \varphi \\\sqrt{1-\mu^2} & \sin \varphi\end{array}\right]_{(\hat{r}, \hat{\Theta}, \hat{\Phi})} \\& =\left[\begin{array}{c}-\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{array}\right]_{(\hat{x}, \hat{y}, \hat{z})}.\end{aligned}\]$(A.3)

Eq. (A.2) is then rB=Rin[sinΘcosΦsinΘsinΦcosΘ](x^,y^,z^)=Rin[2μ1μ2sinφsinΘ(12μ2)2cosΘμ1μ2 cos φ cos Θ(12μ2)+2 sinΘμ1μ2 cos φ](x^,y^,z^).$\[\begin{aligned}\boldsymbol{r}_B & =R_{\mathrm{in}}\left[\begin{array}{c}\sin \Theta^{\prime} \cos \Phi^{\prime} \\\sin \Theta^{\prime} \sin \Phi^{\prime} \\\cos \Theta^{\prime}\end{array}\right]_{(\hat{x}, \hat{y}, \hat{z})} \\& =R_{\mathrm{in}}\left[\begin{array}{c}2 \mu \sqrt{1-\mu^2} \sin \varphi \\\sin \Theta\left(1-2 \mu^2\right)-2 \cos \Theta \mu \sqrt{1-\mu^2} \cos \varphi \\\cos \Theta\left(1-2 \mu^2\right)+2 \sin \Theta \mu \sqrt{1-\mu^2} \cos \varphi\end{array}\right]_{(\hat{x}, \hat{y}, \hat{z})}.\end{aligned}\]$(A.4)

From the z-component of Eq. (A.4) we get, cosΘ=cosΘ(12μ2)+2sinΘμ1μ2 cos φ,$\[\cos \Theta^{\prime}=\cos \Theta\left(1-2 \mu^2\right)+2 \sin \Theta \mu \sqrt{1-\mu^2} \cos \varphi,\]$(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, tanφ=ΩΦ^ΩΘ^,$\[\tan \varphi^{\prime}=\frac{\boldsymbol{\Omega} \cdot \hat{\Phi}^{\prime}}{\boldsymbol{\Omega} \cdot \hat{\Theta}^{\prime}},\]$(A.6)

with, Φ^=[sinΦcosΦ0](x^,y^,z^),Θ^=[cosΘcosΦcosΘsinΦsinΘ](x^,y^,z^),$\[\hat{\Phi}^{\prime}=\left[\begin{array}{c}-\sin \Phi^{\prime} \\\cos \Phi^{\prime} \\0\end{array}\right]_{(\hat{x}, \hat{y}, \hat{z})}, \hat{\Theta}^{\prime}=\left[\begin{array}{c}\cos \Theta^{\prime} \cos \Phi^{\prime} \\\cos \Theta^{\prime} \sin \Phi^{\prime} \\-\sin \Theta^{\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, tanφ=sinφsinΘ(12 μ2)cosφsinΘ2μ 1μ2cosΘ.$\[\tan \varphi^{\prime}=\frac{\sin \varphi \sin \Theta}{\left(1-2 ~\mu^2\right) \cos \varphi \sin \Theta-2 \mu ~\sqrt{1-\mu^2} \cos \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{\Theta}^{\prime}\]$ and Φ^$\[\hat{\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 Di,j,k,l, for each coordinate (r, Θ, μ, φ). For the r, μ, φ 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 Fh$\[\boldsymbol{F}_{h}^{*}\]$ at the element edges. The Gauss-Lobatto quadrature can exactly integrate polynomials up to degree 2 n − 3, with n the number of nodes. For the coordinate Θ, we cannot use the Gauss-Lobatto 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 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 Di,j,k,l with the chosen nodes is shown in Fig. B.1.

thumbnail Fig. B.1

Example of an element Di,j,k,l (dashed lines) using na = nc = nd = 3 and nb = 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 Fh$\[F_{h}^{*}\]$ 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), Di,j,k,l(κextI~hη~)ha,b,c,d d4x=Δxi,j,k,l1611(κextI~hi,j,k,lη~i,j)ha,b,c,d d4x~=Δxi,j,k,l16a,b,c,dWa,b,c,d(κa,bexti,jI~a,b,c,di,j,k,lη~a,bi,j)ha,b,c,d(x~a,b,c,d).$\[\begin{array}{r}\int_{D^{i, j, k, l}}\left(\kappa^{\mathrm{ext}} \tilde{I}_h-\tilde{\eta}\right) h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~d^4 \boldsymbol{x} \\=\frac{\Delta x^{i, j, k, l}}{16} \int_{-1}^1\left(\kappa^{\mathrm{ext}} \tilde{I}_h^{i, j, k, l}-\tilde{\eta}^{i, j}\right) h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~d^4 \tilde{\boldsymbol{x}} \\=\frac{\Delta x^{i, j, k, l}}{16} \sum_{a, b, c, d} W_{a, b, c, d}\left(\kappa_{a, b}^{\mathrm{ext} ^{i, j}} \tilde{I}_{a, b, c, d}^{i, j, k, l}-\tilde{\eta}_{a, b}^{i, j}\right) h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\left(\tilde{\boldsymbol{x}}_{a, b, c, d}\right).\end{array}\]$(B.1)

We note that Δxi,j,k,l = Δri, Δ Θj Δμk Δφl is the 4D volume of the element Di,j,k,l For integration, we rather use the local coordinates (r~,Θ~,μ~,φ~)$\[(\tilde{r}, \tilde{\Theta}, \tilde{\mu}, \tilde{\varphi})\]$, defined as (e.g for the r coordinate), r~=2Δri(rri+1/2)$\[\tilde{r}=\frac{2}{\Delta r^i}\left(r-r^{i+1 / 2}\right) \text {, }\]$(B.2)

with Δri, the element width along the coordinate r and ri+1/2, the radial coordinate of the centre of the element. The same expression holds for the other coordinates. The quantities Wa,b,c,d=Wr(r~a)WΘ(Θ~b)Wμ(μ~c)Wφ(φ~d)$\[W_{a, b, c, d}=W_{r}\left(\tilde{r}_{a}\right) W_{\Theta}\left(\tilde{\Theta}_{b}\right) W_{\mu}\left(\tilde{\mu}_{c}\right) W_{\varphi}\left(\tilde{\varphi}_{d}\right)\]$ are the weights associated with the different quadrature in each direction. Finally, ha,b,c,d(x~a,b,c,d)$\[h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\left(\tilde{\boldsymbol{x}}_{a, b, c, d}\right)\]$ is the 4D Lagrange polynomials, defined in Eq. (14), evaluated at the node x~a,b,c,d=(r~a,Θ~b,μ~c,φ~d)$\[\tilde{\boldsymbol{x}}_{a, b, c, d}=\left(\tilde{r}_{a}, \tilde{\Theta}_{b}, \tilde{\mu}_{c}, \tilde{\varphi}_{d}\right)\]$. By definition of the Lagrange polynomials, we have, ha,b,c,d(r~a,Θ~b, μ~c, φ~d)=δa,a δb,b δc,c δd,d,$\[h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\left(\tilde{r}_a, \tilde{\Theta}_b, ~\tilde{\mu}_c, ~\tilde{\varphi}_d\right)=\delta_{a^{\prime}, a} ~\delta_{b^{\prime}, b} ~\delta_{c^{\prime}, c} ~\delta_{d^{\prime}, d},\]$(B.3)

with δa′,a the usual delta Kronecker. The other volume term in Eq. (18) is expressed as, Di,j,k,lFh.xha,b,c,d d4x=Di,j,k,lar I~h rha,b,c,d d4x+Di,j,k,laΘ I~h Θha,b,c,d d4x+Di,j,k,laμ I~h μha,b,c,dd4x+Di,j,k,laφI~h φha,b,c,dd4x=Δxi,j,k,l8a,b,c,dWa,b,c,d I~a,b,c,di,j,k,l(ar~a,b,c,di,j,k,l r~ha,b,c,dl|x~a,b,c,dΔri+aΘ~a,b,c,di,j,k,l Θ~ha,b,c,d|x~a,b,c,dΔΘj+aμ~a,b,c,di,j,k,l μ~ha,b,c,d|x~a,b,c,dΔμk+aφ~a,b,c,di,j,k,l φ~ha,b,c,d|x~a,b,c,dΔφl).$\[\begin{aligned}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad \int_{D^{i, j, k, l}} \boldsymbol{F}_h. \nabla_x h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~d^4 \boldsymbol{x} \\& =\int_{D^{i, j, k, l}} a_r ~\tilde{\boldsymbol{I}}_h ~\partial_r h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~d^4 \boldsymbol{x}+\int_{D^{i, j, k, l}} a_{\Theta} ~\tilde{\boldsymbol{I}}_h ~\partial_{\Theta} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~d^4 \boldsymbol{x} \\& +\int_{D^{i, j, k, l}} a_\mu ~\tilde{\boldsymbol{I}}_h ~\partial_\mu h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} d^4 \boldsymbol{x}+\int_{D^{i, j, k, l}} a_{\varphi} \tilde{\boldsymbol{I}}_h ~\partial_{\varphi} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} d^4 \boldsymbol{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{\left.a_{\tilde{r}_{a, b, c, d}}^{i, j, k, l} ~\partial_{\tilde{r}} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{l^{\prime}}}\right|{\tilde{x}_{a, b, c, d}}}{\Delta r^i}\right. \\& +\frac{\left.a_{\tilde{\Theta}_{a, b, c, d}}^{i, j, k, l} ~\partial_{\tilde{\Theta}} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\right|{\tilde{x}_{a, b, c, d}}}{\Delta \Theta^j}+\frac{\left.a_{\tilde{\mu}_{a, b, c, d}}^{i, j, k, l} ~\partial_{\tilde{\mu}} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\right|{\tilde{x}_{a, b, c, d}}}{\Delta \mu^k} \\&\qquad\qquad\qquad\qquad\qquad\left.+\frac{\left.a_{\tilde{\varphi} a, b, c, d}^{i, j, k, l} ~\partial_{\tilde{\varphi}} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\right|{\tilde{x}_{a, b, c, d}}}{\Delta \varphi^l}\right).\end{aligned}\]$(B.4)

In Eq. (B.4), we used the definition of the flux Eq. (12). The quantity r~ha,b,c,d|xa,b,c,d $\[\left.\partial_{\tilde{r}} h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\right|x_{a, b, c, d}\]$ is the partial derivative of the Lagrange polynomial, with respect to the coordinate r, evaluated at the node x~a,b,c,d$\[\tilde{\boldsymbol{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 s^=±r^$\[\hat{s}= \pm \hat{r}\]$ for the radial right and left surfaces (r~=±1)$\[(\tilde{r}= \pm 1)\]$, respectively, Di,j,k,l[Frha,b,c,d]r~=1r~=1dΘ dμ dφ=Δxi,j,k,l8 ΔriWb,c,d[Fr~(r~,Θ~b,μ~c,φ~d)ha(r~)]r~=1r~=1.$\[\begin{array}{r}\oint_{\partial D^{i, j, k, l}}\left[F_r^* h_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}}\right]_{\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^{\prime}, c^{\prime}, d^{\prime}}\left[F_{\tilde{r}}^*\left(\tilde{r}, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right) h_{a^{\prime}}(\tilde{r})\right]_{\tilde{r}=-1}^{\tilde{r}=1}.\end{array}\]$(B.5)

At the cell boundaries, we use the upwind numerical flux, for example, at the right edge, Fr~(r~=1,Θ~b,μ~c,φ~d)=max{ar~i,j,k,l(1,Θ~b,μ~c,φ~d),0}I~hi,j,k,l(1,Θ~b,μ~c,φ~d)+min{ar~i,j,k,l(1,Θ~b,μ~c,φ~d),0}I~hi+1,j,k,l(1,Θ~b,μ~c,φ~d).$\[\begin{array}{r}F_{\tilde{r}}^*\left(\tilde{r}=1, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right)= \\\max \left\{a_{\tilde{r}}^{i, j, k, l}\left(1, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right), 0\right\} \tilde{\boldsymbol{I}}_h^{i, j, k, l}\left(1, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right) \\+\min \left\{a_{\tilde{r}}^{i, j, k, l}\left(1, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right), 0\right\} \tilde{\boldsymbol{I}}_h^{i+1, j, k, l}\left(-1, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right).\end{array}\]$(B.6)

The use of the Gauss-Lobatto quadrature further simplify the computation because the quadrature nodes include the endpoints r~=±1$\[\tilde{r}= \pm 1\]$. The numerical flux can then be expressed as [Fr~(r~,Θ~b,μ~c,φ~d)ha(r~)]r~=1r~=1=δa,na1max{ar~na1,b,c,di,j,k,l,0}I~na1,b,c,di,j,k,l+δa,na1min{ar~na1,b,c,di,j,k,l,0}I~0,b,c,di+1,j,k,lδa,0max{ar~0,b,c,di,j,k,l,0}I^na1,b,c,di1,j,k,lδa,0min{ar~0,b,c,di,j,k,l,0}I^0,b,c,di,j,k,l.$\[\begin{aligned}&\qquad {\left[F_{\tilde{r}}^*\left(\tilde{r}, \tilde{\Theta}_{b^{\prime}}, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right) h_{a^{\prime}}(\tilde{r})\right]_{\tilde{r}=-1}^{\tilde{r}=1}=} \\& \delta_{a^{\prime}, n_a-1} \max \left\{a_{\tilde{r}_{n_a}-1, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j, k, l}, 0\right\} \tilde{I}_{n_a-1, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j, k, l} \\&\quad +\delta_{a^{\prime}, n_a-1} \min \left\{a_{\tilde{r}_{n_a}-1, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j,k,l}, 0\right\} \tilde{I}_{0, b^{\prime}, c^{\prime}, d^{\prime}}^{i+1, j, k, l}\\&\qquad -\delta_{a^{\prime}, 0} \max \left\{a_{\tilde{r}_0, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j, k,l}, 0\right\} \hat{I}_{n_a-1, b^{\prime}, c^{\prime}, d^{\prime}}^{i-1, j, k,l}\\&\qquad\quad -\delta_{a^{\prime}, 0} \min \left\{a_{\tilde{r}_0, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j, k,l}, 0\right\} \hat{I}_{0, b^{\prime}, c^{\prime}, d^{\prime}}^{i, j, k,l}.\end{aligned} \]$(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 ∀ xD, so there are no terms proportional to I~i,j,k+1,l$\[\widetilde{I}^{i, j, k+1, l}\]$ and I~i,j,k,l1$\[\tilde{I}^{i, j, k, l-1}\]$. 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), FΘ~(r~a,Θ~=1,μ~c,φ~d)=max{aΘ~i,j,k,l(r~a,1,μ~c,φ~d),0}bI~a,b,c,di,j,k,l hb(1)+min{aΘ~i,j,k,l(r~a,1,μ~c,φ~d),0}bI~a,b,c,di,j+1,k,l hb(1).$\[\begin{aligned}&\qquad\quad\qquad\qquad\qquad\qquad F_{\tilde{\Theta}}^*\left(\tilde{r}_{a^{\prime}}, \tilde{\Theta}=1, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right)= \\&\qquad \max \left\{a_{\tilde{\Theta}}^{i, j, k, l}\left(\tilde{r}_{a^{\prime}}, 1, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right), 0\right\} \sum_b \tilde{I}_{a^{\prime}, b, c^{\prime}, d^{\prime}}^{i, j, k, l} ~h_b(1) \\& +\min \left\{a_{\tilde{\Theta}}^{i, j, k, l}\left(\tilde{r}_{a^{\prime}}, 1, \tilde{\mu}_{c^{\prime}}, \tilde{\varphi}_{d^{\prime}}\right), 0\right\} \sum_b \tilde{I}_{a^{\prime}, b, c^{\prime}, d^{\prime}}^{i, j+1, k, l} ~h_b(-1).\end{aligned}\]$(B.8)

All the terms in this section can be put in the form of the system of equations, given by Eq. (22), Ai,j,k,l,I~hi,j,k,l=bi,j,k,l.$\[\mathcal{A}^{i, j, k, l,} \tilde{\boldsymbol{I}}_h^{i, j, k, l}=\boldsymbol{b}^{i, j, k, l}.\]$(B.9)

In Eq. (B.9), we replaced the RHS by bi,j,k,l because we do not need to formally write the non-diagonal matrices A$\[\mathcal{A}\]$i±1,j±1,k−1,l+1. To assemble A$\[\mathcal{A}\]$i,j,k,l, we make use of the global index α = a nbncnd + b ncnd + c nd + d. The elements of A$\[\mathcal{A}\]$i,j,k,l and bi,j,k,l are then, Aααi,j,k,l=Wa,b,c,d κa,bexti,jδa,a δb,b δc,c δd,d +2Δriδb,bδc,cδd,dWb,c,d(max{ar~a,b,c,di,j,k,l,0}δa,na1δa,na1min{ar~a,b,c,di,j,k,l, 0}δa,0δa,0Waar~a,b,c,di,j,k,lr~ha|r~a)+2ΔΘjδa,aδb,bδd,dWa,c,d(max{aΘ~a,c,di,j,k,l|Θ~=1,0}hb|Θ~=1hb|Θ~=1min{aΘ~a,c,di,j,k,l|Θ~=1,0}hb|Θ~=1hb|Θ~=1Wb aΘ~a,b,c,di,j,k,lΘ~hb|Θ~b)+2Δμkδa,aδb,bδd,dWa,b,daμ~a,b,c,di,j,k,l(δc,nc1δc,nc1Wcμ~hc|μ~c)2Δφlδa,aδb,bδc,cWa,b,caφ~a,b,c,di,j,k,l(δd,0δd,0+Wdφ~hd|φ~d).$\[\begin{aligned}&\qquad\qquad\qquad\qquad \mathcal{A}_{\alpha^{\prime} \alpha}^{i, j, k, l}=W_{a, b, c, d} ~\kappa_{a, b}^{\mathrm{ext}^{ i, j}} \delta_{a^{\prime}, a} ~\delta_{b^{\prime}, b} ~\delta_{c^{\prime}, c} ~\delta_{d^{\prime}, d} \\&~\quad +\frac{2}{\Delta r^i} \delta_{b^{\prime}, b} \delta_{c^{\prime}, c} \delta_{d^{\prime}, d} W_{b, c, d}\left(\max \left\{a_{\tilde{r}_{a, b, c, d}}^{i, j, k, l}, 0\right\} \delta_{a^{\prime}, n_a-1} \delta_{a, n_a-1}\right.\\&\qquad\qquad\qquad-\min \left\{a_{\tilde{r}_{a, b, c, d}^{i, j, k, l}, ~0\}} \delta_{a^{\prime}, 0} \delta_{a, 0}-\left.W_a a_{\tilde{r} a, b, c, d}^{i, j, k, l} \partial_{\tilde{r}} h_{a^{\prime}}\right|_{\tilde{r}_a}\right)\\& +\frac{2}{\Delta\Theta^j} \delta_{a^{\prime}, a} \delta_{b^{\prime}, b} \delta_{d^{\prime}, d} W_{a, c, d}\left(\left.\left.\max \left\{\left.a_{\tilde{\Theta}_{a, c, d}}^{i, j, k,l}\right|_{\tilde{\Theta}=1}, 0\right\} h_{b^{\prime}}\right|_{\tilde{\Theta}=1} h_b\right|_{\tilde{\Theta}=1}\right. \\- & \left.\left.\left.\min \left\{\left.a_{\tilde{\Theta}_{a, c, d}}^{i, j, k, l}\right|_{\tilde{\Theta}=-1}, 0\right\} h_{b^{\prime}}\right|_{\tilde{\Theta}=-1} h_b\right|_{\tilde{\Theta}=-1}-\left.W_b ~a_{\tilde{\Theta}_{a, b, c, d}}^{i, j, k, l}{ } \partial_{\tilde{\Theta}} h_{b^{\prime}}\right|_{\tilde{\Theta}_b}\right)\\&+\frac{2}{\Delta \mu^k} \delta_{a^{\prime}, a} \delta_{b^{\prime}, b} \delta_{d^{\prime}, d} W_{a, b, d} a_{\tilde{\mu}_{a, b, c, d}}^{i, j, k, l} \left(\delta_{c^{\prime}, n_c-1} \delta_{c, n_c-1}-W_c \partial_{\tilde{\mu}} h_{c^{\prime}} |\tilde{{\mu}}_c\right)\\&\qquad\quad-\frac{2}{\Delta \varphi^{l}} \delta_{a^{\prime}, a} \delta_{b^{\prime}, b} \delta_{c^{\prime}, c} W_{a, b, c} a_{\tilde{\varphi} a, b, c, d}^{i, j, k, l}\left(\delta_{d^{\prime}, 0} \delta_{d, 0}+W_d \partial_{\tilde{\varphi}} h_{d^{\prime} |\tilde{\varphi}_{_d}}\right).\end{aligned}\]$(B.10) bαi,j,k,l=Wa,b,c,d η~a,bi,j+2ΔriWb,c,d(max{ar~na1,b,c,di1,j,k,l,0}I~na1,b,c,di1,j,k,lδa,0min{ar~0,b,c,di+1,j,k,l,0}I~0,b,c,di+1,j,k,lδa,na1)+2ΔΘjWa,c,d(max{aΘ~a,c,di,j,1,k,l|Θ~=1,0}hb|Θ~=1bI~a,b,c,di,j1,k,lhb|Θ~=1min{aΘ~acdi,j+1,k,l|Θ~=1,0}hb|Θ~=1bI~a,b,c,di,j+1,k,lhb|Θ~=1)+2ΔμkWa,b,daμ~a,b,nc1,di,j,k1,lI~a,b,nc1,di,j,k1,lδc,02ΔφlWa,b,caφ~a,b,c,0i,j,k,l+1,I~a,b,c,0i,j,k,l+1δd,nd1.$\[\begin{aligned}b_{\alpha^{\prime}}^{i, j, k, l}=W_{a^{\prime}, b^{\prime}, c^{\prime}, d^{\prime}} ~\tilde{\eta}_{a^{\prime}, b^{\prime}}^{i, j}\\+\frac{2}{\Delta r^i} W_{b^{\prime}, c^{\prime}, d^{\prime}}\left(\max \left\{a_{\tilde{r}_{n_a-1, b^{\prime}, c^{\prime}, d^{\prime}}}^{i-1, j, k, l}, 0\right\} \tilde{I}_{n_a-1, b^{\prime}, c^{\prime}, d^{\prime}}^{i-1, j, k, l} \delta_{a^{\prime}, 0}\right.\\\left.-\min \left\{a_{\tilde{r}_0, b^{\prime}, c^{\prime}, d^{\prime}}^{i+1, j,k,l}, 0\right\} \tilde{I}_{0, b^{\prime}, c^{\prime}, d^{\prime}}^{i+1, j, k, l} \delta_{a^{\prime}, n_a-1}\right)\\+\frac{2}{\Delta \Theta^j} W_{a^{\prime}, c^{\prime}, d^{\prime}}\left(\left.\left.\max \left\{a_{\tilde{\Theta}_a{^\prime},c{^\prime},d{^\prime}}^{i, j,-1, k, l}{ }|_{\tilde{\Theta}=1}, 0\right\} h_{b^{\prime}}\right|_{\tilde{\Theta}=-1} \sum_b \tilde{I}_{a^{\prime}, b, c^{\prime}, d^{\prime}}^{i, j-1, k, l} h_b\right|_{\tilde{\Theta}=1}\right.\\\left.-\left.\left.\min \left\{\left.a_{\tilde{\Theta}a{^\prime}c{^\prime}d{^\prime}}^{i, j+1,k, l}\right|_{\tilde{\Theta}=-1}, 0\right\} h_{b^{\prime}}\right|_{\tilde{\Theta}=1} \sum_b \tilde{I}_{a^{\prime}, b, c^{\prime}, d^{\prime}}^{i, j+1, k, l} h_b\right|_{\tilde{\Theta}=-1}\right)\\+\frac{2}{\Delta \mu^k} W_{a^{\prime}, b^{\prime}, d^{\prime}} a_{\tilde{\mu}_{a^{\prime}, b^{\prime}, n_c-1, d^{\prime}}}^{i, j, k-1, l} \tilde{I}_{a^{\prime}, b^{\prime}, n_c-1, d^{\prime}}^{i, j, k-1, l} \delta_{c^{\prime}, 0}\\-\frac{2}{\Delta \varphi^l} W_{a^{\prime}, b^{\prime}, c^{\prime}} a_{\tilde{\varphi}_{a^{\prime}, b^{\prime}, c^{\prime}, 0}}^{i, j, k, l+1,} \tilde{I}_{a^{\prime}, b^{\prime}, c^{\prime}, 0}^{i, j, k, l+1} \delta_{d^{\prime}, n_d-1}.\end{aligned}\]$(B.11)

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 fvobs$\[f_{v}^{\text {obs }}\]$ that an observer receive from the object, situated at a distance dRout 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, fvobs (i,d)=fvobs, +fvobs,env .$\[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 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*/d)2, fvobs, (i,d)=π(Rd)2Bv(T) exp {τ(i)} with    τ(i)=Rin Rout κvext(r,i) dr$\[\begin{aligned}& f_v^{\text {obs, } \star}(i, d)=\pi\left(\frac{R_{\star}}{d}\right)^2 B_v\left(T_{\star}\right) ~\exp~ \{-\tau(i)\} \\& \text { with }~~~ \tau(i)=\int_{R_{\text {in }}}^{R_{\text {out }}} \kappa_v^{\mathrm{ext}}(r, i) ~d r\end{aligned}\]$(C.2)

thumbnail 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 κvext$\[\kappa_{v}^{\text {ext }}\]$ and Sv 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 (x^,y^)$\[(\hat{x}, \hat{y})\]$, at a distance dRout 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 (r^,Θ^,Φ^)$\[(\hat{r}, \hat{\Theta}, \hat{\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 dRout, the flux inside the image can formally be written, fvobs,env (i)=L2L2L2L2Ivenv(x,y,r^)dx dyd2, or fvobs,env(i)=02π0RIvenv(r,ω,r^)r dω drd2.$\[\begin{aligned}f_v^{\text {obs,env }}(i) & =\int_{-\frac{L}{2}}^{\frac{L}{2}} \int_{-\frac{L}{2}}^{\frac{L}{2}} I_v^{\mathrm{env}}(x, y, \hat{r}) \frac{d x ~d y}{d^2}, \\\text { or } \quad f_v^{\mathrm{obs}, \mathrm{env}}(i) & =\int_0^{2 \pi} \int_0^R I_v^{\mathrm{env}}(r, \omega, \hat{r}) \frac{r ~d \omega ~d r}{d^2}.\end{aligned}\]$(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 (Nr × Nω) pixels and the flux in the image can then be rewritten, fvobs,env (i)i=0N1j=0N1Ivenv(xi,yj,r^)Δxi Δyid2, or fvobs,env (i)i=0Nrj=0NωIvenv(ri,ωj,r^)ri Δωj Δrid2,$\[\begin{aligned}f_v^{\text {obs,env }}(i) & \approx \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} I_v^{\mathrm{env}}\left(x_i, y_j, \hat{r}\right) \frac{\Delta x_i ~\Delta y_i}{d^2}, \\\text { or } \quad f_v^{\text {obs,env }}(i) & \approx \sum_{i=0}^{N_r} \sum_{j=0}^{N_\omega} I_v^{\mathrm{env}}\left(r_i, \omega_j, \hat{r}\right) \frac{r_i ~\Delta \omega_j ~\Delta r_i}{d^2},\end{aligned}\]$(C.4)

with Δxi Δyi (ri Δwj Δri;) the pixel size of the square (circular) image. We note that the circular image is particularly well-suited for the computation of fvobs,env(i)$\[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 r0=xix^+yjy^+dr^(r0=risinωjx^+ricosωjy^+dr^$\[\boldsymbol{r}_{0}=x_{i} \hat{x}+y_{j} \hat{y}+d \hat{r}\left(\boldsymbol{r}_{0}=r_{i} \sin \omega_{j} \hat{x}+r_{i} \cos \omega_{j} \hat{y}+d \hat{r}\right.\]$ 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, Ivenv(r0,r^)=sminsmaxηvexp{τv(s)} d s, with τv(s)=sminsκvext d s.$\[\begin{aligned}I_v^{\mathrm{env}}\left(\boldsymbol{r}_0, \hat{r}\right) & =\int_{s_{\min }}^{s_{\max }} \eta_v \exp \left\{-\tau_v(s)\right\} ~d ~s, \\\text { with } \quad \tau_v(s) & =\int_{s_{\min }}^s \kappa_v^{\mathrm{ext}} ~d ~s^{\prime}.\end{aligned}\]$(C.5)

We define s as to be the distance from the pixel centre r0 to a given point along the ray. The quantities ηv and κvext$\[\kappa_{v}^{\text {ext }}\]$ are the emissivity and the extinction coefficient, respectively, as defined in Eq. (1). The points smin and smax correspond to the two intersections of the ray with the sphere of radius Rout.

In practice, ηv and κvext$\[\kappa_{v}^{\text {ext }}\]$ are defined on a discrete grid and the previous integral can be rewritten as, Ivenv(r0,r^)=i=0n2sisi+1ηv exp {τv(s)} d s.$\[I_v^{\mathrm{env}}\left(\boldsymbol{r}_0, \hat{r}\right)=\sum_{i=0}^{n-2} \int_{s_i}^{s_{i+1}} \eta_v ~\exp ~\left\{-\tau_v(s)\right\} ~d ~s.\]$(C.6)

The {si} are the n intersections between the ray and the grid (red crosses in Fig. C.1), with s0 = smin and sn−1 = smax. The coordinates of all intersections can be computed, in the Cartesian coordinate system. We can express ΔIvi$\[\Delta I_{v}^{i}\]$, the contribution to the total intensity Ivenv(x,y,r^)$\[I_{v}^{\mathrm{env}}(x, y, \hat{r})\]$ of each portion between two consecutive intersections as, Ivenvv(r0,r^)=i=0n2exp(τvi)ΔIvi with ΔIvi=sisi+1ηvexp(sisκvextds) d s, and τvi=sminsiκvext d s.$\[\begin{aligned}I_v^{\mathrm{env} v}\left(\boldsymbol{r}_0, \hat{r}\right) & =\sum_{i=0}^{n-2} \exp \left(-\tau_v^i\right) \Delta I_v^i \\\text { with } \Delta I_v^i & =\int_{s_i}^{s_{i+1}} \eta_v \exp \left(-\int_{s_i}^s \kappa_v^{\mathrm{ext}} d s^{\prime}\right) ~d~ s, \\\text { and } \tau_v^i & =\int_{s_{\min }}^{s_i} \kappa_v^{\mathrm{ext}} ~d ~s^{\prime}.\end{aligned}\]$(C.7)

Following Olson et al. (1986), we assume that ηv and κvext$\[\kappa_{v}^{\text {ext }}\]$ are linear functions between two consecutive intersections. Each contribution ΔIvi$\[\Delta I_{v}^{i}\]$ is given by, ΔIvi(1exp{Δτvi}β)S v(si)+β S v(si+1), with β=Δτvi1+exp{Δτvi}Δτvi,Δτvi=κvext(si)+κvext(si+1)2(si+1si),$\[\begin{aligned}\Delta I_v^i & \approx\left(1-\exp \left\{-\Delta \tau_v^i\right\}-\beta\right) S~_v\left(s_i\right)+\beta ~S~_v\left(s_{i+1}\right), \\\text { with } \beta & =\frac{\Delta \tau_v^i-1+\exp \left\{-\Delta \tau_v^i\right\}}{\Delta \tau_v^i}, \\\Delta \tau_v^i & =\frac{\kappa_v^{\mathrm{ext}}\left(s_i\right)+\kappa_v^{\mathrm{ext}}\left(s_{i+1}\right)}{2}\left(s_{i+1}-s_i\right),\end{aligned}\]$(C.8)

and Sv=ηv/κvext$\[S_{v}=\eta_{v} / \kappa_{v}^{\text {ext }}\]$. The values of Sv and κvext$\[\kappa_{v}^{\text {ext }}\]$ at the intersections si and si+1 are estimated by linear interpolation from the grid-adjacent values (red dots in Fig. C.1). The optical depth, τvi$\[\tau_{v}^{i}\]$, can be computed recursively, τvi+1=τvi+sisi+1κvext d s=τvi+Δτvi,$\[\tau_v^{i+1}=\tau_v^i+\int_{s_i}^{s_{i+1}} \kappa_v^{\mathrm{ext}} ~d~ s^{\prime}=\tau_v^i+\Delta \tau_v^i,\]$(C.9)

with τv0=0$\[\tau_{v}^{0}=0\]$ and Δτvi$\[\Delta \tau_{v}^{i}\]$ defined in Eq. (C.8).

References

  1. Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering by a Sphere (John Wiley & Sons, Ltd), 82 [Google Scholar]
  2. Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
  3. Cockburn, B. 2003, Z. Angew. Math. Mech., 83, 731 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cui, X., & Li, B. Q. 2005, JQSRT, 96, 383 [NASA ADS] [CrossRef] [Google Scholar]
  5. Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]
  6. Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
  7. Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: a multipurpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
  8. Dykema, P. G., Klein, R. I., & Castor, J. I. 1996, ApJ, 457, 892 [NASA ADS] [CrossRef] [Google Scholar]
  9. Frisken, S., & Perry, R. 2002, J. Graph. Tools, 7, 1 [CrossRef] [Google Scholar]
  10. Gravity Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Hesthaven, J. S., & Warburton, T. 2007, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (Springer Science & Business Media), 19 [Google Scholar]
  12. Ivezic, Z., Groenewegen, M. A. T., Men’shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] [CrossRef] [Google Scholar]
  13. 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]
  14. Kitzmann, D., Bolte, J., & Patzer, A. B. C. 2016, A&A, 595, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Kopriva, D. A., & Gassner, G. 2010, J. Sci. Comput., 44, 136 [CrossRef] [Google Scholar]
  16. Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
  18. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
  19. Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Machorro, E. 2007, J. Computat. Phys., 223, 67 [CrossRef] [Google Scholar]
  21. Mihalas, D., & Weibel-Mihalas, B. 1999, Foundations of Radiation Hydrodynamics, Dover Books on Physics (Dover Publications) [Google Scholar]
  22. Nair, R. D., Levy, M. N., & Lauritzen, P. H. 2011, Emerging Numerical Methods for Atmospheric Modeling (Berlin, Heidelberg: Springer Berlin Heidelberg), 251 [Google Scholar]
  23. Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
  24. Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Perdigon, J., Niccolini, G., & Faurobert, M. 2021, A&A, 653, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
  28. 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]
  29. Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63 [CrossRef] [Google Scholar]
  32. Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Wolf, S. 2003, Comput. Phys. Commun., 150, 99 [CrossRef] [Google Scholar]

1

Note that (r^,Θ^)$\[(\hat{r}, \hat{\Theta})\]$ is a plane of symmetry when axis-symmetry is assumed.

All Tables

Table 1

Summary of the main results.

All Figures

thumbnail Fig. 1

Illustration showing the coordinate system. We note that (r^,Θ^,Φ^)$\[(\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 Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates.

In the text
thumbnail Fig. 2

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

In the text
thumbnail Fig. 3

Example of the sparse structure of A$\[\mathcal{A}\]$ with Nr = 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 Aα,α,Aα,α+1,Aα,αNφ,Aα,α±NμNφ$\[\mathcal{A}^{\alpha, \alpha}, \mathcal{A}^{\alpha, \alpha+1}, \mathcal{A}^{\alpha, \alpha-N_{\varphi}}, \mathcal{A}^{\alpha, \alpha \pm N_{\mu} N_{\varphi}}\]$, and Aα,α±NΘNμNφ$\[\mathcal{A}^{\alpha, \alpha \pm N_{\Theta} N_{\mu} N_{\varphi}}\]$.

In the text
thumbnail Fig. 4

Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with τv0=1$\[\tau_{v_{0}}=1\]$ (blue curve) and τv0=102$\[\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
thumbnail Fig. 5

Temperature profiles for the axis-symmetric envelope with τv0=101$\[\tau_{v_{0}}=10^{-1}\]$ (blue curve) and τv0=102$\[\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 RADMC-3D. The lower panels show the relative differences between the two codes.

In the text
thumbnail Fig. 6

SED profiles for the axis-symmetric envelope with τv0=101$\[\tau_{v_{0}}=10^{-1}\]$ (blue curve) and τv0=102$\[\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 RADMC-3D. The lower panels show the relative differences between the two codes.

In the text
thumbnail 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 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 Ivenv$\[I_{v}^{\text {env }}\]$, inside each pixel. The solid white lines show the iso-contours of Ivenv$\[I_{v}^{\text {env }}\]$.

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

In the text
thumbnail Fig. 9

L2 norm of the relative differences in temperature between this study and RADMC-3D versus the number of elements in the computational grid (N = Nr × NΘ × Nμ × Nφ. The L2 norm was computed for the benchmark problem presented in Sec. 5.2.

In the text
thumbnail 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 (τ = 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
thumbnail 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
thumbnail Fig. B.1

Example of an element Di,j,k,l (dashed lines) using na = nc = nd = 3 and nb = 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 Fh$\[F_{h}^{*}\]$ at the interface, along the Θ coordinate.

In the text
thumbnail 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 κvext$\[\kappa_{v}^{\text {ext }}\]$ and Sv 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.