Atomic line radiative transfer with MCFOST I. Code description and benchmarking
AAstronomy & Astrophysics manuscript no. bt_spidi1 © ESO 2021January 14, 2021
Atomic line radiative transfer with MCFOST
I. Code description and benchmarking
B. T essore , C. P inte , , J. B ouvier , and F. M´ enard Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France School of Physics and Astronomy, Monash University, VIC 3800, AustraliaReceived 16 / / / / ABSTRACT
Aims.
We present MCFOST-art, a new non-local 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 isversatile 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 LambdaIteration (MALI) method. We tested MCFOST-art on spherically symmetric models of stellar photospheres as well as on a standardmodel of the solar atmosphere. We computed atomic level populations and outgoing fluxes and compared these values with the resultsof the TURBOspectrum and RH codes. Calculations including expansion and rotation of the atmosphere were also performed. Wetested both the pure local thermodynamic equilibrium and the out-of-equilibrium problems.
Results.
In all cases, the results from all codes agree within a few percent at all wavelengths and reach the sub-percent level betweenRH and MCFOST-art. We still note a few marginal discrepancies between MCFOST-art and TURBOspectrum as a result of di ff erenttreatments of background opacities at some critical wavelength ranges. Key words.
Radiative transfer - Methods: numerical
1. Introduction
Radiation processes are critical in many aspects of the evolutionof astrophysical objects such as stellar atmospheres, magneto-spheres and discs. Electromagnetic radiation carries the historyof the emitting object to the observer at which point the objectis analysed. Through spectroscopy we can deduce fundamentalparameters, abundances and unveil the complex atmospheric dy-namics of many stellar objects, often via comparison with nu-merical models.These numerical models require the solution of the radiativetransport equation in complex and di ff erent plasma conditions.In many astrophysical applications matter and radiation are cou-pled : This is the so-called non-local thermodynamic equilibrium(non-LTE) line formation problem. At non-LTE, the atomic levelpopulations are computed from a set of statistical equilibriumequations where radiation plays the central role. In turn, knowl-edge of emission and absorption processes, which depends onthe level populations, is needed to solve for radiation. Therefore,the radiative transport equation must be solved simultaneouslywith statistical equilibrium equations. This makes the solution ofthe non-LTE problem challenging. Identifying the best solutionto that problem depends mainly on the nature of the simulatedplasma and many types of solutions exist (e.g. van Zadelho ff et al. 2002).In this paper, we present MCFOST-art, a code for the solu-tion of the non-LTE problem for multilevel atomic systems. Thecode is embedded in the 3D radiative transfer code MCFOST(Pinte et al. 2006, 2009) and can be applied to a wide range ofastrophysical problems in di ff erent geometries. Our main moti-vation for developing such a code is to address the line forma- tion in the close environment of young stellar objects (YSO).In particular, the evolution of classical T Tauri stars depends onthe interaction between the young star and its accretion disc, ondistances of a few stellar radii. This complex interaction leadsto stellar winds and magnetospheric accretion and ejection pro-cesses. These processes have a strong impact on spectral linesand disentangling their typical radiative signatures is a challeng-ing task.Previous modelling of the environment of T Tauri stars, theirmagnetosphere (Hartmann et al. 1994; Muzerolle et al. 2001),stellar wind, and disc wind regions (Lima et al. 2010; Kuro-sawa et al. 2011) show the di ffi culty to account for all obser-vational signatures of these stars, from their variability to theshape of emission lines. Even a detailed modelling of a specificT Tauri star, with the most updated magneto-hydrodynamics(MHD) models, results in a marginal agreement between syn-thetic observables and observations (Alencar et al. 2012).Most of magnetospheric emission lines synthesis used theSobolev approximation (Sobolev 1957; Rybicki & Hummer1978) to solve the radiative transfer equation. This approxima-tion is applicable when large velocity gradients (i.e. large com-pared to the intrinsic width of atomic lines) are encountered in aplasma. It greatly simplifies the solution of the transfer equation.However, in the case of magnetospheres, the required conditionsfor the Sobolev approximation to hold are not always met, ulti-mately impacting the level populations. In our code, we use anaccelerated Λ -iteration method (Olson et al. 1986) with a pre-conditioning of the statistical equilibrium equations (Rybicki &Hummer 1991) to solve the coupled radiative transfer and statis-tical equilibrium equations. Unlike the Sobolev approximation, Article number, page 1 of 13 a r X i v : . [ a s t r o - ph . S R ] J a n & A proofs: manuscript no. bt_spidi1 our approach does not make any assumptions about the localconditions in the plasma.We tested MCFOST-art on 1D spherically symmetric mod-els including a range of thermodynamic quantities to capture avariety of physical conditions existing in stellar plasma. In §2we give the details of the physical problem and how it is solved.In §3 we describe the implementation of our code in MCFOST,and in §4 we test the accuracy and precision of our code againstbenchmark cases. Finally, we give our conclusions in §5.
2. Non-NLTE radiative transfer problem
In multilevel non-LTE problems, the main di ffi culty is to find aself-consistent solution of the equations for the atomic level pop-ulations with a solution of the transport equation for radiation. Inthis section we describe the coupling between these equations. The unpolarised radiative transfer equation along a ray of length ds in the direction n is written as dI ( ν, n ) ds = − χ ( ν, n ) I ( ν, n ) + η ( ν, n ) (1)where, I ( ν, n ) is the specific intensity at a frequency ν in the di-rection n , χ is the total (line and continuum) opacity (true ab-sorption and scattering σ ), and η is the emissivity.We use the formulation of Rybicki & Hummer (1992) to ex-press the radiative transfer equation, and we define the radiativetransfer coe ffi cients U (cid:96)(cid:96) (cid:48) and V (cid:96)(cid:96) (cid:48) as U (cid:96)(cid:96) (cid:48) ( ν, n ) = h ν π A (cid:96)(cid:96) (cid:48) ψ ( ν, n ) if (cid:96) > (cid:96) (cid:48) V (cid:96)(cid:96) (cid:48) ( ν, n ) = h ν π B (cid:96)(cid:96) (cid:48) φ ( ν, n ) , (2)where A (cid:96)(cid:96) (cid:48) and B (cid:96)(cid:96) (cid:48) are the Einstein’s coe ffi cients and φ ( ν, n ) and ψ ( ν, n ) the absorption and the emission profiles, respectively, forthe transition i → j . We further assume complete frequencyredistribution, implying an equivalence between the absorptioncoe ffi cient appearing in V (cid:96)(cid:96) (cid:48) and the emission coe ffi cient appear-ing in U (cid:96)(cid:96) (cid:48) , that is ψ ( ν, n ) = φ ( ν, n ).For a continuum transition the radiative transfer coe ffi cientsare given by U (cid:96)(cid:96) (cid:48) = α (cid:96)(cid:96) (cid:48) ( ν )( n ∗ (cid:96) (cid:48) n ∗ (cid:96) ) 2 h ν c e − h ν/ kT if (cid:96) > (cid:96) (cid:48) V (cid:96)(cid:96) (cid:48) = α (cid:96)(cid:96) (cid:48) ( ν )( n ∗ (cid:96) (cid:48) n ∗ (cid:96) ) e − h ν/ kT if (cid:96) > (cid:96) (cid:48) V (cid:96)(cid:96) (cid:48) = α (cid:96)(cid:96) (cid:48) if (cid:96) < (cid:96) (cid:48) , (3)where α ( ν ) is the photoionisation cross-section, T the gas tem-perature, and n ∗ (cid:96) the population (number density) of level (cid:96) , eval-uated at local thermodynamic equilibrium (LTE). We adopt the following convention: j and i refer to the upper andlower levels of a transition, respectively. For atomic levels ordered byincreasing energies, E (cid:96) , this implies E j > E i . In other cases, when theordering of levels is unimportant, we use the indexes (cid:96), (cid:96) (cid:48) , (cid:96) (cid:48)(cid:48) ... Finally, the emissivity η i j and the opacity χ i j for a transitionbetween an upper level j and a lower level i are written as χ i j ( ν, n ) = n i V i j ( ν, n ) − n j V i j ( ν, n ) η i j ( ν, n ) = n j U ji ( ν, n ) , (4)where n i and n j are the lower and upper level populations, re-spectively. The total opacity and emissivity for all transitions ofall atoms and other sources of opacity add up to the total opacity χ ( ν, n ) and emissivity η ( ν, n ), χ ( ν, n ) = (cid:88) j , i < j n i V i j ( ν, n ) − n j V i j ( ν, n ) + χ c η ( ν, n ) = (cid:88) j , i < j n j U ji ( ν, n ) + η c , (5)where η c and χ c are the sources of continuous emissivity andopacity evaluated at LTE (the background opacity, see §3.4).The solution of Eq. 1 along a ray of length ds is straight for-ward if χ , η , the populations n , and the properties of the atmo-sphere (temperature, velocity fields, etc...) are known. Definingthe optical depth τ (cid:66) τ ( ν, n ; s ) as d τ ( ν, n ) = − χ ( ν, n ; s ) ds itreads: I ( ν, n ; τ ) = I ( ν, n ; τ ) e − ( τ − τ ) + τ (cid:90) τ S ( ν, n ; τ ) e − ( τ − τ ) d τ , (6)where S ( ν, n ; τ ) = η ( ν, n ; τ ) χ ( ν, n ; τ ) is the source function at opticaldepth τ , at the frequency ν in the direction n .In general, the emissivity and opacity depend on the levelpopulations which in turn depend on the intensity. Therefore, aself-consistent solution of the radiative transfer equation with asolution of statistical equilibrium equations is required, the so-called non-LTE problem. Λ − iteration The solution of Eq. 1 can be recast as I ( ν, n ) = Ψ ( ν, n )[ η ( ν, n )], (7)where Ψ is a matrix operator whose elements depend on thelevel populations (for more details see e.g. Hubeny & Mihalas2014, hereafter, HM14).This equation gives the solution of the intensity if all ele-ments of the Ψ operator, which are functions of the atomic pop-ulations, are known. However, in practice, building the full Ψ operator is never a ff ordable and the emissivity appearing in Eq.7 is substituted by an old value η † evaluated from a known esti-mate of the level populations (i.e. old values obtained with a pre-vious iteration). The new intensity obtained by the application ofthe operator on the old emissivity is then used to compute a newvalue of the level populations and of the emissivity, which areused to compute a new value of the intensity. This iterative pro-cess, between old and new values of the populations, is knownas Λ − iteration and it is repeated until convergence. Although Λ − iteration is very e ffi cient in an optically thin region, it suf-fers convergence problems in very optically thick regions. This It is common in radiative transfer problems to see the Λ operatorwhen it comes to " Λ − iteration". The relation between the Ψ and Λ op-erators is given in Rybicki & Hummer (1991) as Ψ = Λ /χ .Article number, page 2 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost drawback of classical Λ − iteration can be overcome by a slightmodification of Eq. 7 (see Cannon 1973a,b; Scharmer 1981; Ol-son et al. 1986), that is recasting Eq. 7 in I ( ν, n ) = ( Ψ ( ν, n ) − Ψ ∗ ( ν, n ))[ η † ( ν, n )] + Ψ ∗ ( ν, n )[ η ( ν, n )], (8)with an approximate operator Ψ ∗ ( ν, n ) built from a subset of theoriginal operator. Equation 8 improves the convergence of theclassical Λ -iteration method in optically thick regions and hasbeen named ALI for accelerated Λ -iteration (after the work ofHamann (1985) and Werner & Husfeld (1985)).In most applications, the diagonal of the full operator is used,although other studies recommended a tri-diagonal approximateoperator (see for instance Hennicker et al. 2018). Higher orderapproximate operators can improve the convergence of the ra-diative transfer problem at the price of computational time andmemory storage. For 2D radiative transfer problems, Auer et al.(1994) suggested to use the diagonal part of the operator in com-bination with an algorithm to accelerate the convergence (see§2.7). In our code, we use the diagonal part of the full operatoras a compromise between speed and accuracy. The level populations of an atom in a stellar plasma is a func-tion of the rates of transitions populating or depopulating a givenlevel. These rates in turn are a function of level populations, col-lisions, and radiation.The statistical equilibrium equations (SEE) for level (cid:96) isgiven by (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) ( C (cid:96) (cid:48) (cid:96) + R (cid:96) (cid:48) (cid:96) ) = n (cid:96) (cid:88) (cid:96) (cid:48)(cid:48) ( C (cid:96)(cid:96) (cid:48)(cid:48) + R (cid:96)(cid:96) (cid:48)(cid:48) ), (9)where n (cid:96) is the population of level (cid:96) , and C (cid:96) (cid:48) (cid:96) , and R (cid:96) (cid:48) (cid:96) the colli-sional and radiative rates for a transition between level (cid:96) (cid:48) and (cid:96) ,respectively. Unlike collisional rates, radiative rates are a product of the radia-tive transfer code. These rates can be expressed in terms of theradiative transfer coe ffi cients and are written as R (cid:96)(cid:96) (cid:48) = (cid:90) d Ω (cid:90) d ν h ν { U (cid:96)(cid:96) (cid:48) ( ν, n ) + V (cid:96)(cid:96) (cid:48) ( ν, n ) I ( ν, n ) } , (10)where the integration is carried out over frequency ν and all solidangles Ω . For line transitions these rates are written as R i j = B i j J R ji = A ji + B ji J , (11)where B i j , B ji , and A ji are the Einstein’s coe ffi cients, and J themean intensity integrated over the line given by J = (cid:90) d Ω π (cid:90) d ν I ( ν, n ) φ ( ν, n ), (12)where φ ( ν, n ) is the line absorption profile, also appearing inEq. 2. For continuum transitions they are as follows: R i j = π h (cid:90) d ν α i j ( ν ) J ( ν ) R ji = π h n ∗ i n ∗ j (cid:90) d ν α i j ( ν ) e − h ν/ kT { h ν c + J ( ν ) } , (13) We neglect both the advective and non-stationary terms as they arein general negligible for non-relativistic flows. where J ( ν ) is the mean radiation field given by J ( ν ) = (cid:90) d Ω π I ( ν, n ) (14) Starting from Eq. 9, passing the right-hand side term on the left-hand side and noticing that n l = (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) δ (cid:96)(cid:96) (cid:48) with δ (cid:96)(cid:96) (cid:48) non-zeroonly for (cid:96) (cid:48) = (cid:96) , we can factorise out (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) and write the SEE as (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) Γ (cid:96) (cid:48) (cid:96) =
0, (15)where the rate matrix Γ (cid:96) (cid:48) (cid:96) is Γ (cid:96) (cid:48) (cid:96) = R (cid:96) (cid:48) (cid:96) + C (cid:96) (cid:48) (cid:96) − δ (cid:96)(cid:96) (cid:48) (cid:88) (cid:96) (cid:48)(cid:48) { R (cid:96)(cid:96) (cid:48)(cid:48) + C (cid:96)(cid:96) (cid:48)(cid:48) } (16) Multilevel radiative transfer is a complex problem to solve be-cause of the dependence of the radiation field on the populations,and of the populations on the radiation field. Using Eqs.8 and 10the rate matrix in Eq. 16 for a multilevel atom reads: Γ (cid:96) (cid:48) (cid:96) = C (cid:96) (cid:48) (cid:96) + (cid:90) d Ω (cid:90) d ν h ν { U (cid:96) (cid:48) (cid:96) + V (cid:96) (cid:48) (cid:96) I e ff ( ν, n ) }− δ (cid:96)(cid:96) (cid:48) (cid:88) (cid:96) (cid:48)(cid:48) [ C (cid:96)(cid:96) (cid:48)(cid:48) + (cid:90) d Ω (cid:90) d ν h ν { U (cid:96)(cid:96) (cid:48)(cid:48) + V (cid:96)(cid:96) (cid:48)(cid:48) I e ff ( ν, n ) } ] − (cid:90) d Ω (cid:90) d ν h ν (cid:88) j > l ( n l V l j − n j V jl ) (cid:88) i <(cid:96) (cid:48) Ψ ∗ ( ν, n ) U (cid:96) (cid:48) i , (17)where I e ff ( ν, n ) = I † ( ν, n ) − Ψ ∗ ( ν, n ) η † ( ν, n ) is an e ff ective radi-ation field built with known quantities (old values). In Eq. 17we assumed that background opacities are constant within sub-sequent iterations.The system of equations 15 is nonlinear in the new popu-lations as it involves the product of the form ∝ n (cid:96) × n (cid:96) (cid:48) in Eq.17. The multilevel accelerated lambda iteration method (here-after, MALI) proposed by Rybicki & Hummer (1991), latter im-proved by Uitenbroek (2001), uses the operator splitting methodand a full preconditioning to make Eqs. 15 and 17 linear in thenew populations. The coupled equations of statistical equilib-rium with the radiation field are finally written as follows: (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) C (cid:96) (cid:48) (cid:96) + (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) (cid:90) d Ω (cid:90) d ν h ν { U † (cid:96) (cid:48) (cid:96) + V † (cid:96) (cid:48) (cid:96) I e ff ( ν, n ) }− (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) δ (cid:96)(cid:96) (cid:48) (cid:88) (cid:96) (cid:48)(cid:48) [ C (cid:96)(cid:96) (cid:48)(cid:48) + (cid:90) d Ω (cid:90) d ν h ν { U † (cid:96)(cid:96) (cid:48)(cid:48) + V † (cid:96)(cid:96) (cid:48)(cid:48) I e ff ( ν, n ) } ] − (cid:88) (cid:96) (cid:48) n (cid:48) (cid:96) (cid:90) d Ω (cid:90) d ν h ν (cid:88) j ( n † l V † l j − n † j V † jl ) (cid:88) i <(cid:96) (cid:48) Ψ ∗ ( ν, n ) U † (cid:96) (cid:48) i = Λ -iteration (Rybicki & Hummer 1992; Uitenbroek 2001) is called Article number, page 3 of 13 & A proofs: manuscript no. bt_spidi1 the cross-coupling term. From Eq. 4 the term (cid:88) j ( n † l V † l j − n † j V † jl )represents a summation over all opacities of all transitions with aplus sign if l < j or a minus sign if l > j , that is (cid:88) j > l χ † l j − (cid:88) j < l χ † jl . In some cases electron scattering is a dominant source of opac-ity and its e ff ect on both the continuum and spectral lines has tobe taken into account (Hillier 1991). However, the evaluation ofthe scattering emissivity involves the calculation of a scatteringintegral whose expression in the observer’s frame is not trivial(Rybicki & Hummer 1994). We evaluate the electron scatteringemissivity by setting it to the mean intensity (coherent approxi-mation), which is in turn determined iteratively through classical Λ -iterations, as suggested by Rybicki & Hummer (1992). Thismethod is accurate enough for the main range of applications ofthe code, but becomes inaccurate for hot stellar wind where theelectron density is very high and the ratio of scattering to trueabsorption is large. To do reliable radiative transfer simulations it is necessary toallow for line overlap and merging close to series limits. TheMALI formulation of the non-LTE problem is already capableof dealing with overlapping lines, and we include the formal-ism of occupation probabilities of Hummer & Mihalas (1988)to allow for lines to merge at the series limits. Because of thepresence of perturbers, each bound state i of an atom has a prob-ability w ( i ), with respect to the same state of a similar isolatedatom, to remain bound. Conversely, 1 − w ( i ) represents the prob-ability of state i to be dissolved, in other words, to belong tothe continuum states. Dissolved states contribute to a pseudo-continuum opacity beyond the series limits. Our implementationof the occupation probabilities formalism is similar to that ofHubeny 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 oflevel i to be undissolved and w ( j ) this probability for level j ,then if j is undissolved, i is necessarily undissolved as it lies ina lower energy state. Therefore, all upward rates in the statis-tical equilibrium equations have to be multiplied by w ( j ) / w ( i ),the probability that the upper level j is undissolved given i isundissolved. In §4.1.1 the impact of level dissolution on stellarcontinua is shown. To speed up the MALI method we implemented the method ofacceleration of Ng (1974) with general orders as defined in Aueret al. (1994) and Uitenbroek (2001). This method uses N sol + N start MALI iterations (i.e. withoutacceleration), we start accumulating N sol + N pending MALI iter-ations. In our tests, we use N sol = N start = N pending = N pending to 0 mainly results in extra over-head due to the matrix inversion required in the minimisationprocedure (see for instance Uitenbroek 2001).The main drawback of Ng’s acceleration for multidimen-sional models is the memory required to store N sol +
3. Implementation in MCFOST
The formulation of the non-LTE line transfer of Rybicki & Hum-mer (1992), although developed for 1D grids, can be easilyadapted for multidimensional models. Auer et al. (1994) appliedthe MALI method to a 2D Cartesian grid, while Hauschildt &Baron (2014) applied it to a 3D spherical grid. Recently, DeCeuster et al. (2020) applied the method to an unstructured grid.MCFOST-art, which stands for MCFOST atomic radiativetransitions, is an ensemble of modules implemented in MCFOST(Pinte et al. 2006, 2009). A 3D radiative transfer code, MCFOSTis written in modern Fortran, dedicated to the modelling of dustemission and molecular lines in the environment of young stars(e.g. protoplanetary discs, circumstellar envelopes). This codehas been coupled with SPH codes (e.g. Price et al. 2018) andMHD codes (e.g. Riols et al. 2020). Currently MCFOST handlesthree types of grids: a cylindrical grid dedicated to the modellingof discs, a spherical grid, and an unstructured grid based on aVoronoi tessellation (see e.g. Camps et al. 2013). The solutionof the line transfer equation is done by ray-tracing: several raysare propagated in specific directions along which Eq. 1 is solvedfor each wavelength simultaneously.In our new modules, we propagated rays from each cell cen-tre (spatial grid units in 2D and 3D geometries) in directions de-fined by the angular quadrature scheme (see below). Cells repre-sent the smallest resolution elements of the code and all quanti-ties are constant within them. If the model has a non-zero veloc-ity field, the velocity is projected in the direction of propagationof the ray. When velocity fields are present, it is important tomaintain a proper spatial discretisation to prevent nonphysicalfeatures due to large velocity gradients (Ibgui et al. 2013). TheMCFOST code detects if the projected velocity between two gridpoints is too large and linearly interpolates the projected velocityaccordingly.To produce images, we solved the intensity out of each pixelby integrating Eq. 1 along rays centred on each pixel.The solution of Eq. 1 along a ray, specified by a vector direc-tion n , is written as I ( ν, n ) = I b ( ν, n ) e − τ tot ( ν, n ) + (cid:88) k S k ( ν, n ) (1 − e − d τ k ( ν, n ) ) e − τ k ( ν, n ) ,(19)where I ( ν, n ) is the emerging specific intensity in a given direc-tion; 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 Article number, page 4 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost inner (entrance) boundary of a cell. In Eq. 19, the sum representsthe contribution of each cell, along the ray path, to the outgoingradiation.The term d τ k = τ k + − τ k = χ k s is the optical depth accumu-lated 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 + − s ( k ) in the cell is written as I k + ( ν, n ) = I k ( ν, n ) + S k ( ν, n ) (1 − e − d τ k ( ν, n ) ) e − τ k ( ν, n ) , (20)The approximate operator Ψ ∗ appearing in eqs. 8 and 18 isevaluated 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 Ψ ∗ k ( ν, n ) = (1 − e − d τ k ( ν, n ) ) / χ † k ( ν, n ), (21)where the numerator is the diagonal of the Λ operator, that is Λ ∗ k ( ν, n ) = (1 − e − d τ k ( ν, n ) ) (22)This satisfies the following conditions in the two extremecases: in optically thin regions it is 0, which corresponds to clas-sical Λ -iterations, and it approaches 1 in optically thick regions.Although this method allows for fast integration of the radiativetransfer equation, it depends on the numerical resolution, com-pared to other methods that use linear or cubic interpolationsbetween grid points.The part of the code that handles the propagation within thegrid is presented in Pinte et al. (2006). In the following, wepresent in some detail the solution of the non-LTE atomic lineformation problem in the code. In MCFOST-art, there are currently two choices for the initial-isation of the level populations of the non-LTE problem: pop-ulations are evaluated at LTE and populations are obtained bysolving the SEE with the radiation field set to zero. The choiceof the initial solution is the critical point of the non-LTE prob-lem. If the initial guess is too far from the solution, the con-vergence can be slowed down or even fail. In the eventualitythat none of these choices lead to convergence, we implementedthe collisional-radiative switching method of Hummer & Voels(1988). This method allows for a smooth transition between acollision dominated initial solution (i.e. at LTE) to a non-LTEinitial solution. Alternatively, the populations from a previouscalculation are given.
The absorption profile for each atomic bound-bound transitioncan either have a Gaussian or a Voigt shape. The choice of theshape of the absorption profile depends on the line and the prob-lem.The width of the absorption profile, v D , is given by the sumof the line thermal width and the microturbulence ξ as follows: v D = (cid:113) k b T / m + ξ , (23)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 dampingand Van Der Waals and Stark broadenings. The radiative damp-ing of a line transition j → i , Γ j , is computed by summing theEinstein A ji coe ffi cients for all (allowed) transitions for which j is an upper level as follows: Γ j = (cid:88) i < j A ji (24)In the latter expression, we neglected the e ff ect of stimulatedemission.Pressure broadening takes into account van der Waals andStark (linear and quadratic) broadenings computed in the Lind-holm theory. We computed interaction constants using the Un-söld formula (Unsöld & Weidemann 1955) and evaluated thequadratic Stark broadening following Gray (2008, Eq. 11.33).Linear Stark broadening was computed following Sutton (1978).The treatment of the Stark broadening, especially for hydrogenlines is approximate and is only accurate for the first lines of eachseries. In the future, we plan to add a more accurate treatment ofthe Stark broadening (see e.g. Stehlé & Hutcheon 1999). The wavelength grid used for non-LTE calculations was builtfrom individual atomic transitions. Firstly, we set up a grid ofcontinuum transitions at a moderate resolution. This continuumgrid was used to solve the continuum radiative transfer equation.Secondly, lines were gathered in groups depending on theirbounds (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 withthe radiation field. Since in moving atmosphere a line feels theradiation from di ff erent regions, the line bounds include the max-imum shift in frequency due to the velocity fields. However, thelocal (i.e. in the co-moving frame) line profile is not a ff ected byvelocity fields.Each group of lines was sampled with a constant resolutionin kilometers per second. Once wavelength grids for each groupwere determined, they were merged and continuum points fromthe continuum grid were added outside line groups. Opacity and emissivity of transitions were computed for eachatom using Eqs. 4, 2, and 3. The photoionisation cross-sectionof hydrogen-like ions was computed using Kramer’s formulawith correction from quantum mechanics (HM14, Eq. 7.92). Fornon-hydrogen-like ions, the photoionisation cross-section wasdirectly read from the atomic file and interpolated on the wave-length grid of MCFOST. We took the photoionisation data fromTOPbase (Cunto & Mendoza 1992).We also included continuum transitions from other sourceswhich add up in χ c . These are Thomson scattering on free elec-trons, hydrogen free-free (HM14, Eq. 7.100), and H − free-freeand bound-free transitions from John (1988). The main module of our code deals with solving Eq. 15 simul-taneously with Eq. 1 via the non-LTE loop. Figure 1 shows themain steps leading to the solution of the non-LTE problem with http://cdsweb.u-strasbg.fr/topbase/topbase.html Article number, page 5 of 13 & A proofs: manuscript no. bt_spidi1
MCFOST-art. We stress that even if Eq. 15 is solved locally (i.e.within each cell) it is coupled with other cells through Eq. 1.Firstly, the code reads abundances for the atomic species un-der consideration, the corresponding atomic data (e.g. energylevels E i and the oscillator strength f i j ), and an atmosphericmodel (temperature, densities, and velocity fields). The elec-tron density is evaluated at LTE if not present in the model (formore detail on electron density loop, see sect. 3.6). Secondly,MCFOST-art initialises the level populations of each atom asdescribed above.Finally, we evaluate continuous background opacities. Inmodels in which electron scattering is an important sourceof opacity, the continuum mean intensity and the correspond-ing Thomson scattering emissivity are computed with the ALImethod (for more details, see chapter 13 of HM14). Once back-ground opacities have been evaluated, the non-LTE loop starts.It computes a self-consistent solution for level populations andradiative transfer equations. For each cell, the rate matrix Γ (cid:96) (cid:48) (cid:96) isbuilt for a set of rays, for all wavelengths simultaneously. Oncethe rate matrix is known, a new value of the level populationsis determined. Electron densities can be evaluated with the newvalues of the non-LTE populations. To avoid convergence issues,the evaluation of the electron density can be done every N itera-tions of the non-LTE loop.The new populations are then compared to the previous pop-ulations using the following criterion of convergence: δ nn = (cid:107) − n † n (cid:107) < (cid:15) , (25)where (cid:15) is a user defined convergence threshold. The threshold (cid:15) is usually chosen between 10 − - 10 − . If the maximum rel-ative changes for each cell and each atom is below the conver-gence threshold the non-LTE loop stops. Otherwise, a new itera-tion starts with the previously computed populations. Excitationand ionisation temperatures for each transition of each atom arealso checked for convergence. In general, transition temperaturesconverge faster than level populations. Electron densities for a gas mixture of hydrogen and metalscannot be determined analytically. The populations of atomiclevels rely on the knowledge of the electron densities (for in-stance, through the Saha-Boltzmann equation at LTE). In turn,the electron density is a function of the population of atomiclevels through their ionisation fraction. Therefore, electron den-sities have to be determined iteratively. Starting from an initialguess for the electron density, the populations of atomic lev-els 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 arecomputed. This iterative procedure is repeated until the relativechange in electron densities drops below 10 − .In our code we use a linearisation procedure similar to Eq.17.88 of HM14. The initial guess for the electron density is givenby the ionisation of hydrogen plus one metal (M) n e = n e (H) + n e (M). Hydrogen ionisation yields n e (H) = ( (cid:112) n H φ H + − /φ H , (26)and the metal ionisation n e (M) = ( (cid:114) α M n H φ M +
14 (1 + α M ) −
12 (1 + α M )) /φ M , (27) where n H is the total number of hydrogen atoms, φ x the Saha-Boltzmann factor for the ion x, and α M is the abundance of themetal relative to hydrogen. The evaluation of the collisional rates appearing in the SEE (Eq.9) is cumbersome. It generally requires collision calculationswith quantum mechanics for which look-up tables of collisionalrates are available now (Barklem 2016). However, analytical andsemi-analytical recipes exist for faster evaluation of the collisionmatrix C ji (which is proportional to the collisional rates withina factor). We consider collisional excitation and ionisation byelectrons for hydrogen atoms using Johnson (1972). For colli-sional excitation of metals (and helium) by electrons we use vanRegemorter’s fomula (van Regemorter 1962) and the impact pa-rameter method (Seaton 1962). For ions, we follow Eq. 34 ofAller et al. (1982) and, for neutrals, we follow Seaton (1962).For collisional ionisation by electrons, we use Eq. 9.60 of HM14and for collisions with hydrogen atoms, we use Drawin’s for-mula (Eq. 9.62 of HM14). One of the cornerstones in solving the SEE is the evaluation ofthe integrals appearing in Eq. 10. In the following, we describehow we performed wavelength and angular integrations.Angular quadrature for molecular lines transfer in MCFOSTis done with the Monte Carlo method. For atomic lines, we useda set of fixed directions and weights uniformly distributed on theunit sphere such that integration over the sphere gives 4 π . Weused 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ˇepán et al. (2020) for un-polarised radiation. These two types of quadrature points givesimilar results. By default we used type A quadrature of Carlsonwith 80 points sampling 4 π steradians (Ibgui et al. 2013).De Ceuster et al. (2020) used a similar angular quadraturemethod but with the quadrature points and weights determinedby the HEALPix algorithm (Górski et al. 2005). For each cellof a model, starting from the centre of the cell (hereafter onepoint quadrature), the transfer equation is solved for all direc-tions, while the radiative rates and the mean intensity are accu-mulated. Integrating the transfer equation from the centre of eachcell may underestimate angular means and introduce some inac-curacy in the level populations because it only gives the meanintensity at that position. However, we do not know beforehandwhether this inaccuracy is significant or not. Thus, we imple-mented two options to deal with this issue: first, several startingpoints for the angular quadrature scheme can be selected for eachcell; and, second, Monte Carlo iterations can be done to properlysample each cell, in addition to the one point angular quadraturedescribed above. In §4.3, we discuss these two options.We replace the double integral appearing in Eq. 12 by a sumover wavelength points and directions as follows: (cid:90) d Ω π (cid:90) d ν I ( ν, n ) φ ( ν, n ) ≡ (cid:88) ray ω ray (cid:88) λ ω λ I ( λ, ray) φ ( λ, ray),(28)where ω ray is the weight of the angular quadrature for a givendirection (or ray) and ω λ the weight of a trapezoidal wavelengthquadrature. In the case of a line transition, ω λ includes the nor-malisation of the line profile. Article number, page 6 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost
Non-LTE loopSolve for RTE Solve for SEE convergence
Compute outgoing flux !"" < ε ? % & ' loop ( ), +; % - .(), +; % - ) 1 % Γ = 0 Read model • Atmospheric structure ( , T, etc…) • Abundances (log10(A - 12)) • atomic files ( : ,; :< , etc…) Set-up initial solution • LTE populations ( % ∗ ) • Wavelength grid • Background opacities (> ? ,@ ? ) Fig. 1: Flow chart of MCFOST-art. See §3.5 for more details. The green arrows indicate the main iterative loop. The double greyarrow indicates that the electron density loop is iterated within the main iterative loop and with the current estimates of the non-LTEpopulations.
4. Benchmarking the code
The MCFOST code is intrinsically 3D, but in this paper we fo-cus on 1D models that test the core of the method and allow usto validate our implementation with well-established 1D codes.While this is not the main application intended for our code, italso illustrates that MCFOST can now be used to generate a vari-ety of stellar spectra. We will devote a forthcoming paper to theapplications of our code to 2D and 3D geometries, with a par-ticular emphasis on stellar magnetospheres. We used the codesTURBOspectrum (Plez 2012) and RH (Uitenbroek 2001) as ref-erences for testing.The code MCFOST does not include a dedicated grid for 1Dmodels, and we used the spherical grid instead. However, thespherical grid of MCFOST is log-scaled in the radial directionand cannot e ffi ciently map the non-uniform grid from the 1Dmodels. Therefore, for our benchmarks, we directly read the gridfrom the input models as RH and TURBOspectrum do. In thatcase, a cell is defined as the distance between two consecutivegrid points. Because Eq. 1 is solved in MCFOST by summing upthe intensity of each cell (i.e. spatial grid unit) weighted by theoptical depth (see Eq. 19), the solution is more dependent on thegrid resolution (i.e. on the size of cells). On the contrary, RH andTURBOspectrum formal solvers typically interpolate the sourcefunction in the initial grid, resulting in a solution that is highlyaccurate. To mitigate the impact of a lack of resolution on theaccuracy of the solution, we decreased the size of the cells byadding more points between each point of the initial models.The TURBOspectrum code is used to model spectra fromcool evolved stars (Alvarez & Plez 1998; de Laverny et al. 2012) in synergy with the model atmosphere code MARCS (Gustafs-son et al. 2008). This code works only at LTE but includes a widelibrary of molecules, the occupation probabilities formalism, anda proper treatment of hydrogen line broadening. To match thecapabilities of our code, we slightly modified TURBOspectrumand removed molecular opacities, although, at the coolest tem-perature, neutral hydrogen density is dependent on the formationof molecules, such as H .The radiative transfer code RH is mainly used to model lineformation in non-LTE conditions using the MALI method. Re-cent versions of RH have been used to model spectral lines tak-ing into account the Zeeman e ff ect and partial frequency re-distribution in 3D cubes of the solar photosphere (de la CruzRodríguez et al. 2019). Recently, Criscuoli et al. (2020) bench-marked the RH code and showed its ability to model the solar ir-radiance with a good agreement with observations. The RH andMCFOST codes use similar treatment of background opacitiesand both employ the MALI method to solve for level popula-tions. We compared our treatment of background opacities andlevel dissolution at LTE with TURBOspectrum. We used RH forbenchmarking the MALI method on the non-LTE formation ofhydrogen lines. We used the (1D) spherically symmetric grid ofTURBOspectrum and RH for all tests. To test our opacity modules and background opacity calcula-tions, we used LTE models of stellar photospheres from the
Article number, page 7 of 13 & A proofs: manuscript no. bt_spidi1
Model ID T e ff (K) log g (cm s − ) radius (R (cid:12) ) ξ (km s − ) LTEsun 5777 4.44 1.00 1 yess8000 8000 2.50 9.34 1 yess3500 3500 3.50 2.95 2 yesC 5777 4.44 1.00 variable noTable 1: Stellar atmosphere models used for the benchmarking of our code. Models are taken from the MARCS database andcorrespond to the photosphere, except the model C taken from Fontenla et al. (1999). The first column indicates the designation ofthe model; the second, third, and fourth are the e ff ective temperature, log surface gravity, and radius of the star, respectively. Thefifth column gives the constant microturbulence, except for the FAL-C model, which has a depth varying microturbulence. The lastcolumn indicates if the model is computed assuming LTE.MARCS database (Gustafsson et al. 2008). The summary ofeach model is given in Table 1. The first model, sun , is thestandard solar photosphere model from the MARCS code. Thelast two models, s s We compared the stellar continuum computed by the three codesfor the stellar photosphere models (see Table 1). Figure 2 showsthe continuum flux and the discrepancy function δ F / F for eachmodel.The agreement between RH and MCFOST is excellent, be-ing at or below 1% at most wavelengths. For models sun (Fig.2a) and s s ff erences are larger (again owing to the steep slope of the emis-sion), but remain below 30%. Discrepancies between MCFOSTand TURBOspectrum are of the same order as those betweenMCFOST and RH. The largest discrepancy for the model s − continuum. In TURBOspec-trum, the H − density is self-consistently determined by solvingchemical equilibrium equations including hydrogen in its atomicand molecular forms, leading to a slightly di ff erent value thanin RH and MCFOST. For the model s , but not in MC-FOST, thereby impacting the shape of the continuum.The agreement between MCFOST and TURBOspectrum re-mains excellent when level dissolution is taken into account (Fig.2 right panels). The impact of level dissolution is clearly visi-ble at the Balmer jump where it smoothes the abrupt discontinu-ity towards redder wavelengths. The shallower spectrum results https://marcs.astro.uu.se/ Continuum scattering cannot be removed easily in TURBOspectrumto benchmark further the codes, but it is negligible for all models butthe hotter model. in an improved agreement between TURBOspectram and MC-FOST for model s Although the impact of velocity fields on line formation will beaddressed in a forthcoming paper, in this section we show anexample of a solar photosphere with a radial velocity V r givenby V r = v + ( V inf − v ) (1 − r / r ) β , (29)where v is the velocity at the inner point (bottom of thephotosphere), V inf a velocity such that max(V r ) =
500 km s − ,and β the exponent of the velocity law. In our calculations,v = − and β = . ff erence between MCFOST and RH forthe solar photosphere when a velocity field is included. Whilethe continuum radiation is similar to that of Fig. 2a, lines are im-pacted significantly by the large velocity gradient, as expected.The di ff erences between the two codes are negligible, althoughthe lack of resolution in RH is noticeable. [nm] F l u x [ W m m ] [nm] F l u x [ W m m ] [nm] F l u x [ W m m ] Fig. 3: Flux for a solar photospheric model with velocity gradientgiven by Eq. 29. The flux from MCFOST is shown in black andfrom RH in cyan. A zoom-in around Ly β and H α is also shown. Article number, page 8 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost F l u x [ W m m ] RHTURBOspectrumMCFOST
200 300 400 500 600 700 800 [nm] F / F ( % ) (a) Solar photosphere without level dissolution. F l u x [ W m m ] TURBOspectrumMCFOST
200 300 400 500 600 700 800 [nm] F / F ( % ) (b) Solar photosphere with level dissolution. F l u x [ W m m ] RHTURBOspectrumMCFOST
250 500 750 1000 1250 1500 1750 2000 [nm] F / F ( % ) (c) s F l u x [ W m m ] TURBOspectrumMCFOST
250 500 750 1000 1250 1500 1750 2000 [nm] F / F ( % ) (d) s F l u x [ W m m ] RHTURBOspectrumMCFOST
500 1000 1500 2000 2500 3000 [nm] F / F ( % ) (e) s F l u x [ W m m ] TURBOspectrumMCFOST
500 1000 1500 2000 2500 3000 [nm] F / F ( % ) (f) s Fig. 2: Continuum flux as a function of wavelengths for the photospheric models without (left panels) and with (right panels) leveldissolution. 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 di ff erence between RH and MCFOST (thick dark greylines), and between TURBOspectrum and MCFOST (dashed light grey lines). As RH does not include level dissolution, it is notshown in the right figures. In figure 2c, the continuum flux computed by RH stops at about 800 nm, leading to incorrect evaluationof the error at this edge point.We also tested the treatment of the same model with no radialvelocity but a constant rotational velocity of 300 km s − instead.Figure 4 shows the H α line flux as a function of inclination. As expected, the impact of rotation on the line shape increases withinclination, from no e ff ect when the star is seen pole-on to max-imum broadening when the star is equator-on. Article number, page 9 of 13 & A proofs: manuscript no. bt_spidi1 [nm] F l u x [ W m m ] Fig. 4: H α flux variation with inclination for a solar photosphericmodel with rotation. The flux at inclinations of 0 ◦ , 45 ◦ , and 90 ◦ are shown in black, green, and red, respectively. An inclinationof 0 ◦ means that the star is seen pole on. Although a LTE model of the solar photosphere is adequate tomodel most solar lines, it fails to reproduce lines formed beyondthe photosphere. The upper solar atmosphere is a high temper-ature and low-density region above the photosphere where non-LTE e ff ects dominate the line formation. In the following wefocus on the Lyman α (Ly α ) and H α lines as they form in theupper atmosphere (e.g. Figure 1 of Vernazza et al. 1981) and arethus impacted by non-LTE e ff ects.We used a 6-level hydrogen atom, with five bound levelsand one continuum level, which is the ground state of H II. En-ergy levels and transition frequencies are from Johnson (1972).Bound-free cross-sections are computed using Kramer’s formulaand collisional rates are evaluated following Johnson (1972). Forthe Lyman α and H α lines, the absorption profiles are given bya Voigt function. For the other lines, we used a Gaussian profile.We used model C of Fontenla et al. (1999) as a standard model ofthe solar atmosphere including the chromosphere and transitionregion.Overall, the agreement between RH and MCFOST is excel-lent. Figure 5 shows the ratio b between non-LTE and LTE pop-ulations in model C. In the photosphere (below the temperatureminimum), the populations are mostly at LTE with a b ratio of afew. In the deep photosphere, b is close to 1. From the temper-ature minimum to the transition region (the chromosphere), thepopulations start to depart from their LTE values. In the lowercorona, the populations are very far from LTE, with a b factorgreater than 10 for the ground state of HI (level 1 in Fig. 5).The formation regions of the two lines are consistent with theearlier work of Vernazza et al. (1981). Contribution functions(i.e. regions of formation) of the Ly α and H α lines are shown inFigure 6.Figure 7 shows the solar flux computed in model C with RHand MCFOST, at non-LTE (a) and LTE (b). The impact of non-LTE e ff ects on spectral lines is clearly visible. At LTE all linesare in emission, whereas at non-LTE most of the lines are seenin absorption. For comparison, disc-integrated observations ofthe Sun show that most of the hydrogen lines are in absorption,except for strong chromospheric lines (e.g. Ly α ). F l u x [ W m m ] RHMCFOST
200 400 600 800 1000 1200 1400 [nm] F / F ( % ) (a) Non-LTE flux. F l u x [ W m m ] RHMCFOST
200 400 600 800 1000 1200 1400 [nm] F / F ( % ) (b) LTE flux. Fig. 7: Solar flux for the model C. RH is shown with the thincyan and MCFOST with the thick black line. Below each panelthe discrepancy between the two codes is shown.
Article number, page 10 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost [kg m ] n o n - L T E t o L T E p o p . r a t i o b s o l a r r a d i u s h = 0 km t e m p e r a t u r e m i n i m u m h = 463 km t r a n s i t i o n r e g i o n h = 1877 km level Fig. 5: Ratio of non-LTE to LTE populations (departure coe ffi cient b) of a 6-level hydrogen atom as a function of density ( ρ ). RHpopulations are indicated with cyan dots. MCFOST populations are shown with full colour lines, where darker (lighter) red shadesmean lower (higher) energy levels. For clarity, the first three atomic levels and the continuum level (ground state of HII) are labelledexplicitly in the insert. The di ff erent layers of the solar atmosphere are shown with vertical dashed lines. Typical heights above thesolar radius are also indicated. [nm] [ k g m ] solar radiustemperature minimumtransition region d I / d h [ W m H z s r m ] (a) H α line. [nm] [ k g m ] solar radiustemperature minimumtransition region d I / d h [ W m H z s r m ] (b) Ly α line. Fig. 6: Line contribution function, dIdh (response function to height variation), at disc centre. The contribution function is drawn as afunction of the density in the atmosphere and distance from line centre (in nm). Lighter areas show where the contribution functionpeaks. We indicate the di ff erent layers of the solar atmosphere in each panel. Article number, page 11 of 13 & A proofs: manuscript no. bt_spidi1 F l u x [ W m m ] [nm] F / F ( % ) (a) Ly α line. F l u x [ W m m ] [nm] F / F ( % ) (b) H α line. Fig. 8: Zoom-in on solar lines flux. The discrepancy between RHand MCFOST is shown below each panel.At all wavelengths, the discrepancy is below 1%, except inthe core of some lines where it reaches 5%. In the figures, in-terpolation errors create glitches in the discrepancy functions.However, real di ff erences exist in the core of strong lines, suchas the Ly α line, where the mismatch can reach 20%-25%. Fig-ure 8 shows the non-LTE flux of the Ly α line (top panel) and theflux of the H α line (bottom panel). Although these di ff erencesare not critical, using one of the two options discussed in §3.8to improve the angular quadrature leads to a better agreementbetween the Lyman α line computed by RH and MCFOST. In §3.8, we presented the method we used to perform angularquadratures. These quadratures are important since they are usedto evaluate the radiative rates, which in turn determine the levelpopulations. As stated earlier, angular quadratures are performedfor a set of fixed rays (i.e. directions) starting from each cell cen-tre. In the following, we test option (i) of §3.8, that is performingthe angular quadrature not only from cell centres but also at Nrandom positions in the cells. As in section 4.2, we compute thenon-LTE Ly α line flux in model C. Further, we use the origi-nal grid of the 1D model, that is without additional points. Werandomly chose 99 additional positions in each cell. We thenkept these positions fixed during the non-LTE loop to avoid in-troducing Monte Carlo noise in the solution. For each of these 100 positions, the angular quadrature was performed for 80 raysresulting in a total of 8 000 rays. Thus, for each cell, we ob-tained 100 values of the mean intensity and radiative rates, eachat a di ff erent location inside the cell. The mean intensity and theradiative rates in each cell is just the arithmetic mean of thesevalues.Finally, we interpolated the level populations on the same gridused in §4.2 (i.e. with a higher resolution). Because now we havea better estimation of the radiative rates for each cell, the errorin the numerical integration of Eq. 1 dominates the discrepancybetween RH and MCFOST.In figure 9 we show the Ly α flux computed with this methodcompared to the profile shown in in Fig. 8a. The agreement be-tween RH and MCFOST is substantially improved, with a max-imum discrepancy below 5%, five times better than with onlyone position for the angular quadrature. Still we note that, evenif the core is better reproduced, the wings are not, although themaximum error is small.Similar results are obtained by performing a purely MonteCarlo step after the one point quadrature, that is we randomlyselected 1000 positions and directions (§3.8, option (ii)). Themain di ff erence between §3.8 option (i) and §3.8 option (ii) liesin the execution time of the code. The latter approach is eighttimes faster than the former for an equivalent accuracy. However,the populations computed with §3.8 option (ii) have Monte Carlonoise. F l u x [ W m m ] option (ii)option (i)RHMCFOST [nm] F / F ( % ) Fig. 9: Evolution of Ly α flux with angular quadrature methods.The black and cyan lines are the flux computed by MCFOST andRH, respectively, from Fig. 8a. The dotted red line and the pinktriangles are the flux computed using option (i) and (ii) of §3.8,respectively. The bottom panel shows the discrepancy betweenRH and MCFOST, from Fig. 8a (dark grey line), option (i) (lightgrey dotted line), and option (ii) (grey triangles).
5. Conclusions
In this paper, we tested our code against spherically symmetricmodels of stellar atmospheres, representing a range of physicalplasma conditions. We tested background opacity calculationswith lines merging close to series limits using the occupationprobabilities formalism. The agreement between MCFOST-artand both TURBOspectrum and RH is excellent, of the order of1% or below at most wavelengths. Still, a few discrepancies arepresent, arising from the di ff erent treatment of opacities at a fewcritical wavelengths (e.g. at the Balmer jump or at the H − min-imum). The code has also been successfully tested with vertical Article number, page 12 of 13. T essore , C. P inte , J. B ouvier , and F. M´ enard : atomic line mcfost and rotational velocities. The proper behaviour of all spectrallines is recovered in models with non-zero velocity fields.We performed non-LTE calculations of the solar Lyman α and H α lines formation. These calculations were compared withcalculations performed with the RH code. Line formation re-gions and non-LTE e ff ects are reproduced well by our code com-pared to earlier studies and observations. While the discrepancyin the flux between RH and MCFOST can reach 25% in the coreof Ly α , the discrepancy drops below 1% around the H α line andfor most other wavelengths in the spectrum. Furthermore, weshowed that improving the accuracy of the angular quadraturescheme decreases this discrepancy by a factor of five.Although we only tested the code against 1D models, theopacity modules and the non-LTE loop are geometry indepen-dent. These modules are therefore applicable to multidimen-sional models. The MCFOST code is 3D in essence and the so-lution of the radiative transfer equation is fully performed in 3D.Further applications of this new code to 2D and 3D geometrieswill be discussed in a forthcoming paper with a particular focuson stellar magnetospheres. Acknowledgements.
We are thankful to Dr. Bertrand Plez for providing andmaintaining a publicly available version of TURBOspectrum. B Tessore is grate-ful to Dr. Han Uitenbroek for providing a version of RH. We thank the anony-mous referee for his comments, which helped to improve the manuscript. B.Tessore is also thankful to the french minister of Europe and foreign a ff airsand the minister of superior education, research and innovation (MEAE andMESRI) for research funding through FASIC partnership. This project has re-ceived funding from the European Research Council (ERC) under the EuropeanUnion’s Horizon 2020 research and innovation programme (grant agreement No742095; SPIDI : Star-Planets-Inner Disk-Interactions, ). C. Pinte acknowledges funding from the Australian Research Council viaFT170100040 and DP180104235.
References
Alencar, S. H. P., Bouvier, J., Walter, F. M., et al. 2012, A&A, 541, A116Aller, L. H., Appenzeller, I., Baschek, B., et al., eds. 1982, Landolt-Börnstein:Numerical Data and Functional Relationships in Science and Technology- New Series “ Gruppe / Group 6 Astronomy and Astrophysics ” Volume 2Schaifers / Voigt: Astronomy and Astrophysics / Astronomie und Astrophysik“ Stars and Star Clusters / Sterne und Sternhaufen, 54Alvarez, R. & Plez, B. 1998, A&A, 330, 1109Auer, L., Bendicho, P. F., & Trujillo Bueno, J. 1994, A&A, 292, 599Barklem, P. S. 2016, Phys. Rev. A, 93, 042705Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233Camps, P., Baes, M., & Saftly, W. 2013, A&A, 560, A35Cannon, C. J. 1973a, J. Quant. Spectr. Rad. Transf., 13, 627Cannon, C. J. 1973b, ApJ, 185, 621Criscuoli, S., Rempel, M., Haberreiter, M., et al. 2020, Sol. Phys., 295, 50Cunto, W. & Mendoza, C. 1992, Rev. Mexicana Astron. Astrofis., 23, 107Dappen, W., Anderson, L., & Mihalas, D. 1987, ApJ, 319, 195De Ceuster, F., Homan, W., Yates, J., et al. 2020, MNRAS, 492, 1812de la Cruz Rodríguez, J., Leenaarts, J., Danilovic, S., & Uitenbroek, H. 2019,A&A, 623, A74de Laverny, P., Recio-Blanco, A., Worley, C. C., & Plez, B. 2012, A&A, 544,A126Fontenla, J., White, O. R., Fox, P. A., Avrett, E. H., & Kurucz, R. L. 1999, ApJ,518, 480Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759Gray, D. F. 2008, The Observation and Analysis of Stellar PhotospheresGustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951Hamann, W. R. 1985, A&A, 148, 364Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669Hauschildt, P. H. & Baron, E. 2014, A&A, 566, A89Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140Hillier, D. J. 1991, A&A, 247, 455Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151Hubeny, I. & Mihalas, D. 2014, Theory of Stellar AtmospheresHummer, D. G. & Mihalas, D. 1988, ApJ, 331, 794Hummer, D. G. & Voels, S. A. 1988, A&A, 192, 279Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 John, T. L. 1988, A&A, 193, 189Johnson, L. C. 1972, ApJ, 174, 227Kurosawa, R., Romanova, M. M., & Harries, T. J. 2011, MNRAS, 416, 2623Lima, G. H. R. A., Alencar, S. H. P., Calvet, N., Hartmann, L., & Muzerolle, J.2010, A&A, 522, A104Muzerolle, J., Calvet, N., & Hartmann, L. 2001, ApJ, 550, 944Ng, K. C. 1974, J. Chem. Phys., 61, 2680Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spectr. Rad. Transf.,35, 431Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797Plez, B. 2012, Turbospectrum: Code for spectral synthesisPrice, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95Rybicki, G. B. & Hummer, D. G. 1978, ApJ, 219, 654Rybicki, G. B. & Hummer, D. G. 1991, A&A, 245, 171Rybicki, G. B. & Hummer, D. G. 1992, A&A, 262, 209Rybicki, G. B. & Hummer, D. G. 1994, A&A, 290, 553Scharmer, G. B. 1981, ApJ, 249, 720Seaton, M. J. 1962, Proceedings of the Physical Society, 79, 1105Sobolev, V. V. 1957, Soviet Ast., 1, 678Stehlé, C. & Hutcheon, R. 1999, A&AS, 140, 93Sutton, K. 1978, J. Quant. Spectr. Rad. Transf., 20, 333Uitenbroek, H. 2001, ApJ, 557, 389Unsöld, A. & Weidemann, V. 1955, Vistas in Astronomy, 1, 249van Regemorter, H. 1962, ApJ, 136, 906van Zadelho ff , G. J., Dullemond, C. P., van der Tak, F. F. S., et al. 2002, A&A,395, 373Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635Štˇepán, J., Jaume Bestard, J., & Trujillo Bueno, J. 2020, A&A, 636, A24Werner, K. & Husfeld, D. 1985, A&A, 148, 417, G. J., Dullemond, C. P., van der Tak, F. F. S., et al. 2002, A&A,395, 373Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635Štˇepán, J., Jaume Bestard, J., & Trujillo Bueno, J. 2020, A&A, 636, A24Werner, K. & Husfeld, D. 1985, A&A, 148, 417