Issue 
A&A
Volume 647, March 2021



Article Number  A27  
Number of page(s)  12  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202039697  
Published online  01 March 2021 
Atomic line radiative transfer with MCFOST
I. Code description and benchmarking
^{1}
Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: benjamin.tessore@univgrenoblealpes.fr
^{2}
School of Physics and Astronomy, Monash University, VIC 3800, Australia
Received:
16
October
2020
Accepted:
19
December
2020
Aims. We present MCFOSTart, a new nonlocal thermodynamic equilibrium radiative transfer solver for multilevel atomic systems. The code is embedded in the 3D radiative transfer code MCFOST and is compatible with most of the MCFOST modules. The code is versatile and designed to model the close environment of stars in 3D.
Methods. The code solves for the statistical equilibrium and radiative transfer equations using the Multilevel Accelerated Lambda Iteration method. We tested MCFOSTart on spherically symmetric models of stellar photospheres as well as on a standard model of the solar atmosphere. We computed atomic level populations and outgoing fluxes and compared these values with the results of the TURBOspectrum and RH codes. Calculations including expansion and rotation of the atmosphere were also performed. We tested both the pure local thermodynamic equilibrium and the outofequilibrium problems.
Results. In all cases, the results from all codes agree within a few percent at all wavelengths and reach the subpercent level between RH and MCFOSTart. We still note a few marginal discrepancies between MCFOSTart and TURBOspectrum as a result of different treatments of background opacities at some critical wavelength ranges.
Key words: radiative transfer / methods: numerical
© B. Tessore et al. 2021
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.
1. Introduction
Radiation processes are critical in many aspects of the evolution of astrophysical objects such as stellar atmospheres, magnetospheres and discs. Electromagnetic radiation carries the history of the emitting object to the observer at which point the object is analysed. Through spectroscopy we can deduce fundamental parameters, abundances and unveil the complex atmospheric dynamics of many stellar objects, often via comparison with numerical models.
These numerical models require the solution of the radiative transport equation in complex and different plasma conditions. In many astrophysical applications matter and radiation are coupled : This is the socalled nonlocal thermodynamic equilibrium (nonLTE) line formation problem. At nonLTE, the atomic level populations are computed from a set of statistical equilibrium equations where radiation plays the central role. In turn, knowledge of emission and absorption processes, which depends on the level populations, is needed to solve for radiation. Therefore, the radiative transport equation must be solved simultaneously with statistical equilibrium equations. This makes the solution of the nonLTE problem challenging. Identifying the best solution to that problem depends mainly on the nature of the simulated plasma and many types of solutions exist (e.g. van Zadelhoff et al. 2002).
In this paper, we present MCFOSTart, a code for the solution of the nonLTE problem for multilevel atomic systems. The code is embedded in the 3D radiative transfer code MCFOST (Pinte et al. 2006, 2009) and can be applied to a wide range of astrophysical problems in different geometries. Our main motivation for developing such a code is to address the line formation in the close environment of young stellar objects (YSO). In particular, the evolution of classical T Tauri stars depends on the interaction between the young star and its accretion disc, on distances of a few stellar radii. This complex interaction leads to stellar winds and magnetospheric accretion and ejection processes. These processes have a strong impact on spectral lines and disentangling their typical radiative signatures is a challenging task.
Previous modelling of the environment of T Tauri stars, their magnetosphere (Hartmann et al. 1994; Muzerolle et al. 2001), stellar wind, and disc wind regions (Lima et al. 2010; Kurosawa et al. 2011) show the difficulty to account for all observational signatures of these stars, from their variability to the shape of emission lines. Even a detailed modelling of a specific T Tauri star, with the most updated magnetohydrodynamics (MHD) models, results in a marginal agreement between synthetic observables and observations (Alencar et al. 2012).
Most of magnetospheric emission lines synthesis used the Sobolev approximation (Sobolev 1957; Rybicki & Hummer 1978) to solve the radiative transfer equation. This approximation is applicable when large velocity gradients (i.e. large compared to the intrinsic width of atomic lines) are encountered in a plasma. It greatly simplifies the solution of the transfer equation. However, in the case of magnetospheres, the required conditions for the Sobolev approximation to hold are not always met, ultimately impacting the level populations. In our code, we use an accelerated Λiteration method (Olson et al. 1986) with a preconditioning of the statistical equilibrium equations (Rybicki & Hummer 1991) to solve the coupled radiative transfer and statistical equilibrium equations. Unlike the Sobolev approximation, our approach does not make any assumptions about the local conditions in the plasma.
We tested MCFOSTart on 1D spherically symmetric models including a range of thermodynamic quantities to capture a variety of physical conditions existing in stellar plasma. In Sect. 2 we give the details of the physical problem and how it is solved. In Sect. 3 we describe the implementation of our code in MCFOST, and in Sect. 4 we test the accuracy and precision of our code against benchmark cases. Finally, we give our conclusions in Sect. 5.
2. NonNLTE radiative transfer problem
In multilevel nonLTE problems, the main difficulty is to find a selfconsistent solution of the equations for the atomic level populations with a solution of the transport equation for radiation. In this section we describe the coupling between these equations.
2.1. Radiative transfer equation
The unpolarised radiative transfer equation along a ray of length ds in the direction n is written as
where, I(ν, n) is the specific intensity at a frequency ν in the direction n, χ is the total (line and continuum) opacity (true absorption and scattering σ), and η is the emissivity.
We use the formulation of Rybicki & Hummer (1992) to express the radiative transfer equation, and we define the radiative transfer coefficients U_{ℓℓ′} and V_{ℓℓ′} as
where A_{ℓℓ′} and B_{ℓℓ′} are the Einstein’s coefficients and ϕ(ν, n) and ψ(ν, n) the absorption and the emission profiles, respectively, for the transition i → j^{1}. We further assume complete frequency redistribution, implying an equivalence between the absorption coefficient appearing in V_{ℓℓ′} and the emission coefficient appearing in U_{ℓℓ′}, that is ψ(ν, n) = ϕ(ν, n).
For a continuum transition the radiative transfer coefficients are given by
where α(ν) is the photoionisation crosssection, T the gas temperature, and the population (number density) of level ℓ, evaluated LTE.
Finally, the emissivity η_{ij} and the opacity χ_{ij} for a transition between an upper level j and a lower level i are written as
where n_{i} and n_{j} are the lower and upper level populations, respectively. The total opacity and emissivity for all transitions of all atoms and other sources of opacity add up to the total opacity χ(ν, n) and emissivity η(ν, n),
where η_{c} and χ_{c} are the sources of continuous emissivity and opacity evaluated at LTE (the background opacity, see Sect. 3.4).
The solution of Eq. (1) along a ray of length ds is straight forward if χ, η, the populations n, and the properties of the atmosphere (temperature, velocity fields, etc.) are known. Defining the optical depth τ ≔ τ(ν, n; s) as dτ(ν, n) = − χ(ν, n; s) ds it reads:
where is the source function at optical depth τ, at the frequency ν in the direction n.
In general, the emissivity and opacity depend on the level populations which in turn depend on the intensity. Therefore, a selfconsistent solution of the radiative transfer equation with a solution of statistical equilibrium equations is required, the socalled nonLTE problem.
2.2. Λiteration
The solution of Eq. (1) can be recast as
where Ψ^{2} is a matrix operator whose elements depend on the level populations (for more details see e.g. Hubeny & Mihalas 2014, hereafter, HM14).
This equation gives the solution of the intensity if all elements of the Ψ operator, which are functions of the atomic populations, are known. However, in practice, building the full Ψ operator is never affordable and the emissivity appearing in Eq. (7) is substituted by an old value η^{†} evaluated from a known estimate of the level populations (i.e. old values obtained with a previous iteration). The new intensity obtained by the application of the operator on the old emissivity is then used to compute a new value of the level populations and of the emissivity, which are used to compute a new value of the intensity. This iterative process, between old and new values of the populations, is known as Λiteration and it is repeated until convergence. Although Λiteration is very efficient in an optically thin region, it suffers convergence problems in very optically thick regions. This drawback of classical Λiteration can be overcome by a slight modification of Eq. (7) (see Cannon 1973a,b; Scharmer 1981; Olson et al. 1986), that is recasting Eq. (7) in
with an approximate operator Ψ^{*}(ν, n) built from a subset of the original operator. Equation (8) improves the convergence of the classical Λiteration method in optically thick regions and has been named ALI for accelerated Λiteration (after the work of Hamann 1985 and Werner & Husfeld 1985).
In most applications, the diagonal of the full operator is used, although other studies recommended a tridiagonal approximate operator (see for instance Hennicker et al. 2018). Higher order approximate operators can improve the convergence of the radiative transfer problem at the price of computational time and memory storage. For 2D radiative transfer problems, Auer et al. (1994) suggested to use the diagonal part of the operator in combination with an algorithm to accelerate the convergence (see Sect. 2.7). In our code, we use the diagonal part of the full operator as a compromise between speed and accuracy.
2.3. Statistical equilibrium equations
The level populations of an atom in a stellar plasma is a function of the rates of transitions populating or depopulating a given level. These rates in turn are a function of level populations, collisions, and radiation.
The statistical equilibrium equations (SEE) for level ℓ is given by^{3}
where n_{ℓ} is the population of level ℓ, and C_{ℓ′ℓ}, and R_{ℓ′ℓ} the collisional and radiative rates for a transition between level ℓ′ and ℓ, respectively.
2.3.1. Radiative rates
Unlike collisional rates, radiative rates are a product of the radiative transfer code. These rates can be expressed in terms of the radiative transfer coefficients and are written as
where the integration is carried out over frequency ν and all solid angles Ω. For line transitions these rates are written as
where B_{ij}, B_{ji}, and A_{ji} are the Einstein’s coefficients, and the mean intensity integrated over the line given by
where ϕ(ν, n) is the line absorption profile, also appearing in Eq. (2). For continuum transitions they are as follows:
where 𝒥(ν) is the mean radiation field given by
2.3.2. Matrix form of the SEE
Starting from Eq. (9), passing the righthand side term on the lefthand side and noticing that with δ_{ℓℓ′} nonzero only for ℓ′ = ℓ, we can factorise out and write the SEE as
where the rate matrix Γ_{ℓ′ℓ} is
2.4. MALI method
Multilevel radiative transfer is a complex problem to solve because of the dependence of the radiation field on the populations, and of the populations on the radiation field. Using Eqs. (8) and (10) the rate matrix in Eq. (16) for a multilevel atom reads:
where I_{eff}(ν, n) = I^{†}(ν, n)−Ψ^{*}(ν, n)η^{†}(ν, n) is an effective radiation field built with known quantities (old values). In Eq. (17) we assumed that background opacities are constant within subsequent iterations.
The system of Eq. (15) is nonlinear in the new populations as it involves the product of the form ∝n_{ℓ} × n_{ℓ′} in Eq. (17). The multilevel accelerated lambda iteration method (hereafter, MALI) proposed by Rybicki & Hummer (1991), latter improved by Uitenbroek (2001), uses the operator splitting method and a full preconditioning to make Eqs. (15) and (17) linear in the new populations. The coupled equations of statistical equilibrium with the radiation field are finally written as follows:
The last term in Eq. (18), which disappears in the classic Λiteration (Rybicki & Hummer 1992; Uitenbroek 2001) is called the crosscoupling term. From Eq. (4) the term represents a summation over all opacities of all transitions with a plus sign if l < j or a minus sign if l > j, that is .
2.5. Coherent electron scattering
In some cases electron scattering is a dominant source of opacity and its effect on both the continuum and spectral lines has to be taken into account (Hillier 1991). However, the evaluation of the scattering emissivity involves the calculation of a scattering integral whose expression in the observer’s frame is not trivial (Rybicki & Hummer 1994). We evaluate the electron scattering emissivity by setting it to the mean intensity (coherent approximation), which is in turn determined iteratively through classical Λiterations, as suggested by Rybicki & Hummer (1992). This method is accurate enough for the main range of applications of the code, but becomes inaccurate for hot stellar wind where the electron density is very high and the ratio of scattering to true absorption is large.
2.6. Level dissolution and occupation probability
To do reliable radiative transfer simulations it is necessary to allow for line overlap and merging close to series limits. The MALI formulation of the nonLTE problem is already capable of dealing with overlapping lines, and we include the formalism of occupation probabilities of Hummer & Mihalas (1988) to allow for lines to merge at the series limits. Because of the presence of perturbers, each bound state i of an atom has a probability w(i), with respect to the same state of a similar isolated atom, to remain bound. Conversely, 1 − w(i) represents the probability of state i to be dissolved, in other words, to belong to the continuum states. Dissolved states contribute to a pseudocontinuum opacity beyond the series limits. Our implementation of the occupation probabilities formalism is similar to that of Hubeny et al. (1994), based on the work of Dappen et al. (1987), and summarised in Hillier & Miller (1998). We further assume, as in Hillier & Miller (1998), that if w(i) is the probability of level i to be undissolved and w(j) this probability for level j, then if j is undissolved, i is necessarily undissolved as it lies in a lower energy state. Therefore, all upward rates in the statistical equilibrium equations have to be multiplied by w(j)/w(i), the probability that the upper level j is undissolved given i is undissolved. In Sect. 4.1.1 the impact of level dissolution on stellar continua is shown.
2.7. Acceleration of the convergence
To speed up the MALI method we implemented the method of acceleration of Ng (1974) with general orders as defined in Auer et al. (1994) and Uitenbroek (2001). This method uses N_{sol} + 2 previous solutions and accumulated iteration after iteration to extrapolate the new value of the solution by minimising the residual between the previous iterations. The latter step is called Ng’s iteration. When using Ng’s accelerations, it is important to impose a delay before accumulating solutions, as extrapolating the new solution too early may result in negative or inconsistent values. Following Auer et al. (1994), we recommend doing some MALI iterations between each extrapolation to let the solution settle down. Therefore, after N_{start} MALI iterations (i.e. without acceleration), we start accumulating N_{sol} + 2 solutions up to the Ng’s iteration. We repeat this process every N_{pending} MALI iterations. In our tests, we use N_{sol} = 2, N_{start} = 5 and N_{pending} = 5. We note that setting N_{pending} to 0 mainly results in extra overhead due to the matrix inversion required in the minimisation procedure (see for instance Uitenbroek 2001).
The main drawback of Ng’s acceleration for multidimensional models is the memory required to store N_{sol} + 2 solutions for each atom treated in nonLTE. Following Auer et al. (1994), using a diagonal operator in addition to Ng’s iterations is a good compromise between speed (of convergence) and memory requirements in multidimensional models. Whereas in 1D geometry memory storage is usually not a limitation, we discuss the practicability of Ng’s acceleration scheme in multidimensional models in a future paper.
3. Implementation in MCFOST
The formulation of the nonLTE line transfer of Rybicki & Hummer (1992), although developed for 1D grids, can be easily adapted for multidimensional models. Auer et al. (1994) applied the MALI method to a 2D Cartesian grid, while Hauschildt & Baron (2014) applied it to a 3D spherical grid. Recently, De Ceuster et al. (2020) applied the method to an unstructured grid.
MCFOSTart, which stands for MCFOST atomic radiative transitions, is an ensemble of modules implemented in MCFOST (Pinte et al. 2006, 2009). A 3D radiative transfer code, MCFOST is written in modern Fortran, dedicated to the modelling of dust emission and molecular lines in the environment of young stars (e.g. protoplanetary discs, circumstellar envelopes). This code has been coupled with SPH codes (e.g. Price et al. 2018) and MHD codes (e.g. Riols et al. 2020). Currently MCFOST handles three types of grids: a cylindrical grid dedicated to the modelling of discs, a spherical grid, and an unstructured grid based on a Voronoi tessellation (see e.g. Camps et al. 2013). The solution of the line transfer equation is done by raytracing: several rays are propagated in specific directions along which Eq. (1) is solved for each wavelength simultaneously.
In our new modules, we propagated rays from each cell centre (spatial grid units in 2D and 3D geometries) in directions defined by the angular quadrature scheme (see below). Cells represent the smallest resolution elements of the code and all quantities are constant within them. If the model has a nonzero velocity field, the velocity is projected in the direction of propagation of the ray. When velocity fields are present, it is important to maintain a proper spatial discretisation to prevent nonphysical features due to large velocity gradients (Ibgui et al. 2013). The MCFOST code detects if the projected velocity between two grid points is too large and linearly interpolates the projected velocity accordingly.
To produce images, we solved the intensity out of each pixel by integrating Eq. (1) along rays centred on each pixel.
The solution of Eq. (1) along a ray, specified by a vector direction n, is written as
where I_{0}(ν, n) is the emerging specific intensity in a given direction; I_{b}(ν, n) is the intensity at the inner boundary of the model, typically this is the stellar radiation if a ray intersects the star; τ_{tot} is the total optical depth from the boundary to the observer; S_{k} the source function for cell k; and τ_{k} the optical depth at the inner (entrance) boundary of a cell. In Eq. (19), the sum represents the contribution of each cell, along the ray path, to the outgoing radiation.
The term dτ_{k} = τ_{k + 1} − τ_{k} = χ_{k} s is the optical depth accumulated by the ray as it travels a distance s inside a cell of opacity χ_{k}. The intensity outgoing in direction n from a specific cell, after travelling a distance s(k + 1) − s(k) in the cell is written as
The approximate operator Ψ^{*} appearing in Eqs. (8) and (18) is evaluated by solving the transfer equation with a source function, which is 1 at cell k and 0 elsewhere, that is S_{k} = δ(τ − τ_{k}). According to Eq. (20), this is given by
where the numerator is the diagonal of the Λ operator, that is
This satisfies the following conditions in the two extreme cases: in optically thin regions it is 0, which corresponds to classical Λiterations, and it approaches 1 in optically thick regions. Although this method allows for fast integration of the radiative transfer equation, it depends on the numerical resolution, compared to other methods that use linear or cubic interpolations between grid points.
The part of the code that handles the propagation within the grid is presented in Pinte et al. (2006). In the following, we present in some detail the solution of the nonLTE atomic line formation problem in the code.
3.1. Initial solution
In MCFOSTart, there are currently two choices for the initialisation of the level populations of the nonLTE problem: populations are evaluated at LTE and populations are obtained by solving the SEE with the radiation field set to zero. The choice of the initial solution is the critical point of the nonLTE problem. If the initial guess is too far from the solution, the convergence can be slowed down or even fail. In the eventuality that none of these choices lead to convergence, we implemented the collisionalradiative switching method of Hummer & Voels (1988). This method allows for a smooth transition between a collision dominated initial solution (i.e. at LTE) to a nonLTE initial solution. Alternatively, the populations from a previous calculation are given.
3.2. Line profiles
The absorption profile for each atomic boundbound transition can either have a Gaussian or a Voigt shape. The choice of the shape of the absorption profile depends on the line and the problem.
The width of the absorption profile, v_{D}, is given by the sum of the line thermal width and the microturbulence ξ as follows:
where k_{b} is the Boltzmann constant, T the temperature, and m the mass of the atomic species.
The damping parameter is the sum of the radiative damping and Van Der Waals and Stark broadenings. The radiative damping of a line transition j → i, Γ_{j}, is computed by summing the Einstein A_{ji} coefficients for all (allowed) transitions for which j is an upper level as follows:
In the latter expression, we neglected the effect of stimulated emission.
Pressure broadening takes into account van der Waals and Stark (linear and quadratic) broadenings computed in the Lindholm theory. We computed interaction constants using the Unsöld formula (Unsöld & Weidemann 1955) and evaluated the quadratic Stark broadening following Gray (2008, Eq. (11.33)). Linear Stark broadening was computed following Sutton (1978). The treatment of the Stark broadening, especially for hydrogen lines is approximate and is only accurate for the first lines of each series. In the future, we plan to add a more accurate treatment of the Stark broadening (see e.g. Stehlé & Hutcheon 1999).
3.3. Wavelength grid
The wavelength grid used for nonLTE calculations was built from individual atomic transitions. Firstly, we set up a grid of continuum transitions at a moderate resolution. This continuum grid was used to solve the continuum radiative transfer equation.
Secondly, lines were gathered in groups depending on their bounds (i.e. overlapping regions). Line bounds range from −V to +V, where V is a multiple of the line Doppler width (Eq. (23)). The line bounds represent the interaction zone of the line with the radiation field. Since in moving atmosphere a line feels the radiation from different regions, the line bounds include the maximum shift in frequency due to the velocity fields. However, the local (i.e. in the comoving frame) line profile is not affected by velocity fields.
Each group of lines was sampled with a constant resolution in kilometers per second. Once wavelength grids for each group were determined, they were merged and continuum points from the continuum grid were added outside line groups.
3.4. Background opacity
Opacity and emissivity of transitions were computed for each atom using Eqs. (2)–(4). The photoionisation crosssection of hydrogenlike ions was computed using Kramer’s formula with correction from quantum mechanics (HM14, Eq. (7.92)). For nonhydrogenlike ions, the photoionisation crosssection was directly read from the atomic file and interpolated on the wavelength grid of MCFOST. We took the photoionisation data from TOPbase^{4} (Cunto & Mendoza 1992).
We also included continuum transitions from other sources which add up in χ_{c}. These are Thomson scattering on free electrons, hydrogen freefree (HM14, Eq. (7.100)), and H^{−} freefree and boundfree transitions from John (1988).
3.5. NonLTE loop
The main module of our code deals with solving Eq. (15) simultaneously with Eq. (1) via the nonLTE loop. Figure 1 shows the main steps leading to the solution of the nonLTE problem with MCFOSTart. We stress that even if Eq. (15) is solved locally (i.e. within each cell) it is coupled with other cells through Eq. (1).
Fig. 1.
Flow chart of MCFOSTart. See Sect. 3.5 for more details. The green arrows indicate the main iterative loop. The double grey arrow indicates that the electron density loop is iterated within the main iterative loop and with the current estimates of the nonLTE populations. 
Firstly, the code reads abundances for the atomic species under consideration, the corresponding atomic data (e.g. energy levels E_{i} and the oscillator strength f_{ij}), and an atmospheric model (temperature, densities, and velocity fields). The electron density is evaluated at LTE if not present in the model (for more detail on electron density loop, see Sect. 3.6). Secondly, MCFOSTart initialises the level populations of each atom as described above.
Finally, we evaluate continuous background opacities. In models in which electron scattering is an important source of opacity, the continuum mean intensity and the corresponding Thomson scattering emissivity are computed with the ALI method (for more details, see chapter 13 of HM14). Once background opacities have been evaluated, the nonLTE loop starts. It computes a selfconsistent solution for level populations and radiative transfer equations. For each cell, the rate matrix Γ_{ℓ′ℓ} is built for a set of rays, for all wavelengths simultaneously. Once the rate matrix is known, a new value of the level populations is determined. Electron densities can be evaluated with the new values of the nonLTE populations. To avoid convergence issues, the evaluation of the electron density can be done every N iterations of the nonLTE loop.
The new populations are then compared to the previous populations using the following criterion of convergence:
where ϵ is a user defined convergence threshold. The threshold ϵ is usually chosen between 10^{−3}–10^{−4}. If the maximum relative changes for each cell and each atom is below the convergence threshold the nonLTE loop stops. Otherwise, a new iteration starts with the previously computed populations. Excitation and ionisation temperatures for each transition of each atom are also checked for convergence. In general, transition temperatures converge faster than level populations.
3.6. Electron density
Electron densities for a gas mixture of hydrogen and metals cannot be determined analytically. The populations of atomic levels rely on the knowledge of the electron densities (for instance, through the SahaBoltzmann equation at LTE). In turn, the electron density is a function of the population of atomic levels through their ionisation fraction. Therefore, electron densities have to be determined iteratively. Starting from an initial guess for the electron density, the populations of atomic levels are derived and the new ionisation fractions f_{is} for level i in the ionisation stage s for each electron donors are evaluated. Then, from the new ionisation fractions, electron densities are computed. This iterative procedure is repeated until the relative change in electron densities drops below 10^{−6}.
In our code we use a linearisation procedure similar to Eq. (17.88) of HM14. The initial guess for the electron density is given by the ionisation of hydrogen plus one metal (M) . Hydrogen ionisation yields
and the metal ionisation
where n_{H} is the total number of hydrogen atoms, ϕ_{x} the SahaBoltzmann factor for the ion x, and α_{M} is the abundance of the metal relative to hydrogen.
3.7. Collisional rates
The evaluation of the collisional rates appearing in the SEE (Eq. (9)) is cumbersome. It generally requires collision calculations with quantum mechanics for which lookup tables of collisional rates are available now (Barklem 2016). However, analytical and semianalytical recipes exist for faster evaluation of the collision matrix C_{ji} (which is proportional to the collisional rates within a factor). We consider collisional excitation and ionisation by electrons for hydrogen atoms using Johnson (1972). For collisional excitation of metals (and helium) by electrons we use van Regemorter’s fomula (van Regemorter 1962) and the impact parameter method (Seaton 1962). For ions, we follow Eq. (34) of Aller et al. (1982) and, for neutrals, we follow Seaton (1962). For collisional ionisation by electrons, we use Eq. (9.60) of HM14 and for collisions with hydrogen atoms, we use Drawin’s formula (Eq. (9.62) of HM14).
3.8. Quadratures
One of the cornerstones in solving the SEE is the evaluation of the integrals appearing in Eq. (10). In the following, we describe how we performed wavelength and angular integrations.
Angular quadrature for molecular lines transfer in MCFOST is done with the Monte Carlo method. For atomic lines, we used a set of fixed directions and weights uniformly distributed on the unit sphere such that integration over the sphere gives 4π. We used the directions and weights of type A quadratures from B. Carlson (whose calculations are described in Bruls et al. (1999)) and the directions and weights from Štěpán et al. (2020) for unpolarised radiation. These two types of quadrature points give similar results. By default we used type A quadrature of Carlson with 80 points sampling 4π steradians (Ibgui et al. 2013).
De Ceuster et al. (2020) used a similar angular quadrature method but with the quadrature points and weights determined by the HEALPix algorithm (Górski et al. 2005). For each cell of a model, starting from the centre of the cell (hereafter one point quadrature), the transfer equation is solved for all directions, while the radiative rates and the mean intensity are accumulated. Integrating the transfer equation from the centre of each cell may underestimate angular means and introduce some inaccuracy in the level populations because it only gives the mean intensity at that position. However, we do not know beforehand whether this inaccuracy is significant or not. Thus, we implemented two options to deal with this issue: first, several starting points for the angular quadrature scheme can be selected for each cell; and, second, Monte Carlo iterations can be done to properly sample each cell, in addition to the one point angular quadrature described above. In Sect. 4.3, we discuss these two options.
We replace the double integral appearing in Eq. (12) by a sum over wavelength points and directions as follows:
where ω_{ray} is the weight of the angular quadrature for a given direction (or ray) and ω_{λ} the weight of a trapezoidal wavelength quadrature. In the case of a line transition, ω_{λ} includes the normalisation of the line profile.
4. Benchmarking the code
The MCFOST code is intrinsically 3D, but in this paper we focus on 1D models that test the core of the method and allow us to validate our implementation with wellestablished 1D codes. While this is not the main application intended for our code, it also illustrates that MCFOST can now be used to generate a variety of stellar spectra. We will devote a forthcoming paper to the applications of our code to 2D and 3D geometries, with a particular emphasis on stellar magnetospheres. We used the codes TURBOspectrum (Plez 2012) and RH (Uitenbroek 2001) as references for testing.
The code MCFOST does not include a dedicated grid for 1D models, and we used the spherical grid instead. However, the spherical grid of MCFOST is logscaled in the radial direction and cannot efficiently map the nonuniform grid from the 1D models. Therefore, for our benchmarks, we directly read the grid from the input models as RH and TURBOspectrum do. In that case, a cell is defined as the distance between two consecutive grid points. Because Eq. (1) is solved in MCFOST by summing up the intensity of each cell (i.e. spatial grid unit) weighted by the optical depth (see Eq. (19)), the solution is more dependent on the grid resolution (i.e. on the size of cells). On the contrary, RH and TURBOspectrum formal solvers typically interpolate the source function in the initial grid, resulting in a solution that is highly accurate. To mitigate the impact of a lack of resolution on the accuracy of the solution, we decreased the size of the cells by adding more points between each point of the initial models.
The TURBOspectrum code is used to model spectra from cool evolved stars (Alvarez & Plez 1998; de Laverny et al. 2012) in synergy with the model atmosphere code MARCS (Gustafsson et al. 2008). This code works only at LTE but includes a wide library of molecules, the occupation probabilities formalism, and a proper treatment of hydrogen line broadening. To match the capabilities of our code, we slightly modified TURBOspectrum and removed molecular opacities, although, at the coolest temperature, neutral hydrogen density is dependent on the formation of molecules, such as H_{2}.
The radiative transfer code RH is mainly used to model line formation in nonLTE conditions using the MALI method. Recent versions of RH have been used to model spectral lines taking into account the Zeeman effect and partial frequency redistribution in 3D cubes of the solar photosphere (de la Cruz Rodríguez et al. 2019). Recently, Criscuoli et al. (2020) benchmarked the RH code and showed its ability to model the solar irradiance with a good agreement with observations. The RH and MCFOST codes use similar treatment of background opacities and both employ the MALI method to solve for level populations. We compared our treatment of background opacities and level dissolution at LTE with TURBOspectrum. We used RH for benchmarking the MALI method on the nonLTE formation of hydrogen lines. We used the (1D) spherically symmetric grid of TURBOspectrum and RH for all tests.
4.1. Stellar photospheres
To test our opacity modules and background opacity calculations, we used LTE models of stellar photospheres from the MARCS database^{5} (Gustafsson et al. 2008). The summary of each model is given in Table 1. The first model, sun, is the standard solar photosphere model from the MARCS code. The last two models, s8000 and s3500, are representative of the photospheres of hot and cool giants, respectively. All models have standard chemical composition and are computed assuming spherical symmetry, except for the solar model, which was computed using the planeparallel approximation.
Stellar atmosphere models used for the benchmarking of our code.
4.1.1. Continuum radiation
We compared the stellar continuum computed by the three codes for the stellar photosphere models (see Table 1). Figure 2 shows the continuum flux and the discrepancy function δF/F for each model.
Fig. 2.
Continuum flux as a function of wavelengths for the photospheric models without (left panels) and with (right panels) level dissolution. In each panel, the upper window shows the flux from RH, TURBOspectrum, and MCFOST in cyan (thin line), orange (dashed), and black (thick line), respectively. The lower window shows the difference between RH and MCFOST (thick dark grey lines), and between TURBOspectrum and MCFOST (dashed light grey lines). As RH does not include level dissolution, it is not shown in the right figures. In Fig. 2c, the continuum flux computed by RH stops at about 800 nm, leading to incorrect evaluation of the error at this edge point. (a) Solar photosphere without level dissolution. (b) Solar photosphere with level dissolution. (c) s8000 model without level dissolution. (d) s8000 model with level dissolution. (e) s3500 model without level dissolution. (f) s3500 model with level dissolution. 
The agreement between RH and MCFOST is excellent, being at or below 1% at most wavelengths. For models sun (Fig. 2a) and s3500 (Fig. 2e), the largest discrepancies are due to interpolation errors at the shortest wavelengths. For the model s8000 (Fig. 2c) the agreement between RH and MCFOST is also at the percent level, except at the Balmer jump (about 350 nm), where differences are larger (again owing to the steep slope of the emission), but remain below 30%. Discrepancies between MCFOST and TURBOspectrum are of the same order as those between MCFOST and RH. The largest discrepancy for the model s3500, of the order of 10%, is due to the H^{−} continuum. In TURBOspectrum, the H^{−} density is selfconsistently determined by solving chemical equilibrium equations including hydrogen in its atomic and molecular forms, leading to a slightly different value than in RH and MCFOST. For the model s8000, the largest discrepancy occurs at the Balmer jump and is of the order of 40%. As the largest discrepancy with RH occurs at the same location we might not exclude numerical resolution problems. The agreement between RH and TURBOspectrum is of the same order as MCFOST and TURBOspectrum. Furthermore, we note that for this model, continuum scattering becomes nonnegligible and is automatically included in TURBOspectrum^{6}, but not in MCFOST, thereby impacting the shape of the continuum.
The agreement between MCFOST and TURBOspectrum remains excellent when level dissolution is taken into account (Fig. 2 right panels). The impact of level dissolution is clearly visible at the Balmer jump where it smoothes the abrupt discontinuity towards redder wavelengths. The shallower spectrum results in an improved agreement between TURBOspectram and MCFOST for model s8000 around the Balmer jump.
4.1.2. Velocity fields
Although the impact of velocity fields on line formation will be addressed in a forthcoming paper, in this section we show an example of a solar photosphere with a radial velocity V_{r} given by
where v_{0} is the velocity at the inner point (bottom of the photosphere), V_{inf} a velocity such that max(V_{r}) = 500 km s^{−1}, and β the exponent of the velocity law. In our calculations, v_{0} = 9 km s^{−1} and β = 0.5.
We compared the result of our code with RH since TURBOspectrum does not deal with macroscopic velocity fields. For comparison, we implemented a solver for the velocity field in RH because the spherically symmetric version does not handle velocity fields. However, unlike MCFOST, the velocity fields between two points are not interpolated.
Figure 3 shows the difference between MCFOST and RH for the solar photosphere when a velocity field is included. While the continuum radiation is similar to that of Fig. 2a, lines are impacted significantly by the large velocity gradient, as expected. The differences between the two codes are negligible, although the lack of resolution in RH is noticeable.
Fig. 3.
Flux for a solar photospheric model with velocity gradient given by Eq. (29). The flux from MCFOST is shown in black and from RH in cyan. A zoomin around Lyβ and Hα is also shown. 
We also tested the treatment of the same model with no radial velocity but a constant rotational velocity of 300 km s^{−1} instead. Figure 4 shows the Hα line flux as a function of inclination. As expected, the impact of rotation on the line shape increases with inclination, from no effect when the star is seen poleon to maximum broadening when the star is equatoron.
Fig. 4.
Hα flux variation with inclination for a solar photospheric model with rotation. The flux at inclinations of 0°, 45°, and 90° are shown in black, green, and red, respectively. An inclination of 0° means that the star is seen pole on. 
4.2. Solar atmosphere model
Although a LTE model of the solar photosphere is adequate to model most solar lines, it fails to reproduce lines formed beyond the photosphere. The upper solar atmosphere is a high temperature and low density region above the photosphere where nonLTE effects dominate the line formation. In the following we focus on the Lyman α (Lyα) and Hα lines as they form in the upper atmosphere (e.g. Fig. 1 of Vernazza et al. 1981) and are thus impacted by nonLTE effects.
We used a 6level hydrogen atom, with five bound levels and one continuum level, which is the ground state of H II. Energy levels and transition frequencies are from Johnson (1972). Boundfree crosssections are computed using Kramer’s formula and collisional rates are evaluated following Johnson (1972). For the Lyman α and Hα lines, the absorption profiles are given by a Voigt function. For the other lines, we used a Gaussian profile. We used model C of Fontenla et al. (1999) as a standard model of the solar atmosphere including the chromosphere and transition region.
Overall, the agreement between RH and MCFOST is excellent. Figure 5 shows the ratio b between nonLTE and LTE populations in model C. In the photosphere (below the temperature minimum), the populations are mostly at LTE with a b ratio of a few. In the deep photosphere, b is close to 1. From the temperature minimum to the transition region (the chromosphere), the populations start to depart from their LTE values. In the lower corona, the populations are very far from LTE, with a b factor greater than 10^{6} for the ground state of HI (level 1 in Fig. 5). The formation regions of the two lines are consistent with the earlier work of Vernazza et al. (1981). Contribution functions (i.e. regions of formation) of the Lyα and Hα lines are shown in Fig. 6.
Fig. 5.
Ratio of nonLTE to LTE populations (departure coefficient b) of a 6level hydrogen atom as a function of density (ρ). RH populations are indicated with cyan dots. MCFOST populations are shown with full colour lines, where darker (lighter) red shades mean lower (higher) energy levels. For clarity, the first three atomic levels and the continuum level (ground state of HII) are labelled explicitly in the insert. The different layers of the solar atmosphere are shown with vertical dashed lines. Typical heights above the solar radius are also indicated. 
Fig. 6.
Line contribution function, (response function to height variation), at disc centre. The contribution function is drawn as a function of the density in the atmosphere and distance from line centre (in nm). Lighter areas show where the contribution function peaks. We indicate the different layers of the solar atmosphere in each panel. (a) Hα line. (b) Lyα line. 
Figure 7 shows the solar flux computed in model C with RH and MCFOST, at nonLTE (a) and LTE (b). The impact of nonLTE effects on spectral lines is clearly visible. At LTE all lines are in emission, whereas at nonLTE most of the lines are seen in absorption. For comparison, discintegrated observations of the Sun show that most of the hydrogen lines are in absorption, except for strong chromospheric lines (e.g. Lyα).
Fig. 7.
Solar flux for the model C. RH is shown with the thin cyan and MCFOST with the thick black line. Below each panel the discrepancy between the two codes is shown. (a) NonLTE flux. (b) LTE flux. 
At all wavelengths, the discrepancy is below 1%, except in the core of some lines where it reaches 5%. In the figures, interpolation errors create glitches in the discrepancy functions. However, real differences exist in the core of strong lines, such as the Lyα line, where the mismatch can reach 20%–25%. Figure 8 shows the nonLTE flux of the Lyα line (top panel) and the flux of the Hα line (bottom panel). Although these differences are not critical, using one of the two options discussed in Sect. 3.8 to improve the angular quadrature leads to a better agreement between the Lyman α line computed by RH and MCFOST.
Fig. 8.
Zoomin on solar lines flux. The discrepancy between RH and MCFOST is shown below each panel. (a) Lyα line. (b) Hα line. 
4.3. Discussion on the angular quadrature
In Sect. 3.8, we presented the method we used to perform angular quadratures. These quadratures are important since they are used to evaluate the radiative rates, which in turn determine the level populations. As stated earlier, angular quadratures are performed for a set of fixed rays (i.e. directions) starting from each cell centre. In the following, we test option (i) of Sect. 3.8, that is performing the angular quadrature not only from cell centres but also at N random positions in the cells. As in Sect. 4.2, we compute the nonLTE Lyα line flux in model C. Further, we use the original grid of the 1D model, that is without additional points. We randomly chose 99 additional positions in each cell. We then kept these positions fixed during the nonLTE loop to avoid introducing Monte Carlo noise in the solution. For each of these 100 positions, the angular quadrature was performed for 80 rays resulting in a total of 8000 rays. Thus, for each cell, we obtained 100 values of the mean intensity and radiative rates, each at a different location inside the cell. The mean intensity and the radiative rates in each cell is just the arithmetic mean of these values. Finally, we interpolated the level populations on the same grid used in Sect. 4.2 (i.e. with a higher resolution). Because now we have a better estimation of the radiative rates for each cell, the error in the numerical integration of Eq. (1) dominates the discrepancy between RH and MCFOST.
In Fig. 9 we show the Lyα flux computed with this method compared to the profile shown in in Fig. 8a. The agreement between RH and MCFOST is substantially improved, with a maximum discrepancy below 5%, five times better than with only one position for the angular quadrature. Still we note that, even if the core is better reproduced, the wings are not, although the maximum error is small.
Fig. 9.
Evolution of Lyα flux with angular quadrature methods. The black and cyan lines are the flux computed by MCFOST and RH, respectively, from Fig. 8a. The dotted red line and the pink triangles are the flux computed using option (i) and (ii) of Sect. 3.8, respectively. The bottom panel shows the discrepancy between RH and MCFOST, from Fig. 8a (dark grey line), option (i) (light grey dotted line), and option (ii) (grey triangles). 
Similar results are obtained by performing a purely Monte Carlo step after the one point quadrature, that is we randomly selected 1000 positions and directions (Sect. 3.8, option (ii)). The main difference between Sect. 3.8 option (i) and Sect. 3.8 option (ii) lies in the execution time of the code. The latter approach is eight times faster than the former for an equivalent accuracy. However, the populations computed with Sect. 3.8 option (ii) have Monte Carlo noise.
5. Conclusions
In this paper, we tested our code against spherically symmetric models of stellar atmospheres, representing a range of physical plasma conditions. We tested background opacity calculations with lines merging close to series limits using the occupation probabilities formalism. The agreement between MCFOSTart and both TURBOspectrum and RH is excellent, of the order of 1% or below at most wavelengths. Still, a few discrepancies are present, arising from the different treatment of opacities at a few critical wavelengths (e.g. at the Balmer jump or at the H^{−} minimum). The code has also been successfully tested with vertical and rotational velocities. The proper behaviour of all spectral lines is recovered in models with nonzero velocity fields.
We performed nonLTE calculations of the solar Lyman α and Hα lines formation. These calculations were compared with calculations performed with the RH code. Line formation regions and nonLTE effects are reproduced well by our code compared to earlier studies and observations. While the discrepancy in the flux between RH and MCFOST can reach 25% in the core of Lyα, the discrepancy drops below 1% around the Hα line and for most other wavelengths in the spectrum. Furthermore, we showed that improving the accuracy of the angular quadrature scheme decreases this discrepancy by a factor of five.
Although we only tested the code against 1D models, the opacity modules and the nonLTE loop are geometry independent. These modules are therefore applicable to multidimensional models. The MCFOST code is 3D in essence and the solution of the radiative transfer equation is fully performed in 3D. Further applications of this new code to 2D and 3D geometries will be discussed in a forthcoming paper with a particular focus on stellar magnetospheres.
It is common in radiative transfer problems to see the Λ operator when it comes to “Λiteration”. The relation between the Ψ and Λ operators is given in Rybicki & Hummer (1991) as Ψ = Λ/χ.
Acknowledgments
We are thankful to Dr. Bertrand Plez for providing and maintaining a publicly available version of TURBOspectrum. B Tessore is grateful to Dr. Han Uitenbroek for providing a version of RH. We thank the anonymous referee for his comments, which helped to improve the manuscript. B. Tessore is also thankful to the french minister of Europe and foreign affairs and the minister of superior education, research and innovation (MEAE and MESRI) for research funding through FASIC partnership. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 742095; SPIDI: StarPlanetsInner DiskInteractions, http://www.spidieu.org). C. Pinte acknowledges funding from the Australian Research Council via FT170100040 and DP180104235.
References
 Alencar, S. H. P., Bouvier, J., Walter, F. M., et al. 2012, A&A, 541, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aller, L. H., Appenzeller, I., Baschek, B., et al. 1982, LandoltBörnstein: Numerical Data and Functional Relationshipsin Science and Technology– New Series Gruppe/Group 6 Astronomy andAstrophysics, Volume 2 Schaifers/Voigt: Astronomy and Astrophysics /Astronomie und Astrophysik Stars and Star Clusters / Sterne undSternhaufen, (Berlin Heidelberg New York: SpringerVerlag) [Google Scholar]
 Alvarez, R., & Plez, B. 1998, A&A, 330, 1109 [NASA ADS] [Google Scholar]
 Auer, L., Bendicho, P. F., & Trujillo Bueno, J. 1994, A&A, 292, 599 [Google Scholar]
 Barklem, P. S. 2016, Phys. Rev. A, 93, 042705 [NASA ADS] [CrossRef] [Google Scholar]
 Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 [NASA ADS] [Google Scholar]
 Camps, P., Baes, M., & Saftly, W. 2013, A&A, 560, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cannon, C. J. 1973a, J. Quant. Spec. Radiat. Transf., 13, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Cannon, C. J. 1973b, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Criscuoli, S., Rempel, M., Haberreiter, M., et al. 2020, Sol. Phys., 295, 50 [CrossRef] [Google Scholar]
 Cunto, W., & Mendoza, C. 1992, Rev. Mex. Astron. Astrofis., 23, 107 [NASA ADS] [Google Scholar]
 Dappen, W., Anderson, L., & Mihalas, D. 1987, ApJ, 319, 195 [NASA ADS] [CrossRef] [Google Scholar]
 De Ceuster, F., Homan, W., Yates, J., et al. 2020, MNRAS, 492, 1812 [NASA ADS] [CrossRef] [Google Scholar]
 de la Cruz Rodríguez, J., Leenaarts, J., Danilovic, S., & Uitenbroek, H. 2019, A&A, 623, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Laverny, P., RecioBlanco, A., Worley, C. C., & Plez, B. 2012, A&A, 544, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fontenla, J., White, O. R., Fox, P. A., Avrett, E. H., & Kurucz, R. L. 1999, ApJ, 518, 480 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W. R. 1985, A&A, 148, 364 [Google Scholar]
 Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2014, A&A, 566, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J. 1991, A&A, 247, 455 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151 [NASA ADS] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres [Google Scholar]
 Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Hummer, D. G., & Voels, S. A. 1988, A&A, 192, 279 [NASA ADS] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 John, T. L. 1988, A&A, 193, 189 [NASA ADS] [Google Scholar]
 Johnson, L. C. 1972, ApJ, 174, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Kurosawa, R., Romanova, M. M., & Harries, T. J. 2011, MNRAS, 416, 2623 [NASA ADS] [CrossRef] [Google Scholar]
 Lima, G. H. R. A., Alencar, S. H. P., Calvet, N., Hartmann, L., & Muzerolle, J. 2010, A&A, 522, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muzerolle, J., Calvet, N., & Hartmann, L. 2001, ApJ, 550, 944 [NASA ADS] [CrossRef] [Google Scholar]
 Ng, K. C. 1974, J. Chem. Phys., 61, 2680 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797 [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]
 Plez, B. 2012, Turbospectrum: Code for Spectral Synthesis [Google Scholar]
 Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270 [Google Scholar]
 Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95 [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1978, ApJ, 219, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1992, A&A, 262, 209 [NASA ADS] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1994, A&A, 290, 553 [NASA ADS] [Google Scholar]
 Scharmer, G. B. 1981, ApJ, 249, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J. 1962, Proc. Phys. Soc., 79, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1957, Sov. Astron., 1, 678 [NASA ADS] [Google Scholar]
 Stehlé, C., & Hutcheon, R. 1999, A&A, 140, 93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sutton, K. 1978, J. Quant. Spec. Rad. Transf., 20, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Unsöld, A., & Weidemann, V. 1955, Vist. Astron., 1, 249 [CrossRef] [Google Scholar]
 van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
 van Zadelhoff, G. J., Dullemond, C. P., van der Tak, F. F. S., et al. 2002, A&A, 395, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [Google Scholar]
 Štěpán, J., Jaume Bestard, J., & Trujillo Bueno, J. 2020, A&A, 636, A24 [CrossRef] [EDP Sciences] [Google Scholar]
 Werner, K., & Husfeld, D. 1985, A&A, 148, 417 [NASA ADS] [Google Scholar]
All Tables
All Figures
Fig. 1.
Flow chart of MCFOSTart. See Sect. 3.5 for more details. The green arrows indicate the main iterative loop. The double grey arrow indicates that the electron density loop is iterated within the main iterative loop and with the current estimates of the nonLTE populations. 

In the text 
Fig. 2.
Continuum flux as a function of wavelengths for the photospheric models without (left panels) and with (right panels) level dissolution. In each panel, the upper window shows the flux from RH, TURBOspectrum, and MCFOST in cyan (thin line), orange (dashed), and black (thick line), respectively. The lower window shows the difference between RH and MCFOST (thick dark grey lines), and between TURBOspectrum and MCFOST (dashed light grey lines). As RH does not include level dissolution, it is not shown in the right figures. In Fig. 2c, the continuum flux computed by RH stops at about 800 nm, leading to incorrect evaluation of the error at this edge point. (a) Solar photosphere without level dissolution. (b) Solar photosphere with level dissolution. (c) s8000 model without level dissolution. (d) s8000 model with level dissolution. (e) s3500 model without level dissolution. (f) s3500 model with level dissolution. 

In the text 
Fig. 3.
Flux for a solar photospheric model with velocity gradient given by Eq. (29). The flux from MCFOST is shown in black and from RH in cyan. A zoomin around Lyβ and Hα is also shown. 

In the text 
Fig. 4.
Hα flux variation with inclination for a solar photospheric model with rotation. The flux at inclinations of 0°, 45°, and 90° are shown in black, green, and red, respectively. An inclination of 0° means that the star is seen pole on. 

In the text 
Fig. 5.
Ratio of nonLTE to LTE populations (departure coefficient b) of a 6level hydrogen atom as a function of density (ρ). RH populations are indicated with cyan dots. MCFOST populations are shown with full colour lines, where darker (lighter) red shades mean lower (higher) energy levels. For clarity, the first three atomic levels and the continuum level (ground state of HII) are labelled explicitly in the insert. The different layers of the solar atmosphere are shown with vertical dashed lines. Typical heights above the solar radius are also indicated. 

In the text 
Fig. 6.
Line contribution function, (response function to height variation), at disc centre. The contribution function is drawn as a function of the density in the atmosphere and distance from line centre (in nm). Lighter areas show where the contribution function peaks. We indicate the different layers of the solar atmosphere in each panel. (a) Hα line. (b) Lyα line. 

In the text 
Fig. 7.
Solar flux for the model C. RH is shown with the thin cyan and MCFOST with the thick black line. Below each panel the discrepancy between the two codes is shown. (a) NonLTE flux. (b) LTE flux. 

In the text 
Fig. 8.
Zoomin on solar lines flux. The discrepancy between RH and MCFOST is shown below each panel. (a) Lyα line. (b) Hα line. 

In the text 
Fig. 9.
Evolution of Lyα flux with angular quadrature methods. The black and cyan lines are the flux computed by MCFOST and RH, respectively, from Fig. 8a. The dotted red line and the pink triangles are the flux computed using option (i) and (ii) of Sect. 3.8, respectively. The bottom panel shows the discrepancy between RH and MCFOST, from Fig. 8a (dark grey line), option (i) (light grey dotted line), and option (ii) (grey triangles). 

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.