Issue 
A&A
Volume 653, September 2021



Article Number  A65  
Number of page(s)  22  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202140275  
Published online  10 September 2021 
MPSATLAS: A fast allinone code for synthesising stellar spectra
^{1}
Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: witzke@mps.mpg.de
^{2}
Blackett Laboratory, Imperial College London, South Kensington Campus, SW7 2AZ London, UK
^{3}
School of Space Research, Kyung Hee University, 701 Yongin, Gyeonggi, Republic of Korea
^{4}
Center for Astrophysics Harvard & Smithsonian, 60 Garden Street, 02138 Cambridge, MA, USA
Received:
1
January
2021
Accepted:
28
May
2021
Context. Stellar spectral synthesis is essential for various applications, ranging from determining stellar parameters to comprehensive stellar variability calculations. New observational resources as well as advanced stellar atmosphere modelling, taking three dimensional effects from radiative magnetohydrodynamics calculations into account, require a more efficient radiative transfer.
Aims. For accurate, fast and flexible calculations of opacity distribution functions (ODFs), stellar atmospheres, and stellar spectra, we developed an efficient code building on the wellestablished ATLAS9 code. The new code also paves the way for easy and fast access to different elemental compositions in stellar calculations.
Methods. For the generation of ODF tables, we further developed the wellestablished DFSYNTHE code by implementing additional functionality and a speedup by employing a parallel computation scheme. In addition, the line lists used can be changed from Kurucz’s recent lists. In particular, we implemented the VALD3 line list.
Results. A new code, the Merged Parallelised Simplified ATLAS, is presented. It combines the efficient generation of ODF, atmosphere modelling, and spectral synthesis in local thermodynamic equilibrium, therefore being an allinone code. This allinone code provides more numerical functionality and is substantially faster compared to other available codes. The fully portable MPSATLAS code is validated against previous ATLAS9 calculations, the PHOENIX code calculations, and highquality observations.
Key words: stars: atmospheres / stars: latetype / radiative transfer / opacity / convection
© V. Witzke 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.
Open Access funding provided by Max Planck Society.
1. Introduction
Electromagnetic radiation emitted from the stellar photosphere is one of the key sources of information about a star. To advance our understanding of stars, accurate spectroscopic and photometric measurements, as well as accurate modelling that allows us to connect the measured stellar electromagnetic radiation to properties of the stellar atmospheres and stellar interiors are essential.
The data gathered by the numerous groundbased and spaceborne telescopes that started observing within the last decade have highlighted the need for accurate modelling of stellar atmospheres and their spectra. For example, the Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO, see Pepe et al. 2010) and the High Accuracy Radial Velocity Planet Searcher (HARPS, see Mayor et al. 2003) made highresolution spectral data for thousands of stars available, while the Large Sky Area MultiObject Fibre Spectroscopic Telescope (LAMOST, see Ai et al. 2016) provides lowresolution spectra for millions of stars. The interpretation of these spectral measurements requires the modelling of stellar atmospheres on a fine grid of stellar fundamental parameters, such as effective temperature, surface gravity, and chemical composition.
Furthermore, the advent of planetary hunting missions, for example Kepler (Borucki et al. 2010), TESS (Ricker et al. 2014), CHEOPS (Benz et al. 2021), and WASP (Pollacco et al. 2006), has brought measurements of the photometric variability for several hundred thousand stars. Even more data are expected from the forthcoming PLATO mission (Rauer et al. 2014). The main source for photometric variability of cool stars, such as the Sun, are surface magnetic fields that affect the local structure in stellar atmospheres. Consequently, interpreting stellar photometric data requires an assessment of the effect of the magnetic field in stellar atmospheres on the emergent radiation.
There are different approaches to modelling stellar spectral and photometric fluxes. One of the simplest and most widely used approach relies on onedimensional (1D) modelling of stellar atmospheres under the assumption of radiativeconvective equilibrium (with a simple parameterisation for convective flux and overshooting). While such a 1D approach has a number of shortcomings (see, e.g., Koesterke et al. 2008; Uitenbroek & Criscuoli 2011), it proved itself to be an invaluable tool for various applications (see, e.g., Castelli & Kurucz 1994; Claret 2000; Mészáros et al. 2012; Marfil et al. 2020) and is extensively used in stellar physics.
A more comprehensive approach relies on three dimensional (3D) hydrodynamic and magnetohydrodynamic (HD and MHD, respectively) simulations of nearsurface convection in stars (see, e.g., Nordlund et al. 2009; Stein 2012; Freytag et al. 2012; Magic et al. 2013; Beeck et al. 2015b,a). Using the simulated 3D cubes, the emitted radiation can be calculated following a 1.5D approach, that is along many rays passing through such a 3D cube (see, Asplund et al. 2006; Riethmüller et al. 2014; Norris et al. 2017, for a detailed description of the 1.5D approach). An important advantage of the 3D MHD simulations and 1.5D approach over 1D radiative equilibrium modelling is that it allows us to directly account for the effects of the magnetic field on the emergent radiation and, consequently, to model the stellar spectral and photometric variability.
Altogether, there are two separate but similar challenges that demand comprehensive and fast spectral synthesis on an adjustable frequency grid, and broad band spectral intervals. First, the 1.5 dimensional (1.5D) modelling needs a huge amount of fast radiative transfer (RT) computations on 1D structures. Second, accurate and fast 1D atmosphere modelling with subsequent RT calculations for any stellar parameter is needed. The aim of the present study is to develop a fast and easily applicable RT code to address both of these challenges.
The accurate treatment of the line opacity poses the main challenge to spectral synthesis over broad spectral ranges because of the immense number of spectral lines, where current line lists contain up to hundreds of millions of lines (Kurucz 2005b) that need to be taken into account. A careful treatment of both atomic and molecular lines is imperative because atomic and molecular lines are interspersed in particular in the spectra of cool stars. Spectral lines not only dominate some spectral regions (e.g., the ultraviolet), but also affect the atmospheric structure by blocking photons.
While a forward spectral synthesis on a highresolution wavelength grid, that is typically with a resolving power, R = 500 000, is computationally expensive, not all applications require highresolution spectra. To that end, different techniques were developed to correctly account for line opacity, but reduce the computational cost on coarse resolution grids. The most commonly used methods are the opacity distribution functions (ODF) method and the Opacity Sampling (OS) method (Carbon et al. 1984; Castelli 2005a). Both methods approximate the opacity using different ways of sampling it. We preferred the ODF method for developing the RT code, as it results in significantly faster RT calculations compared to the OS method, thus making it more suitable for the 1.5D approach. Recently, the ODF approach was further optimised to make it more efficient. Cernetic et al. (2019) found the best configurations of the ODF to reach several times faster computations while maintaining accuracy. Moreover, the ODF approach was extended to calculate stellar fluxes as they are observed in various filters (Cernetic et al. 2019; Anusha et al. 2021).
Widely used RT codes used for spectral synthesis of cool stars include the local thermodynamic equilibrium (LTE) codes MARCS (Gustafsson et al. 2008), and MAFAGSOS (Grupp 2004), and the nonLTE codes PHOENIX (Husser et al. 2013), and TLUSTY (Hubeny & Lanz 2017). The PHOENIX code was developed to account for expanding atmospheres and deviations from LTE, and is therefore more complex and CPU intensive. Another successful RT code in LTE, which can calculate both model atmospheres in radiative equilibrium and the emergent spectra, is the ATLAS code by Kurucz (Kurucz 1970). Since for most applications mentioned above the spectral range does not include the extreme and far ultraviolet (UV), it is sufficient to consider LTE codes. Thus, the ATLAS code includes all crucial physics for radiative transfer in the atmospheres of mainsequence stars while keeping the setup simple. Consequently, we chose to build on the ATLAS code. While the ATLAS code comes in two versions, ATLAS12 (Castelli 2005a) and ATLAS9, we prefer the ATLAS9 version as it uses ODF whereas the ATLAS12 works with the OS method. The advantage of the ATLAS9 code is the short time required to compute a single model (a few minutes on a single core using the standard ODF table with 328 bins).
When updating ATLAS9 it is essential to also consider the DFSYNTHE code (Kurucz 2005a; Castelli 2005b). This code computes ODF tables for the ATLAS9 code. The main disadvantage of the ODF so far has been the limitation to the chemical composition and microturbulence velocities for which the ODF tables were pretabulated. The generation of ODFs using the DFSYNTHE code was not suitable for massive computations as several routines had to be successively executed, and the computation time was very long (Kurucz 2005a; Castelli 2005b). To achieve our goal of fast spectral synthesis for arbitrary abundances, we need to eliminate this bottleneck, and in addition make the ODF generation more user friendly. This is achieved by merging the DFSYNTHE code, which calculates highresolution opacities and uses them to obtain ODF, with the ATLAS9 code that can calculate both the atmospheric structures as well as the emergent spectra. Additionally, a more flexible treatment of the ODF calculations was implemented. The frequency resolution and ODF configuration can be changed, and a spectral filter included. Furthermore, we parallelised the code to speed up the highresolution opacity calculations and model calculations. Finally, the line lists used can be exchanged and we give example calculations with both, Kurucz’s line list as well as the most up to date VALD3 line list (Ryabchikova et al. 2015). In this paper we describe the structure and improvements of the resulting code, which we call the Merged Parallelised Simplified ATLAS code (MPSATLAS). Moreover, we validate MPSATLAS against previous ATLAS9 calculations, PHOENIX code calculations, and observations. The MPSATLAS code is available on request and it will be made publicly available soon.
The paper is structured as follows. A brief summary of the physics and its implementation is given in Sect. 2 and Sect. 3, respectively. All code improvements are listed in Sect. 4, where we also discuss limitations of the code. In Sect. 5 we comment on the code performance and test the resulting emergent spectra by a codetocode comparison and a code to observations comparison. The outcome of the work is summarised in Sect. 6, where also conclusions are drawn.
2. Radiative transfer calculations
The energy transport in stars defines the structure of their interior and atmosphere, as well as the emitted electromagnetic radiation. Focusing on the upper layers around the optical surface of a star, the radiative transport of energy is dominant. Thus, the imperative problem for modelling stellar atmospheres is solving the radiation transfer equation (RTE)
along a ray for a timeindependent system, where the subscript λ denotes the wavelength and implies that all quantities are monochromatic. I_{λ} is the intensity, S_{λ} is the source function, and μ = cos θ, where θ is the angle between the viewing direction and the normal to the stellar surface. The optical depth τ_{λ} is determined from dτ_{λ} = −χ_{λ}(s) ds, where s is the height in the atmosphere. The total extinction coefficient per unit volume is χ_{λ} = α_{λ} + σ_{λ}, which contains the absorption coefficient, α_{λ}, and the scattering coefficient, σ_{λ}. The source function is generally defined as the ratio of the total emission coefficient, j_{λ}, to the total extinction coefficient
Solving the RTE becomes straightforward once the emission coefficient and opacity are known along the atmosphere. In LTE and under the assumption of coherent isotropic scattering, the emissivity can be expressed as
where B_{λ} is the Planckfunction and J_{λ} is the mean intensity (Mihalas & Mihalas 1984). Then, calculating the opacity becomes the keystone to solving the RTE (except scattering coefficients).
2.1. Calculating opacity
The opacity κ_{λ} (hereafter, we refer to the opacity normalised per unit of mass κ_{λ} ≡ χ_{λ}/ρ, where ρ is the density) can be decomposed into continuum opacity, κ_{λ, c}, associated with different processes that involve atomic and molecular transitions with nondiscrete wavelength (i.e., boundfree and freefree transitions), and line opacity, κ_{λ, l}, due to discrete transitions, namely atomic and molecular spectral lines. To determine the total opacity, κ_{λ} = κ_{λ, c} + κ_{λ, l}, at each point, the atomic and molecular level populations have to be determined. Under the assumption of LTE, the atomic and molecular level populations in equilibrium can be calculated using the SahaBoltzmann (SB) equation. This implies that they depend only on elemental composition (i.e. abundances), local temperature, and pressure. In addition, in stellar atmospheres the Doppler broadening of lines occurs due to turbulence, which in 1D models is taken into account using the microturbulence parameter, v_{turb}. Hence, for a given composition, opacity is a function of wavelength, local temperature, pressure, and microturbulence κ ≡ κ(λ, T, P, v_{turb}).
In addition to the temperature and pressure values, calculating the level populations using the SB equation requires knowledge of the electron number density, n_{e}, which however is not known a priori. Thus under the assumption of particle conservation, the SB equation has to be solved iteratively. For a detailed description of the solver implemented in MPSATLAS see Appendix A and Kurucz (1970).
Having determined the populations, the continuous opacity, κ_{λ, c}, and the line opacity, κ_{λ, l}, can be calculated. For the continuous opacity, the MPSATLAS code takes the following contributors into account: Freefree (ff) and boundfree (bf) transitions in H^{−}, H_{2}, He, He^{−}, C, N, O, Ne, Mg, Al, Si, Ca, Fe, the molecules CH, OH and NH, and their ions. While for the calculation of the equilibrium number densities and the line opacity also C_{2} and CN are included, their contribution to the continuum opacity via photodissociation is neglected. Moreover, electron scattering and Rayleigh scattering on H I, He I, and H_{2} were considered.
The line opacity is more expensive to compute, simply because of the huge number of lines for which the line absorption coefficient needs to be calculated. For atomic and molecular transitions from an initial energy level, denoted by i to a final energy level j, the line absorption is defined as:
where ν_{0} is the frequency corresponding to the transition, e it the elementary charge, m_{e} is the electron rest mass, c is the speed of light, g_{i} is the statistical weight of the level i, f_{ij} is the dimensionless oscillator strength of the transition , and Δν_{D} is the Doppler width. The gas density is ρ, the temperature is T, and E_{i} is the energy of the initial energy level. The number density over partition functions for the entire ionisation stage is given by N_{k}/U_{k}, where the subscript k indicates the ionisation stage. The term in brackets in Eq. (4) describes the correction for stimulated emission, where the Boltzmann constant, k_{B}, and the Planck constant, h, have their usual designations. For the line broadening the Voigt function H(Δν/Δν_{D}, γ/4πΔν_{D}), where γ is the total damping constant and Δν = ν − ν_{0}, is used. While for most lines the Voigt profile is used, for hydrogen lines more accurate profiles are employed, namely Stark profiles are taken into account (Cowley & Castelli 2002).
2.2. Atmosphere models in radiative equilibrium and emergent spectra
Modelling a stellar atmosphere for a particular set of stellar fundamental parameters and elemental composition involves an iterative process of recalculating the atmospheric structure until radiative equilibrium (RE) is reached (see for example Collins 1989; Hubeny & Mihalas 2015, and references therein). In this process, in each iteration the RTE for a currently estimated atmospheric structure has to be solved in order to determine by how much the model structure needs to be corrected to satisfy the RE (for convergence criteria see Appendix B). The iterations are usually initialised with a starting structure that represents a converged solution for some other set of stellar fundamental parameters and composition. Before the iterative procedure starts, using the effective temperature of this starting structure, the temperature value at each of its depth points is rescaled by applying the ratio of this effective temperature to the required one (for more details see Appendix C).
Then the RE calculations start. At each of the iterations, the wavelength and depth dependent source function (consisting of the thermal and scattering parts) needs to be determined. This can be achieved by various methods. By default MPSATLAS uses a Feautrier method, but a second iterative method is also implemented (more details see Appendix D). Having found the source function, the moments of the intensity are calculated, namely the mean intensity J_{λ}, and the flux, H_{λ}. This process is repeated for each wavelength point, typically in an interval from 9 nm to 10 000 nm, in order to obtain the total wavelength integrated Eddington flux, ℋ.
The aim is to match the total Eddington flux, ℋ, to the flux that corresponds to the desired effective temperature. Considering the flux, ℋ, and keeping in mind that the majority of it comes from the region where 0.1 ≤ τ ≤ 2.0, the temperature corrections δT can be found in the optically thick regions. For that a modified version of the AvrettKrook procedure is applied to take deviations from the RE due to convection and overshooting into account (for more details see Avrett & Krook 1963; Kurucz 1970). In the optically thin regions, the temperature does not greatly affect the overall flux, whereas the derivative of it, dH_{ν}/dτ_{ν} = 0, is sensitive to a temperature change. The boundary condition for the temperature corrections δT in the optically thin part of the atmosphere requires the flux derivative of the Eddington flux to be zero (for a discussion see Mihalas 1978). This is achieved using a Λ correction method (BöhmVitense 1964). Then, the part of the atmosphere where both approaches overlap is smoothed out to match.
After the atmospheric structure is either calculated by the method described above or taken from some other source, for example from a ray through a 3D cube, and the source function is found, the emergent intensity can be calculated by evaluating the integral
The integral in Eq. (5) has to be computed for each wavelength separately to get the whole emergent spectrum. We note that the framework described here is restricted to a plane parallel setup. Thus, when using a computed 1D model, the emergent intensities for different view angles, μ, can be obtained by setting a set of μ angles in the code. To calculate spectra emerging from 3D cubes for different view angles, a 3D cube is rotated and for each view angle a different set of 1D rays is obtained for which the RT is solved along the ray, but keeping the view angle μ = 1.
2.3. The ODF approach
The synthesis of the line opacity is computationally demanding due to the tremendous number of spectral lines to be taken into account. For example, the default MPSATLAS line list contains more than 100 million atomic and molecular lines (Kurucz 2005b). These millions of lines lead to very complex and rich spectra for cool stars. Consequently, a very fine wavelength grid is needed to catch all the details in the spectra. For a number of applications, such detailed spectra are not required and lowresolution calculations are sufficient. However, lines still affect even lowresolution spectra. A straightforward way to include the effect of lines on lowresolution spectra is to simply average highresolution spectra. Let us consider a small wavelength interval (hereafter, bin), for example between 0.1 and 10 nm wide, between λ_{i} and λ_{i + 1}. This bin represents one lowresolution grid point, in which the intensity, I_{bin}, representing the whole bin, has to be calculated
The simplest way of approximating I_{bin} is by obtaining I_{λ} on a fine wavelength grid using a large number of points, N, and taking the sum
where Δλ_{j} is the jth discretised wavelength step. While this method is straightforward, the main issue is that I_{λ} has to be calculated N times for an accurate approximation, where N for example is of order O(10^{4}) for a one nanometre interval in the UV, if the resolving power R is 500 000. An alternative method to avoid solving the RT many times is the ODF approach (see e.g., Hubeny & Mihalas 2015, pp. 625–627), which approximates the complex structure of the opacity.
The main goal of the ODF method is to reduce the number of points N in a bin to a minimum and still approximate the flux accurately. However, the opacity in one bin can abruptly change by several orders of magnitude within very narrow spectral intervals (see Fig. 1a). Consequently, a fine spectral grid is required to accurately describe the opacity profile. One could think about averaging the opacity in the entire bin to avoid the highresolution calculations. In most cases such an averaging would lead to a gross overestimation of the line opacity effect, and would essentially lead to trapping photons that otherwise would escape (see Fig. 2 and its detailed discussion in Cernetic et al. 2019).
Fig. 1. Illustration of ODF generation in one example bin. (a) Detailed highresolution opacity in the bin from 420–422 nm. (b) Sorted opacity without the information of the wavelength, and the corresponding geometric mean values for 12 subbins, which were chosen as in Castelli (2005b). 
This can be circumvented by taking averages over wavelength points with similar opacity values, which can be achieved by grouping points in a certain opacity range together. A straightforward way is to sort the points in one bin in an ascending order by their opacity value. The sorted opacity profile can be described using significantly fewer points than that of the original opacity (compare Figs. 1a and 1b). The opacity profile in each bin is usually divided into several subbins, whereafter the opacity is averaged over each subbin using the geometric mean (see Fig. 1b). The geometric mean works better for the optically thick regime, because it avoids skewing the opacity towards large values. However, in the optically thin regime, the arithmetic mean is better, since the intensity depends on the opacity rather linearly. While both approaches are implemented, we used the geometric mean throughout this work in order to be consistent with the approach in the original version of the DFSYNTHE code.
Subbins might have different widths, which are defined by the fraction of the whole bin size. Generally, the part with the greatest opacity values is subdivided into smaller subbins, while the part with smaller opacity values is split into a few large subbins. The resulting stepfunction in the entire bin [λ_{i}, λ_{i + 1}], where each step is the averaged opacity κ_{s, i} of the sth subbin, is called opacity distribution function. For a more detailed description of ODF and the importance of different subbin configurations, that is the number of subbins and their widths, see Cernetic et al. (2019) and references therein.
Then, the intensity in each of the subbins can be calculated by solving the RT equation only once using the corresponding κ_{s, i} value. For the entire ith bin, the intensity is obtained by summing over the contributions from the subbins, that is
where I_{s} is calculated using κ_{s, i} and Δλ_{s} is the subbin width.
We note that by rearranging the opacity in a wavelength interval, as is done during the sorting in the ODF approach, two important assumptions are made implicitly: (i) the opacity shape in the bin does not change rapidly along the line of sight, in particular, within the region of the atmosphere where the radiation for the corresponding bin is formed (Kurucz et al. 1974) (ii) the wavelength interval of the bins are small enough that changes of the Planck function as well as changes of the continuum opacity can be neglected.
As mentioned in Sect. 2.1, the opacity for a given elemental composition and turbulent velocity, v_{turb}, only depends on pressure and temperature in the LTE case. Thus, the ODF table contains the information of the mean opacity in each subbin for a given elemental composition and microturbulence and can be pretabulated on a temperature and pressure (TP) grid (see e.g., Castelli 2005b, and references therein). For more details on the implementation see Appendix E. Having a pretabulated ODF table, it is now straightforward for any stellar atmosphere to obtain the total opacity for any temperature and pressure using linear interpolation on the TP grid and adding the continuous opacity. We note that when line opacity from the ODF tables is added to the continuous opacity, the assumption that the continuous opacity in a particular bin is independent of wavelength is made implicitly. This is currently a limiting factor for the upper limit of the bin size. One way to circumvent this limitation is to include the continuum opacity in the ODF tables.
3. Code structure
The MPSATLAS code results in one monolithic executable which depending on a control input file can perform a different set of calculations. An overview of the original three routines on which MPSATLAS is based is given in Fig. 2. We note that the original ODF generation routine^{1} consists of four separate codes, resulting in four executable. These three separate routines have been designed to perform the calculations outlined in Sect. 2 and correspond to the three MPSATLAS modules (see Figs. 3 and 4), where a module is a code internal execution mode. The wellestablished DFSYNTHE routine (Kurucz 2005a; Castelli 2005b) corresponds to module I in the MPSATLAS code and generates ODF tables (see left column in Fig. 2). Two slightly different ATLAS9 codes^{2} correspond to module II and module III in MPSATLAS, where one calculates 1D atmosphere models and the other the emergent intensity or flux (see middle and right columns in Fig. 2). The main difference between the two ATLAS9 codes is the main routine. The code for calculating model atmospheres contains routines to perform the temperature corrections in addition to solving the RT.
Fig. 2. Original structure of codes that calculate ODF, model atmospheres, and emergent intensities. Each green bubble indicates a separate code that results in a separate executable. The listed tasks are the main procedures used in each programme. Grey background indicates procedures that are the same for all three calculations. 
The original collection of codes is intricate to run due to several codes that need consecutive execution. Furthermore, the three different codes make use of similar procedures. For example, both the DFSYNTHE and the ATLAS9 codes require calculations of the populations and continuous opacity. While these calculations can be obtained using the same functions, both codes have their own copy of these functions with slightly different implementations. This renders any modification quite difficult. Therefore, we chose to merge the three codes into one code with different execution modes, that we call here modules, that take care of the different calculations. The resulting, merged code is further adjusted, parallelised, and simplified, and thus called MergedParallelisedSimplifiedATLAS. MPSATLAS is mainly written in Fortran 90 (free form), but some parts are left in F77. We will rework them in future releases. We use dynamical memory allocation in several modules, for example in the ODF calculations.
The advantages of the MPSATLAS code are a more user friendly operation, significant computational speedup, and wider functionality (in particular flexible ODF setup, which is discussed in Sect. 4). A schematic diagram in Fig. 3 shows the overall structure of the MPSATLAS code. The new code has three execution modes that correspond to the original structure of the three codes: ODF generation, model calculations, and emergent intensity or flux calculation. These three modules can be either executed consecutively or separately, which is controlled during runtime. Both model calculations and emergent intensity calculations need ODF tables as input, which is indicated by the arrow. As discussed in Sect. 2.2, the model calculation is not necessary if alternative input atmospheres are used instead, for example a ray from a 3D radiative MHD cube (such a bypass of module II is indicated by the arrow in the bottom of Fig. 3).
Fig. 3. Schematic structure of the code modules I–III. 
In the following we give an overview of the code structure in terms of input and output. There are two types of parameters (shown in Fig. 4): (i) overall input parameters (highlighted in green), that are specified once and do not change for a particular calculation, for example the elemental composition (ii) parameters that are set in the input, but which determine a grid on which the output calculations are performed (highlighted in grey). Figure 4 shows that the calculations of an ODF table using module I require the elemental composition, line lists, and the turbulent velocity. In addition, the temperature and pressure grid, as well as frequency range, and the bin, subbin configuration need to be specified (more details in Appendix E). An example input file for the ODF calculation is shown in Appendix G. The output ODF is needed as input for both module II, and module III.
Fig. 4. Input parameters for different code modules together with information on the output. 
Module II calculates a 1D atmosphere model in RE. As input parameters this module needs the elemental composition, effective temperature (T_{eff}), surface gravity (log g), and the mixinglength parameter, which is needed to include convection through mixing length theory (BöhmVitense 1958), for which overshoot can be turned on. In addition, an ODF table is required as input to account for the effect of line blanketing on the atmosphere structure. For the opacity table, either the output from module I, or from a precalculated ODF table grid can be used. Moreover, an initial 1D model atmosphere is read, which is used as an initial starting point, and has the same format as the output model (for more details see Appendix C). The output model is obtained on an a priori specified τ_{Ross}−grid. An atmosphere model consist of column mass, temperature, total gas pressure, electron number density (n_{e}), Rosseland mean opacity (κ_{Ross}), radiation pressure, and microturbulence velocity (v_{turb}) for each depth point. Currently, we only consider cases of height independent v_{turb}, though for heightdependent microturbulence it is possible to use several ODF tables covering a range of v_{turb} values and interpolate between them.
The third module either synthesises the emergent flux or the emergent intensity for different view angles. Here, the only actual free parameters to set are the wavelength range and viewing angles, as the input is already predetermined by the ODF and model atmosphere. While the elemental composition should match the ODF, the wavelength grid on which the ODF is calculated sets the wavelength grid of the emergent intensities. Moreover, the ATLAS output 1D atmosphere model in RE contains all the necessary information to calculate the emergent spectrum, in particular the electron number density, and the radiative pressure at each depth point obtained with the elemental composition for which the model was calculated. On the contrary, if using an external model for example from a ray through a 3D MHD cube, several quantities might be unknown. In particular, the electron number densities, which are needed to find the populations of all other ions, have to be obtained. In order to use external models that only have the column mass, pressure, and temperature, an additional flag was introduced to recalculate the equilibrium number densities in each depth point for the given elemental composition. Example input files are discussed in more detail in Appendix G.
4. Functionality extensions
In the following we describe the main functionality extensions included in the MPSATLAS code.
4.1. Flexible wavelength grid and optimised subbinning
The DFSYNTHE code is limited to two particular wavelengthbin grids onto which ODFs are calculated, both with the same wavelength independent subbin sizes. However, various applications require spectral synthesis with different spectral resolutions. Hence, it is important to have an option for changing the ODF wavelength grid. Furthermore, Cernetic et al. (2019) showed that an optimal choice of the subbin grid leads to significant improvements in both accuracy and speed of the calculations. In particular, they showed that for most of the wavelength range two to four subbins are sufficient to reach the same accuracy as the standard 12 subbins (which are hard coded into DFSYNTHE). Such a threefold speedup is especially important for 1.5D calculations which require an immense number of 1D calculations (about a million per cube with 1000 times 1000 horizontal grid points).
In this context, we implemented a flexible treatment of the bin grid, enabling the user to choose the overall wavelength interval on which the opacity is calculated (keeping the maximum range between 9 nm and 10 000 nm), as well as the binning and subbinning on this interval (more details on the implementation can be found in Appendix G). The underlying highresolution opacity is then only calculated in the chosen wavelength interval to save computational time. The flexible wavelength binning is in particularly useful if synthesised spectra need to be compared to observations of different resolving power. In order to reduce the computational cost, the number of subbins per bin can be reduced significantly if an optimal subbin border distribution is chosen. In addition, having a flexible binning and subbinning grid allows rapid computation of ODF tables for broadband filters, if only the flux through a particular filter is needed (Cernetic et al. 2019; Anusha et al. 2021).
Moreover, we have added an option to pretabulate the opacity on a highresolution grid. For that, the code reads the starting and end wavelength of an interval, and uses a resolving power of R = 500 000 in this interval. The opacity is then written in the same format as an ODF, but using only one subbin, which contains the opacity of the wavelength point. This can be used to calculate emergent intensity at high spectral resolution. Since the line opacity on a highresolution grid might have strong changes with temperature and pressure, the interpolation error is potentially greater compared to an ODF. An exemplary calculation of the Vega flux in the wavelength region of the Balmer jump shows that the error between the flux using the interpolated opacity (using a TP grid that splits the same ranges as used in Castelli 2005b into 100 logarithmically equidistant steps) and the calculated highresolution opacity can reach up to 3% (see Fig. 5). The largest deviations occur in line wings of strong lines (see Fig. 6), which indicates that the temperature steps in the considered TP grid are insufficient to account for the line broadening sensitivity. Therefore, the TP grid can be adjusted depending on the needed accuracy (for more details see Sect. 4.2).
Fig. 5. Highresolution emergent flux for Vega around the Balmer jump region, in the range of 360 nm to 400 nm. (a) Flux calculated using interpolated opacities from pretabulated highresolution opacities on a TP grid and using opacities calculated for each depth point. (b) Ratio of the fluxes from (a). 
Fig. 6. Same as Fig. 5, but for a smaller wavelength interval (395.5 nm to 398.5 nm), corresponding to the Hε line core. 
4.2. Flexible TP grid for pretabulation
The DFSYNTHE code uses a predefined TP grid of 57 temperature and 25 pressure values for pretabulating ODF tables. The pressure covers 12 orders of magnitude ranging from 10^{−4} to 10^{8} and the temperature ranges from 2 × 10^{3} − 2 × 10^{6} K. In MPSATLAS the number of TP points and their values can be specified (as input parameters). This has two advantages. On the one hand, the numerical cost can be reduced, while the resolution of the TP grid is kept, but the ranges of temperatures and pressures are adjusted to certain types of stars. On the other hand, the resolution can be increased if more accurate interpolation is needed. This is especially important if the code is used to generate pretabulated highresolution opacity (see Fig. 5).
4.3. Line preselection criterion
The most extensive line lists such as Kurucz’s original line list^{3} and VALD3 (Ryabchikova et al. 2015) contain many lines that do not contribute to the opacity within the temperature and pressure range of the ODF table. To reduce the computational cost in the line opacity calculations, the lines are preselected based on their line core opacity and the continuous opacity for each TP point. For that, lines which lead to an opacity in the line core of several orders of magnitude (controlled by the cutoff factor) less than the continuous opacity are discarded, see Appendix E for the selection criteria. By default the cutoff factor is set to 10^{−3}. Here, we tested how the solar surface intensity is affected by choosing several smaller cutoff factors (see Fig. 7). Since a smaller cutoff factor increases the computational time significantly, the smallest value we tested is 10^{−12}. The difference between a cutoff factor of 10^{−9} and 10^{−12} is negligible. However, overall the largest decrease in intensity can be observed around 200 nm, where it reaches up to 5%. This implies that even very weak lines have a nonnegligible effect on the flux, in particular if there are many of them.
Fig. 7. Ratios of the solar emergent intensity at the disccentre calculated using ODF tables that were generated with different cutoff factors for the line preselection. Here cutoff factor refers to the opacity of a given line relative to the continuum opacity. If this ratio drops below the cutoff factor, this particular line is not included in the computation of the intensity spectrum. The default value of the cutoff factor is 10^{−3}. 
4.4. Changing line lists
Kurucz’s original line list includes an immense number of theoretically computed lines which have not been measured experimentally. Over the last few decades, the atomic and molecular data has been constantly updated to obtain more accurate data. Recent developments led to various line lists which are updated regularly. In particular, there have been several important updates of the Vienna Atomic Line Database (VALD), which is a database focusing on atomic data relevant for modelling of stellar atmospheres (Ryabchikova et al. 2015).
While the DFSYNTHE code only works with Kurucz’s line list, MPSATLAS can be used with different line lists. For the current tests we exchanged the atomic line list, but we kept the line list for diatomic and H2O lines (Schwenke 1998) as in the original DFSYNTHE version.
To illustrate the effect on the flux caused by the differences in the line lists, we first calculated the model atmosphere for the Sun using the standard ‘little’ grid ODF (Castelli 2005b) with the original line list and the VALD3 line list. Then, using these models, we calculated the intensity at disccentre on a highresolution wavelength grid with resolving power R = 500 000 for each of the two line lists. Figure 8 shows both intensities in the range between 5340 Å–5360 Å together with the Hamburg atlases of the solar spectrum^{4} (Neckel 1999; Doerr et al. 2016). For this comparison, we degraded the resolving power of the calculated flux using a Gaussian kernel. It becomes evident that the intensities obtained using the VALD3 line list agree significantly better with the solar observations, while the intensities obtained using Kurucz’s original line list show quite a few lines that are not present in the solar spectrum. This is because Kurucz’s line list contains significantly more lines than VALD3 and most of these excess lines have never been measured in the laboratory.
Fig. 8. Highresolution solar disccentre intensity in the range 5340 Å–5360 Å computed with MPSATLAS using different line lists together with data from the Hamburg atlas of the solar spectrum (Neckel 1999; Doerr et al. 2016): (a) computed intensity using Kurucz’s line list and the Hamburg atlas data, (b) intensity using the VALD3 line list and the Hamburg atlas data. 
In a second step, we compare the overall energy distribution for the Sun and a Ktype star (T_{eff} = 4000 K) obtained using different line lists. For that, we calculated the model atmospheres using the little grid ODF with different line lists, and subsequently generated the overall disc integrated emergent flux using the standard little grid ODF with the corresponding line list. Figure 9 shows that the fluxes calculated with the original Kurucz’s line list are significantly different from the fluxes obtained using the VALD3 line lists below 450 nm. The effect in the UV is huge for the Sun, where the flux obtained using the original Kurucz’s line list gives a better agreement (see a more detailed discussion in Sect. 5.2). However, for a Ktype star for which line opacity is even more important than for the Sun, there is rather a moderate difference. This indicates that the main difference between the VALD3 and the Kurucz’s line lists are lines whose lower level have higher energies (i.e. mainly lines contributing to opacity at higher temperatures than in Kstars).
Fig. 9. Ratios of the emergent intensity at disccentre calculated using ODF tables with the VALD3 line list to the emergent intensity at disccentre calculated using ODF tables with Kurucz’s original line list. 
4.5. NH photodissociation
Comprehensive comparison of observed and modelled solar spectra showed that molecular photodissociation is an important source of continuum opacity in the UV range (Fontenla et al. 2011). Moreover, Fontenla et al. (2015) proposed that NH photodissociation is the source of the missing opacity in the UV (see, e.g., Busá et al. 2001; Short & Hauschildt 2009; Shapiro et al. 2010, for the detailed discussion of the missing opacity problem). So far only the molecules CH and OH were included in the ATLAS9 continuous opacity calculations (Castelli 2005a). We added the opacity contribution from the NH molecule using crosssections that we obtained from Kurucz (2020, priv. comm.). A more detailed description on how molecules are included in MPSATLAS is given in Appendix A, and how CH, OH, and NH photodissociation opacities are implemented in Appendix F. This additional opacity source changes the specific emergent intensity by at most 0.5% (see Fig. F.1). We conclude that NH is not sufficient to explain the missing opacity in the UV.
4.6. MPI parallelisation for faster computations
While the RT calculations can be sped up by using optimised ODF configurations, the calculations of the ODFs table itself remain timeconsuming. To achieve faster calculations, we parallelised the module for ODF generation (module I in Fig. 3) along temperature (T) values. We note that the number of lines contributing noticeably to opacity depends on temperature. In particular, more lines have to be calculated for lower temperature values so that producing ODFs for them is the most time consuming. To account for this imbalance of computational time, we used a masterslave implementation. Here, one core distributes the T values for which the ODF must be calculated to all other cores as soon as they finish with their previous values. Consequently, while different cores calculate ODFs for different numbers of temperature values, the number of calculations is distributed between the cores more or less equally. This also has the advantage that the number of T values does not need to be divisible by the number of cores used. With this implementation, an example ODF table on the same TP grid as used in the original DFSYNTHE code can be calculated in under 10 min on 36 cores on a modern high performance cluster.
For the RE calculations (module II in Fig. 3), the integrated flux over the whole spectrum is needed, where the number of frequency points can be varied. Calculations with a large number of frequency points, or if the highresolution grid is used, are very time consuming as the RT needs to be solved in each iteration on chosen frequency grid. Thus, we also parallelised the RT calculations along wavelengths in the module for modelling atmosphere structures.
Finally, module III that calculates the emergent intensity (see module III in Fig. 3) uses a parallelisation if the intensity for more than one model has to be calculated. This is only the case for postprocessing 3D calculations. In this case the atmosphere models are distributed among the available cores.
4.7. Speedup and portability
Using the parallelised version of the ODF module, we tested the MPSATLAS scalability. For the ODF calculation, we used the standard TP grid as in Castelli (2005b), but set the wavelength interval from 200 nm–10 000 nm, which significantly reduces the computational time. The low temperature values in the standard grid (below 5000 K) require the highest number of lines, and thus are very timeconsuming compared to larger temperature values.
Figure 10 shows the speedup for ODF generation and RE calculations. It becomes evident that there is no further speedup of the ODF calculations for more than 10 cores. This is because the lowest temperature value (1995 K) forms the bottleneck in this example. Even if all other temperature values are calculated faster using more cores, at the end of the entire ODF calculation the temperature value that takes the longest determines the amount of time needed. The number of cores at which such a bottleneck (due to one particular temperature) appears depends on the TP grid and wavelength range of the ODF table. In general, the fewer P points and more T points a grid has, the higher the number of cores until the speedup saturates.
Fig. 10. Top panel: speedup for ODF table generation as a function of the number of compute cores used (module I from Fig. 3). Bottom panel: speedup for atmosphere model calculation (module II from Fig. 3). Black lines in both panels indicates the speedup averaged over 5 runs. The grey shaded area indicates the spread in the 5 runs, and the red line indicates the ideal speedup. 
To test the parallelisation for the RE calculations, we chose to model a 4000 K star with log g = 4.8, and M/H = 0.3. We used a standard little grid ODF (Castelli 2005b), and reduced the wavelength interval to the first 1210 points on the grid, which corresponds to 9 nm–10 000 nm. The model calculation converges after 37 iterations. The speedup of the RE calculations (module II) shows a different behaviour: the increase in speedup starts to flatten after 8 cores (see Fig. 10b). The flattening happens because only the part of the RT calculations along the wavelength points is parallelised, while all other calculations, especially the temperature corrections are executed on one core. Here, the behaviour of the speedup depends on the number of wavelength points that need to be calculated.
The parallelisation for the emergent intensity calculations (module III) is along atmosphere models, thus it is only useful for calculating along model atmospheres from 3D cubes. Since it only needs a few communications in the beginning to distribute the input settings, we expect an almost ideal speedup.
The code is fully portable. The available git version of the code includes an installation script that downloads all necessary libraries and files, such as a NetCDF library, line lists, and example input files. The code can be compiled using an openMPI and intel compiler, as well as a corresponding gnu compiler. Note that for the gnu version, a different line list format is needed, which we provide on request.
5. Benchmark
In order to put the MPSATLAS code into context, we compared emergent flux obtained by the MPSATLAS code to that from other codes and to observations. For the codetocode comparison, the stellar parameters were chosen to be available on both the PHOENIX and Kurucz model grid of calculated specific intensities. To this end, we considered three hypothetical stars of spectral classes A, F, and K (with corresponding temperatures of 8000 K, 6500 K, and 4000 K). These spectral classes have been chosen to test the performance of the code for cases when opacity is dominated by different sources. For example, in a Ktype star the opacity is mainly due to millions of atomic and molecular lines, while the main opacity source in an Atype star is the continuum and hydrogen lines. Furthermore, for comparison to observations, we used the Sun and Vega (see Table 1 for the summary of the fundamental stellar parameters of the stars we used). We chose Vega for comparison, because accurate measurements of the total flux over a large wavelength range are available.
Stellar parameters for the comparison stars.
5.1. Codetocode comparison
We compare emergent fluxes calculated using the MPSATLAS code to the original fluxes from the ATLAS9 grid calculated by Kurucz^{5} (Kurucz 1993), and to the modelled fluxes on the PHOENIX grid ^{6} for three different stellar types (Ktype, Ftype and Atype; see Table 1 for details). In our calculations we assume the same setting as in the Kurucz grid modelled fluxes: convection is turned on, without overshoot, and the mixing length is set to 1.25. We further consider a constant microturbulence of 2 km s^{−1}. Generally, the PHOENIX code provides a more intricate setup, for example the models are calculated in spherical geometry, condensation is included in the equation of state, and some nonLTE effects can be turned on. However, the model grid we use (Husser et al. 2013) was obtained using LTE and the RT calculations only make use of the nonLTE effects by a special line profile treatment for some species (Li I, Na I, K I, Ca I, Ca II). Moreover, the mixing length and the microturbulence parameters vary on the PHOENIX model grid, and thus can slightly differ from the one used in ATLAS9 and MPSATLAS (for more details see Husser et al. 2013).
We preformed calculations using two different elemental compositions. For a comparable calculation to the Kurucz grid, we used the same element abundances as in Kurucz’s calculations. These abundances are taken from Anders & Grevesse (1989) and will be referred to as ‘Anders composition’. For the comparison with PHOENIX, we used the more uptodate ‘Asplund composition’ taken from Asplund et al. (2009). For the modeltomodel comparisons, both the model atmosphere and the emergent spectra are calculated using the standard little ODF with 1221 frequency points (Castelli 2005b). We smoothed all spectral fluxes by applying the average over a 15 nm interval around each wavelength grid point. The PHOENIX highresolution spectral fluxes are first averaged over the same intervals as given by the ATLAS wavelength grid, and then smoothed with the same procedure over 15 nm intervals.
For an Atype star the emergent spectral flux returned by the MPSATLAS code and those from PHOENIX and ATLAS9 grids are shown in Fig. 11a. For a more detailed comparison we show the deviations of the calculated MPSATLAS fluxes and the PHOENIX to the original ATLAS9 fluxes in Fig 11b. Overall there is a reasonable agreement between the three codes. The largest deviations are in the UV for wavelength shorter than 400 nm. A good example for the importance of the elemental composition can be seen in the range between 210 nm and 290 nm. While the flux obtained by using the MPSATLAS code and the Anders composition is closer to the ATLAS9 calculations (also performed with the Anders composition), the MPSATLAS flux for the Asplund composition is closer to the PHOENIX calculations (performed with the Asplund composition). Large differences can occur especially in the UV due to different line lists used. In the wavelength interval between 300 nm–1000 nm, the PHOENIX calculation gives an overall smaller spectral flux than fluxes calculated with MPSATLAS and ATLAS9. The difference is approximately 5% around 400 nm and then decreases to a few percent for wavelengths greater than 1000 nm. In this wavelength interval the dominant continuum opacity contributors are H^{−} boundfree transitions. Thus, either a slight difference in the equilibrium number densities or a different implementation of the H^{−} boundfree crosssections can potentially lead to these differences.
Fig. 11. Flux comparison between MPSATLAS, KuruczATLAS9, and PHOENIX code for an Atype star with the effective temperature of 8000 K. Flux values are shown (top panel) together with the corresponding flux deviations (bottom panel) in % compared to the original Kurucz calculations. 
In addition, there are three larger deviations around 820 nm, 1458 nm, and 2279 nm which correspond to the Paschen limit, the Brackett limit and the Pfund limit, respectively. This implies that all three codes have a different treatment of the hydrogen continuous transitions of these series. On the contrary, the Balmer series transitions, Hα (at 656 nm) to Hζ (at 389 nm), show only very slight differences between the three codes. The significant deviations around 395 nm are caused by the Ca II H and K doublet, where the MPSATLAS calculation gives more similar results to PHOENIX than to the older ATLAS9 calculations. We note, however, that the proper treatment of the Ca II H and K doublet requires nonLTE modelling which is absent in all three codes.
In Fig. 12 the emergent flux for a Ftype star (see Table 1 for exact fundamental parameters) is displayed. Overall the deviations in the UV are greater and reach up to 20%, while the differences in the visible and infrared (IR) are smaller. The two elemental compositions entering the MPSATLAS lead to the following differences in the spectra in the range 200 nm–500 nm: the flux is larger for the Asplund composition, and it decreases towards 450 nm whereafter the flux is smaller for the Asplund composition than for the Anders composition. This behaviour is most probably caused, directly or indirectly, by the effect of the composition on line blanketing. The line opacity, in particular from iron, is smaller for the Asplund composition and allows more photons to escape. Then, the decrease in the visible can be explained by a compensation effect as the wavelength integrated flux has to be the same for both element compositions. The PHOENIX flux oscillates around the MPSATLAS flux obtained using the Anders composition in the region 200 nm–400 nm, which implies that the line opacity differs significantly due to a different line list. Finally, the wavelike deviation between 500 nm–2600 nm indicates that the photospheric structures differ somewhat from each other. It becomes evident that the PHOENIX, and MPSATLAS model have a similar structure, while the older Kurucz’s model deviates.
The deviations between the codes are larger for a Ktype star (see Fig. 13). In the UV the difference between PHOENIX and other codes reaches 160%, this might be due to a smaller number of UV lines contributing to opacity at lower temperatures (i.e. below 4000 K) in the line list used in the PHOENIX code. The flux calculated using MPSATLAS deviates by up to 25% from Kurucz’s originally computed flux. Focusing on wavelengths longward of 400 nm, the deviations oscillate, but the amplitude of the deviations decreases towards the infrared. Comparing the difference between the Anders composition and the Asplund composition, the flux obtained using the Asplund composition is closer to the PHOENIX calculation for most of the wavelength in the visible. This is also the case in the infrared, where the overall deviations become less than 10% between all four models.
Fig. 13. Flux comparison between MPSATLAS, KuruczATLAS9, and PHOENIX code for a Ktype star with T_{eff} = 4000 K. Flux values are shown (the two top panels) together with the corresponding flux deviations (in the two bottom panels) in % compared to the original Kurucz calculations. 
5.2. Codetoobservation comparison
Besides comparing with the output of other codes, it is important to also check how well MPSATLAS reproduces observations. For the Sun a lot of accurate highresolution data exist, as well as lowresolution spectra. In contrast, for most of the stars either broadband fluxes, intermediate resolution spectra or normalised (e.g., to continuum level) highresolution spectra for a rather small wavelength interval are available. Here we first test how well the solar flux can be matched by the MPSATLAS code, and in a second step we chose Vega as a comparison star.
5.2.1. The Sun
For the calculations presented in this subsection, both the model atmosphere and the emergent spectra are calculated using the standard little ODF with 1221 frequency bins (Castelli 2005b). For all considered cases, we recalculated the model when changing any parameters. We always assumed a microturbulence of 1.5 km s^{−1} and if convection was turned on, the mixinglength is 1.25. We compare our calculations to Solar Irradiance Reference Spectra (SIRS) for the 2008 whole heliosphere interval (WHI; Woods et al. 2009). For a better comparison, all calculated fluxes are smoothed out using a trapezoidal kernel (Harder,priv. comm.) to match the WHI resolution in wavelength regions where the little ODF resolution is higher than that of the WHI. We kept the original resolution otherwise.
Figure 14 shows the solar irradiance at a given wavelength at one AU from the Sun as in the WHI compared to MPSATLAS calculations using the Anders composition with and without convection, and overshoot. One can see that in the mid UV the deviations between the calculated and observed irradiance significantly exceed the uncertainty of the measurements. This is not surprising (and similar deviations are also present in the ATLAS9 and PHOENIX calculations, see below) since it is well known that currently available line lists miss a significant number of weak lines (this excess flux constitutes the famous ‘missing UV opacity’ problem, see e.g., the recent review by Rutten 2019, and citations therein). We note that the solar UV spectrum is also affected by deviations from LTE, which are ignored in MPSATLAS calculations. However, it has been shown that accounting for nonLTE effects will make the deviations between the calculated and observed UV flux level even larger (see, e.g., Short & Hauschildt 2009; Shapiro et al. 2010; Tagirov et al. 2019).
Fig. 14. Irradiance calculated for solar parameters compared to WHI. The light and dark grey shaded area indicates the minimum and maximum measurement error of the SIRS WHI (Solar Irradiance Reference Spectra for the 2008 Whole Heliosphere Interval), respectively (see, Woods et al. 2009, for a detailed description). Solar irradiance (top panels) calculated using MPSATLAS with three different assumptions: no convection, with convection but no overshoot, and with convection and overshoot. The flux deviations between WHI observed irradiance and models in % are shown in the bottom panels. 
In contrast to the UV spectral domain, the modelled irradiance agrees very well with the observations in the visible and infrared. Figure 14 shows that the deviations there are mostly within the minimum measurement uncertainty. While the differences between the three calculations are very small, the agreement between the observations and the MPSATLAS calculation including convection with overshoot is best in the interval 450 nm–650 nm. However, in the region 650 nm–900 nm the fluxes including overshoot show the greatest deviations. The wavelike behaviour of the deviations in the range between 500 nm–1500 nm might indicate that the modelled atmosphere temperature in the region where the continuum is formed does not accurately match the actual temperature in the Sun.
The solar elemental composition has been intensively investigated for the last decades and updated following more accurate modelling (Anders & Grevesse 1989; Grevesse & Sauval 1998; Asplund et al. 2009). Changing ratios between hydrogen, helium, and heavier elements, which significantly contribute to the electron number density, affects the structure of the atmosphere and the spectral synthesis in a competing way. Increasing the concentration of the electron donors without recalculating the atmospheric structure leads to a drop of the flux, and thus a decrease of the effective temperature. However, as soon as the structure is recalculated in RE, it compensates for the decreased flux and the effective temperature should return to its original value.
In the next step, we show the difference in the irradiance for the Anders composition and Asplund composition, both with convection but no overshoot. Figure 15 shows flux model calculations for the different elemental compositions, together with the WHI measurements. The difference in the elemental composition leads to a redistribution of the flux in different wavelength intervals due to a change of equilibrium number densities of the species (which in turn affects the opacity). The effective temperatures obtained from the total flux are almost identical for the two compositions, being T_{eff} = 5778.9 K for the Asplund composition, and T_{eff} = 5778.5 K for the Anders composition. Nevertheless, the irradiance calculated using the Asplund composition has a higher flux in the UV, compared to both the WHI and the flux using the Anders composition, so that the former displays greater deviations from the observations. In contrast, the flux calculated using the Asplund composition is lower in the visible and IR compared to that calculated using the Anders composition. We note that the same behaviour was also observed for the Atype and Ftype star of solar metallicity (see discussion in Sect. 5.1 and Figs. 11 and 12). Both calculations match the observations within the measurement uncertainties for wavelengths greater than 450 nm.
Fig. 15. Irradiance calculated for solar parameters compared to WHI. Same as in Fig. 14, except that the MPSATLAS model calculations were done without overshoot and using the Anders and Asplund elemental compositions. 
Finally, in Fig. 16 we plot the flux obtained by MPSATLAS for the Anders composition without overshoot, together with the original fluxes calculated by Kurucz^{7}, and the PHOENIX flux calculations. Since the PHOENIX grid does not have the flux for the solar effective temperature, and surface gravity, we used linear interpolation between the closest available stellar parameters. Furthermore, before applying the same smoothing procedure using a trapezoidal kernel to match the WHI resolution, we averaged the PHOENIX flux over the same wavelength intervals as given in the standard little ODF grid. The interpolated PHOENIX flux shows a slightly worse agreement with the observations in the UV, but an overall good agreement in the visible and IR. Kurucz’s original calculation agrees better with the observations in the interval 450 nm–650 nm than the PHOENIX flux. All three codes show a wavelike behaviour between 500 nm–1500 nm indicating a slight mismatch of the modelled atmosphere as discussed above.
Fig. 16. Irradiance calculated for solar parameters compared to WHI. Same as in Fig. 14, but without the Asplund composition, while in addition the solar flux provided on Kurucz’s website and the flux from PHOENIX are shown. 
Note that we used Kurucz’s original line list for the above atmosphere model and flux calculations. The modelled flux in the UV is significantly larger than the observed values, and using the VALD3 line list leads to an even greater flux in the UV (see Fig. 9). This indicates that Kurucz’s line list leads to a better agreement of the calculated UV flux to observations. This is because Kurucz’s line list contains significantly more lines than VALD3 and even though most of these lines have never been measured in the laboratory they allow a better representation of the UV line haze than the VALD3 lines.
5.2.2. Vega
The other star for which accurate absolute spectrophotometry is available is Vega. The most uptodate observed absolute flux consists of data from two instruments, the International Ultraviolet Explorer (IUE), and the Space Telescope Imaging Spectrograph (STIS) mounted on the Hubble Space Telescope (Bohlin & Gilliland 2004; Bohlin et al. 2014). The fundamental parameters for Vega still contain a non negligible uncertainty (Castelli & Kurucz 1994). Therefore, we tested a small range of stellar parameters to seek good agreement, where the range was chosen based on previous stellar parameter determinations (Castelli & Kurucz 1994; Catanzaro et al. 2014). Namely, for the model comparison we calculated a small set of models with slightly different surface gravity values (log g = [3.90, 3.95, 4.00]), and different metallicity (M/H = [ − 0.3, −0.5]), but the same effective temperature, T_{eff} = 9550 K. The Anders composition was used for the abundances.
For the search of the best stellar parameters, we used the flux calculations obtained using the standard little ODF with 1221 frequency bins (Castelli 2005b). The spectral flux observed for Vega was averaged on the same wavelength grid, and subsequently we broadened the fluxes using a Gaussian kernel with a full width at half maximum (FWHM) of 3 nm below 350 nm, and with a FWHM of 5 nm above 350 nm. The output of MPSATLAS is the flux at top of the stellar atmosphere. To connect it to the observed flux, one needs a scaling factor, sf = (d_{Vega}/R_{Vega})^{2}, where R_{Vega} is Vega’s radius and d_{Vega} is the distance of Vega to the observer. This factor is still uncertain and some assumptions must be made for the comparison (see Castelli & Kurucz 1994, where is was estimated to be sf = (1.62 ± 0.07)×10^{16}).
Here, we take this factor into account by taking the ratio of the observed flux to the modelled flux in the wavelength region 400–600 nm, where the ratio showed the weakest dependence on the metallicity and surface gravity values. The scaling factor we obtained is sf = 1.59 ⋅ 10^{16}, which is within the uncertainty of the estimation. Overall, the best agreement of the modelled flux to the measured flux was achieved using M/H = −0.5 and log g = 3.90, but very little difference is found between the flux for log g = 3.90 and the flux for larger surface gravity of log g = 4.00.
For the comparison to PHOENIX, we downloaded the PHOENIX emergent intensities for T_{eff} = 9400 K, and T_{eff} = 9600 K, and the surface gravity values log g = 3.50 and log g = 4.00, with the metallicity M/H = −0.5. After calculating the fluxes, we had to interpolate to get to the effective temperature of T_{eff} = 9550 K, and a surface gravity of log g = 3.90. Due to the limited wavelength range for which the PHOENIX emergent intensities are provided, we can not calculate the exact effective temperature from the flux.
The PHOENIX flux is averaged in exactly the same way as the IUE and STIS data, and subsequently all fluxes are broadened using a Gaussian kernel. Figure 17 shows the comparison in the UV (120 nm–350 nm) between the measured flux, the flux calculated using the MPSATLAS code, and the PHOENIX flux. The MPSATLAS calculation has an overall better agreement with the measurements than the PHOENIX modelled flux, but still deviates significantly for wavelength shorter than 260 nm. Note that the elemental composition in the PHOENIX flux is the Asplund composition by default, but scaled to M/H = −0.5.
Fig. 17. Comparison of measured and calculated UV flux for Vega. The grey area indicates systematic and statistical measurement uncertainties. For comparison, MPSATLAS flux and PHOENIX flux are plotted. The modelled fluxes are scaled with the same factor and were obtained for the same effective temperature, metallicity, and surface gravity (Table 1). Top panel: absolute flux, with error. Bottom panel: deviations of modelled fluxes from measured flux in %. 
The comparison in the range 350 nm–1000 nm is displayed in Fig. 18. The MPSATLAS calculations show very good agreement, with only a few spectral intervals with greater deviation than the measurement uncertainty. In contrast, there is an offset in the PHOENIX flux, whose deviation from the observations shows a significant wavelength dependence that continues into the UV range (see Fig. 17). While an overall shift could be explained by the scaling factor (which was taken to be the same as for MPSATLAS), the wavelength dependence cannot be removed by a different normalisation. This indicates that the atmospheric structure or the continuum opacity has a slightly worse agreement. Moreover, the Paschen limit (820.4 nm) deviates significantly. The agreement of the modelled flux using MPSATLAS with the measurement is very good in the visible and nearinfrared.
Fig. 18. Comparison of measured and calculated flux for Vega. The observed flux consists of data from STIS. The grey area indicates systematic and statistical measurement uncertainties. For comparison, MPSATLAS flux and PHOENIX flux is plotted. The modelled fluxes are scaled with the same factor and were obtained for the same effective temperature, metallicity, and surface gravity. Top panel: absolute flux, with error. Bottom panel: deviations of modelled fluxes from measured flux in %. 
Finally, we compared calculated and observed highresolution spectra around the hydrogen Balmer and Paschen lines. For that, we used the same atmospheric structure as before, but we calculated the highresolution flux with resolving power R = 500 000. Subsequently, we used Gaussian broadening to degrade the resolving power to R = 400. In Fig. 19 the observed and calculated flux together with the deviation of the calculation to the observations are shown. The H_{α} line shows the best agreement with at most 4% deviations, while the deviations for the Paschen lines are up to 6% and for higher Balmer lines even up to about 10%. We note that in a previous comparison (Bohlin & Gilliland 2004) a better agreement (with deviations of a few percent) was achieved. However, these differences are mainly attributed to differences in the line depth of the updated STIS data (see Fig. 4 in Bohlin & Gilliland 2004), while some small deviations might be a result of modifications and updates in the MPSATLAS code compared to the original ATLAS9 at that time.
Fig. 19. Comparison of STIS flux with MPSATLAS calculation with R = 400. (a) Hydrogen Balmer series. (b) H_{α} line. (c) Paschen series. (d)–(f) the deviation of the calculated flux to the measured flux in %. 
6. Summary and conclusion
We presented the structure and extended functionality of the MPSATLAS code. The code is based on the ATLAS9 code, but has been extended to allow flexible and faster handling of ODF calculations, as well as emergent flux calculations. This makes the code suitable for radiative transfer calculations along rays from 3D MHD cubes. Furthermore, the atmosphere model calculations were sped up and made more userfriendly. We also improved the equilibrium number density calculations, and included NH photo dissociation opacity. The codetoobservations comparison showed that MPSATLAS gives excellent agreement with the observed solar spectrum and Vega spectrum (better than some other widely used codes).
The source code is available on request along with a set of testing input files. A detailed explanation of the input files is given in Appendix G. Furthermore, an online tool for calculating ODFs, stellar models, and fluxes will soon be released.
A fine grid of stellar models, fluxes, and centretolimb variations is being calculated (Kostogryz et al., in prep.). This grid will cover the range of effective temperature between 3500 K and 9000 K in 100 K steps, a range of surface gravity from log g = 3.0 to log g = 5.0, and metallicities between –5.0 to 1.5 using very fine steps around solar metallicity of 0.05 dex.
Acknowledgments
This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 715947). This work has been partially supported by the BK21 plus programme through the National Research Foundation (NRF) funded by the Ministry of Education of Korea and by the Max Planck Society grant “PLATO Science” and DLR PLATO grant Nr. 50OO1501 and 50OP1902. YCU acknowledges funding through STFC consolidated grants ST/S000372/1. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna.
References
 Ai, Y. L., Wu, X.B., Yang, J., et al. 2016, AJ, 151, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Anusha, L. S., Shapiro, A. I., Witzke, V., et al. 2021, ApJS, 255, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., & Jacques Sauval, A. 2006, Nucl. Phys., 777, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Avrett, E. H., & Krook, M. 1963, ApJ, 137, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Beeck, B., Schüssler, M., Cameron, R. H., & Reiners, A. 2015a, A&A, 581, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beeck, B., Schüssler, M., Cameron, R. H., & Reiners, A. 2015b, A&A, 581, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109 [Google Scholar]
 Bohlin, R. C., & Gilliland, R. L. 2004, AJ, 127, 3508 [NASA ADS] [CrossRef] [Google Scholar]
 Bohlin, R. C., Gordon, K. D., & Tremblay, P. E. 2014, PASP, 126, 711 [NASA ADS] [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 BöhmVitense, E. 1964, SAO Spec. Rep., 167, 99 [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Busá, I., Andretta, V., Gomez, M. T., & Terranegra, L. 2001, A&A, 373, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carbon, D. F. 1984, ed. Kalkofen, W. Methods in Radiative Transfer (Cambridge University Press), 395 [Google Scholar]
 Castelli, F. 1996, in M.A.S.S., Model Atmospheres and Spectrum Synthesis, eds. S. J. Adelman, F. Kupka, & W. W. Weiss, ASP Conf. Ser., 108, 85 [Google Scholar]
 Castelli, F. 2005a, Mem. Soc. Astron. Italiana Suppl., 8, 25 [NASA ADS] [Google Scholar]
 Castelli, F. 2005b, Mem. Soc. Astron. Italiana Suppl., 8, 34 [NASA ADS] [Google Scholar]
 Castelli, F., & Kurucz, R. L. 1994, A&A, 281, 817 [NASA ADS] [Google Scholar]
 Castelli, F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841 [NASA ADS] [Google Scholar]
 Catanzaro, G. 2014, eds. E. Niemczura, B. Smalley, & W. Pych in Determination of Atmospheric Parameters of B, A, F and GType Stars: Lectures from the School of Spectroscopic Data Analyses, (Cham: Springer International Publishing), 97 [Google Scholar]
 Cernetic, M., Shapiro, A. I., Witzke, V., et al. 2019, A&A, 627, A157 [CrossRef] [EDP Sciences] [Google Scholar]
 Chase, M. W. J. E. 1998, J. Phys. Chem. Ref. Data, Monograph, 9 [Google Scholar]
 Claret, A. 2000, A&A, 363, 1081 [NASA ADS] [Google Scholar]
 Collins, G. W. 1989, The fundamentals of stellar astrophysics II (New York: W. H. Freeman) [Google Scholar]
 Cowley, C. R., & Castelli, F. 2002, A&A, 387, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doerr, H. P., Vitas, N., & Fabbian, D. 2016, A&A, 590, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fontenla, J. M., Harder, J., Livingston, W., Snow, M., & Woods, T. 2011, J. Geophys. Res. Atmos., 116, D20108 [Google Scholar]
 Fontenla, J. M., Stancil, P. C., & Landi, E. 2015, ApJ, 809, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [Google Scholar]
 Grupp, F. 2004, A&A, 420, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubeny, I., & Lanz, T. 2017, ArXiv eprints [arXiv:1706.01859] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2015, Theory of stellar atmospheres: an introduction to astrophysical nonequilibrium quantitative spetroscopic analysis (Princeton Univ. Press), [Google Scholar]
 Huber, K. P., & Herzberg, G. 1979, Molecular Spectra and Molecular Structure: IV. Constants of Diatomic Molecules (Boston, MA: Springer) [Google Scholar]
 Husser, T.O., Wendevon Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koesterke, L., Allende Prieto, C., & Lambert, D. L. 2008, ApJ, 680, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 1969, ApJ, 156, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 1970, SAO Spec. Rep., 309 [Google Scholar]
 Kurucz, R. 1993, ATLAS9 Stellar Atmosphere Programs and 2 km/s grid. Kurucz CDROM, 13 [Google Scholar]
 Kurucz, R. L. 2005a, Mem. Soc. Astron. Italiana Suppl., 8, 14 [NASA ADS] [Google Scholar]
 Kurucz, R. L. 2005b, Mem. Soc. Astron. Italiana Suppl., 8, 86 [NASA ADS] [Google Scholar]
 Kurucz, R. L., Peytremann, E., & Avrett, E. H. 1974, Blanketed model atmospheres for earlytype stars, (Smithsonian Inst. Press) [Google Scholar]
 Lester, J. B., & Neilson, H. R. 2008, A&A, 491, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Magic, Z., Collet, R., Asplund, M., et al. 2013, A&A, 557, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marfil, E., Tabernero, H. M., Montes, D., et al. 2020, MNRAS, 492, 5470 [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Mészáros, S., Allende Prieto, C., Edvardsson, B., et al. 2012, AJ, 144, 120 [Google Scholar]
 Mihalas, D. 1978, Stellar atmospheres (San Fransisco: Freeman) [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (Courier Corporation) [Google Scholar]
 Neckel, H. 1999, Sol. Phys., 184, 421 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Solar Phys., 6, 2 [Google Scholar]
 Norris, C. M., Beeck, B., Unruh, Y. C., et al. 2017, A&A, 605, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F. A., Cristiani, S., Rebolo Lopez, R., et al. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, Proc. SPIE, 7735, 77350F [NASA ADS] [Google Scholar]
 Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407 [NASA ADS] [CrossRef] [Google Scholar]
 Ralston, A., & Rabinowitz, P. 1978, A First Course in Numerical Analysis, International series in pure and applied mathematics (McGrawHill: New York, NY) [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 9143, 914320 [Google Scholar]
 Riethmüller, T. L., Solanki, S. K., Berdyugina, S. V., et al. 2014, A&A, 568, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rutten, R. J. 2019, Sol. Phys., 294, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr., 90, 054005 [Google Scholar]
 Schwenke, D. W. 1998, Faraday Discus., 109, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, A. I., Schmutz, W., Schoell, M., Haberreiter, M., & Rozanov, E. 2010, A&A, 517, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Short, C. I., & Hauschildt, P. H. 2009, ApJ, 691, 1634 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F. 2012, Liv. Rev. Solar Phys., 9, 4 [NASA ADS] [Google Scholar]
 Tagirov, R. V., Shapiro, A. I., Krivova, N. A., et al. 2019, A&A, 631, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tatum, J. B. 1966, Publ. Dominion Astrophys. Observatory Victoria, 13, 1 [NASA ADS] [Google Scholar]
 Uitenbroek, H., & Criscuoli, S. 2011, ApJ, 736, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Woods, T. N., Chamberlin, P. C., Harder, J. W., et al. 2009, Geophys. Res. Lett., 36, L01101 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Solving the statistical equilibrium
MPSATLAS assumes LTE, which is a good enough approximation for photospheres in cool mainsequence stars. In LTE the populations of atomic and molecular levels are in thermal equilibrium and, thus, can be evaluated independently of the radiation field with the help of the SahaBoltzmann (SB) and GuldbergWaage equations (with the latter being the analogue of the SahaBoltzmann equation for molecular chemical equilibrium calculations, see, e.g., Tatum 1966).
Following ATLAS9, in the MPSATLAS code a set of equilibrium equations for all species together with conservation constraints is formulated. Then, using the SB equation, the set of equations are expressed in terms of neutral atom and electron number densities. For calculations without molecules, the species taken into account in the set of equilibrium equations are hardcoded and include the most important electron donors: H, He, C, Na, Mg, Al, Si, K, Ca and Fe. For calculations with molecules, the number of species, nmol, is set by the number of elements and molecules listed in the file ‘molecules.dat’, whose default version includes the atoms listed above and significantly more. In this file, not only the species are listed, but also the molecular equilibrium constants are provided together with their dissociation energies, D_{0}. The equilibrium constants are pretabulated as a function of temperature using coefficients from fits to the NISTJANAF tables (Chase 1998), and the dissociation energies were taken from Huber & Herzberg (1979). These pretabulated coefficients are used to evaluate the GuldbergWaage equation.
For both cases, the equilibrium equations can be rearranged, such that for each considered species, i, one has an equation f_{i}(n_{1}, ..., n_{i}, n_{total}) = 0, where n_{i} is the number density of the species i. To find the number densities, the matrix M_{ij} = ∂f_{i}/∂n_{j} is set up and the matrix equation
where Δn is the change in n, is solved using the method of triangular decomposition (Ralston & Rabinowitz 1978, Chap. 9). Using a NewtonRaphson technique, the solutions are iterated until the relative change for each species in number density, Δn_{i}/n_{i}, becomes less than 10^{−4}. A more detailed description of the underlying procedure can be found in Kurucz (1970).
Appendix B: Improvements for model calculations
Since calculating atmosphere models is an iterative process, it is useful to know when to stop. While in ATLAS9 there were no criteria that identified when the atmosphere model is close enough to RE, we implemented the following procedure. For the model convergence criteria, the greatest relative temperature adjustment in all atmosphere points has to be smaller than 10^{−5}. If this criterion is reached before the prescribed maximum number of iterations, it is considered that the atmosphere model has converged. This threshold can be easily changed. For reference, the greatest relative temperature adjustment is written out in the last iteration. Moreover, atmosphere models are calculated on a prescribed Rosseland mean opticaldepth, τ_{Ross}, grid. We improved the setup to make it more user friendly (for an example setting see Appendix G). The treatment of convective flux and overshoot in the code is mainly kept as originally implemented (for more details see Castelli et al. (1996)). The socalled approximate overshoot was tested in Castelli et al. (1997).
Appendix C: Starting model for RE calculations
For the model calculations, the initial 1D atmosphere model is rescaled using the ratio of the desired effective temperature value to its initial value on a Rosseland mean opticaldepth, τ_{Ross} grid. This is achieved using the temperature as a function of τ_{Ross} and applying
where the superscript ‘initial’ indicates the initial model temperature structure, and its effective temperature. Kurucz (1970, Chap. 2.12) showed that the Rosseland optical depth is suitable for such a rescaling. The rescaled temperature structure is converted back onto the column mass grid. This results in a first starting point model for which subsequently the RTE has to be solved.
Appendix D: Solving radiative transfer
The MPSATLAS code has two different radiative solvers implemented. In the original ATLAS9 implementation of the RTE solver by R. L. Kurucz, the scattering part of the source function is found iteratively (hereafter, iterative solver). To minimise the computational time, pretabulated matrices for the evaluation of the mean intensity on a fixed optical depth grid are used (for more details see Kurucz 1969, 1970). The alternative way to account for the scattering part of the source function is to use a Feautrier method as described in Mihalas & Mihalas (1984) (hereafter, Feautrier solver). In this approach, the secondorder transfer equation, derived using Feautrier variables together with upper and lower boundary conditions, is solved. The Feautrier method allows to set the number of viewing angles, μ, which are included in the calculation, to three, four or eight (Lester & Neilson 2008).
The two solvers have a different treatment of the source function at frequencies where strong lines are formed in the upper atmosphere, that is the case when the top atmospheric point has τ_{ν} > 0.2. Since for such frequencies the atmosphere is not sufficiently high in order to obtain an accurate solution of the RT and nonLTE effects are important, the Feautrier solver sets the source function and the flux to zero (see orange line in Fig. D.1). This differs from the iterative method that still solves the RT equation in such cases, but makes an approximation by considering a prolongation of the atmosphere with a source function, , equal to the value. An example calculation of the highresolution flux for a small wavelength interval obtained using the two different solvers is shown in Fig. D.1. The two solutions are both approximations, but lead to different emergent intensities, as can be seen by looking at the Ca II K line (λ = 393.478 nm). We note that the different treatment has a negligible effect when ODFs are used for wavelengths longward of 200 nm (see Fig. D.2).
Fig. D.1. Highresolution discintegrated flux, F_{ν}, calculated using the iterative solver and the Feautrier solver with 3 view angles. 
Fig. D.2. Ratio of the emergent intensities for two different limb positions with the ODF approach calculated using the Feautrier solver to the one using the iterative solver. 
We tested the performance of these two solvers. For that we compare the resulting effective temperature calculated from the emergent flux with the value set in the atmosphere model calculation, and we measure the computation time needed to reach a converged model. We found that, generally, the iterative solver is the most accurate but slower compared to the Feautrier solver. The optimal tradeoff between computational time and accuracy is achieved when using the Feautrier solver with three μ angles. Ideally, to obtain consistent results, the same solver should be used for the specific emergent intensities as is used for the radiative equilibrium calculations.
Appendix E: ODF and opacity calculations
The calculation and pretabulation of ODFs on a TP grid is done in three steps. First, the quantities needed for selecting which atomic and molecular lines to include and for calculating the line opacity are computed (e.g., molecular and atomic number densities and continuum opacity). Second, lines are selected and the highresolution opacity is calculated. Finally, the highresolution opacity is split into bins, sorted and averaged over the subbins.
In the first stage, the TP grid together with the element composition is set up. Subsequently, the equilibrium number densities are found in each TP point as described in Sect. A. Having this information the following quantities are obtained: (i) the number density over partition function, N_{j}/U_{j}, where j indicates the ionisation stage, are obtained up to j = 5 for all atoms, (ii) the total continuous opacity, κ_{c}, on a grid of 858 frequency points in the range from 1 nm–500 000 nm, (iii) the quantity , which is the ratio of the thermal velocity of a given element (with atomic mass m_{el}) to the speed of light, for microturbulence v_{turb} = 0 km s^{−1}, and (iv) the quantity
which gives an estimated measure of the line opacity relative to the minimum continuous opacity min(κ_{c}) in a given frequency bin multiplied by a cutoff, where the cutoff typically has a value of order 10^{−3}. In this ratio ρ is the density and the Doppler shift of the transition at the frequency v_{0} is
These quantities are required in the second stage for two purposes: the first being for the preselection of lines in order to reduce the computational cost, where only lines are selected that pass several conditions. The first condition is that the quantity nf_{max} has to be greater than unity, which ensures that without the knowledge of the oscillator strength there are enough particles in a given ionisation stage. For the final condition, the oscillator strength, f_{ij}, of the transition, the statistical weight of the ith level, g_{i}, and the Boltzmann factor are taken into account. This results in
This conditions ensures that the line core opacity is greater than a thousandth of the continuous opacity.
The second purpose being that the precalculated quantities are needed for calculating the linestrength and broadening, as they depend on the ionisation fraction and the Doppler shift. Namely, the line absorption coefficient as given in Eq. (4) is computed. Here, for all lines, except hydrogen lines, the VoigtProfile is used. For hydrogen lines the Stark broadened profile is used (Cowley & Castelli 2002). With the preselected lines and the lineprofiles, the detailed highresolution opacity is calculated on a wavelength grid from 8.9766 nm to 10 000 nm with a resolving power of R = 500 000. This results in the wavelength points
where n is the index of the grid points.
During the third stage, the highresolution opacity is split into wavelength bins, which are either on the standard Kurucz’s grid, or userdefined. The standard wavelength grid can be set by the keyword ‘binning off’ and is hardcoded. For a userdefined bin grid, an additional file ‘bingridsizes.dat’, that contains the bin borders, has to be provided. After selecting the bin range, the opacity in each bin is sorted and further split into subbins as described in Sect. 2.3. The number of subbins and their sizes are either as in the standard Kurucz’s configuration, which is automatically set if the standard bin grid is used, or a userdefined configuration, which has to be specified in the file ‘subbininfo.dat’. Finally, the opacity is averaged over the subbins, and written out.
Having generated a ODF table, the subbin opacity values can be read and further processed in module II and module III. In both modules the subbin opacity is added to the continuum opacity for each subbin separately. Subsequently, the RTE is solved. The resulting quantities, such as the flux H_{ν}, and its derivative in module II, or the surface flux, F_{ν}, or surface intensity, I_{ν}, in module III, need to be calculated for each frequency bin. Thus, an average over the weighted contributions from the subbins is used. The weights correspond to the wavelength interval of the subbins. As an example, the surface intensity in a bin i, I_{bin, i}, with Δλ_{i} that is centred at λ_{c} = (λ_{i + 1} + λ_{i}) ⋅ 0.5, is obtained by summing the calculated intensity in each subbin, weighted by the subbin width, w_{s}, as follows
Here, I_{λc} is the intensity at λ_{c} obtained using the opacity values along the atmosphere, κ_{tot, s, i}, which are the sum of the continuous opacity at λ_{c} and the averaged subbin opacity of the subbin s in each atmosphere point.
Appendix F: Molecular photodissociation
ATLAS9 includes photodissociation for CH and OH. This is achieved by using pretabulated crosssections and the corresponding partition functions. In addition, we implemented the NH dissociation opacity. Our implementation is, however, slightly different from that of CH and OH. The number densities for NH that are obtained from the equilibrium calculations are multiplied with the crosssections directly. Since a negligible fraction of NH might be in an excited state, this treatment potentially overestimates the opacity.
To test the effect of NH photodissociation opacity on the emergent intensities, we calculated the solar atmosphere model in RE with and without the NH photodissociation and subsequently the corresponding emergent intensities. The ratio of the emergent intensities is shown in Fig. F.1. It indicates that the difference is below 0.5% at disccentre and below 0.3% at the limb (μ = 0.05). Thus, we conclude that including NH is not sufficient to account for the missing opacity in the UV.
Fig. F.1. Ratio between emergent intensities for the Sun obtained including NH photodissociation and without it for two different viewing angle μ. 
Appendix G: Example input
To run MPSATLAS several input files are required. In an overall control file, ‘mpsa.control’, (shown in Fig. G.1) the user has to specify which module or modules should be executed and give the names of the input files. To set a module for calculations, a line with the keyword ‘calc_’ followed by either ‘odf’, ‘model’ or ‘flux’ has to be included. By adding an ‘off’ after the keyword, the module can be switched off again without deleting the line. In the example in Fig. G.1, the code will only calculate the emergent flux module.
Fig. G.1. Example of an ‘mpsa.control’ file. 
For every module there are three types of input files that need to be specified: (i) the input file that contains computational flags, (ii) a model atmosphere file, and (iii) the file that contains the opacity tables, where for the ODF table calculation the name of the output file is specified. The order of the lines in the control file is irrelevant, but to indicate that there are no more lines to read the keyword ‘endfile’ has to be set as the last line.
G.1. Input files
To start a particular module, it has to be activated in the ‘mpsa.control’ file. Subsequently, each module requires different settings in the ‘.input’ file, where some settings are common. The order of the keyword lines in the input file is not important, except for the last two lines which have to be ‘begin’ and ‘end’, to indicate the end of the input file and start the calculation. Example input files for each module are given in Fig. G.2–G.4. All three modules have most of the flags in common, while a few are specific for certain modules.
Fig. G.2. Example input for ODF generation (odf.input). 
Fig. G.3. Example input file for model calculations (model.input). 
Fig. G.4. Example input for specific emergent intensities (flux.input). 
G.1.1. Common settings
In any of the modules, molecules can be included or excluded by the keywords ‘molecules on’ or ‘molecules off’. If the molecules are ‘on’, a routine reads the file ‘molecules.dat’ which should be located in the INPUT folder.
For all modules the chemical composition needs to be specified as shown in Fig. G.2–G.4. The line starting with the keyword ‘abundance scale’ specifies the metallicity. The number following this keyword scales all elements heavier than helium. Consequently, a single scaling factor accounts for changes in metallicity, so that it is not necessary to change all elements individually when changing the metallicity. All element abundances can be set using the keyword ‘abundance change’ followed by the elemental number and the number fraction. While for hydrogen and helium the number faction is given as N_{element}/N_{total}, for all other elements it is given in log10(N_{element}/N_{total}). For consistency the settings in all input files should be kept the same throughout the calculation of a given star.
For the opacity calculations, the continuous opacity sources as listed in Sect. 2.1 are all included by default. If any of them need to be excluded they can be switched off by adding ‘opacity off’ followed by keyword of the corresponding opacity source. The continuous opacity sources are grouped as follows: (1) HI boundfree transitions (bf) and freefree transitions (ff), keyword ‘H1’ (2) bf and ff, keyword ‘H2+’ (3) H^{−} bf and ff, keyword ‘H’ (4) Rayleigh scattering on HI, keyword ‘Hray’ (5) HeI bf and ff, keyword ‘He1’ (6) HeII bf and ff, keyword ‘He2’ (7) He^{−} ff, keyword ‘He’ (8) Rayleigh scattering on He, keyword ‘Heray’ (9)–(11) various bf and ff transitions of C, N, O, Ne, Mg, Al, Si, Ca, Fe, the molecules CH, OH and NH, and their ions, keywords ‘Cool’, ‘Luke’, and ‘Hot’ (12) electron scattering, keyword ‘Elec’ (13) Rayleigh scattering on H_{2}, keyword ‘H2ray’. In addition, if the line opacity should be excluded in either the model calculations or the flux calculations, this can be achieved by including the line ‘opacity off Lines’.
The keywords ‘user defined binning’ switches between the original standard bin and subbin configuration used in Castelli (2005b) if set ‘off’, and a userdefined configuration if set ‘on’. Having set a user defined binning and subbinning in the input files, individual grids have to be provided in separate input files, where the number of grid points, their intervals, and the number of subbins, and their sizes are specified. Then, the code will allocate the requested number of bins and subbin. This switch sets the bin and subbin configuration for module I, but has to be consistent with the format of the ODF used in module II and III.
The two keywords ‘print’ and ‘punch’ control how much information is written out during a run. They are both followed by an integer. For ‘print’ the following values can be set:

0 for no print

1 for summary tables

2 prints the temperature corrections, and surfaces fluxes

3 prints everything listed as for point 2, and also τ_{ν}, S_{ν} and J_{ν} for each frequency

4 prints the quantities listed under point 2 and all opacities for each frequency
For ‘punch’ there are only four options:

0 nothing is written out

1 writes atmosphere model

2 writes punch 1 and in addition surface fluxes or intensities

5 writes punch 2 and also molecular number densities over partition functions
G.1.2. ODF generation specific settings
A specific switch for the ODF calculations can turn filter calculations on or off by using the keyword ‘filter’ (for more details on filter ODF see Anusha et al. 2021). If this switch is not set, the filter calculations are switched off by default.
Module I has additional keywords: in the default setting at least one ODF table using the microturbulence velocity of v_{turb} = 0.0 km s^{−1} is calculated. To obtain additional ODF tables, the keyword ‘velocities values’, followed by an integer indicating how many turbulent velocity values should be calculated, can be set. After the integer, the code expects the corresponding number of float numbers setting the microturbulence velocities. Moreover, by default the opacities in the subbin are averaged using the geometric mean. This can be changed to an arithmetic mean by setting the keyword combination ‘mean arithmetic’.
G.1.3. Model calculation specific settings
For the atmosphere model calculations, the desired stellar parameters and the Rosseland grid need to be specified. This is controlled by using the keyword ‘scale’ which has three different options:

no letter follows the ‘scale’ keyword; the code reads in the number of depth points of the model, the starting Rosseland mean opticaldepth, τ_{Ross}, in log10, the steplength in log10, the effective temperature, and the surface gravity, log g;

if the letter ‘x’ follows, instead of the steplength, the maximum τ_{Ross} is read

if the letter ‘b’ follows, instead of the number of depth points, the starting τ_{Ross}, the steplength, and the maximum τ_{Ross} is read, before T_{eff}, and log g
The keyword ‘iterations’ has to be set for the model calculations indicating the maximal number of iterations to be executed. We note that the calculations will be stopped if the model converges before this number is reached.
For this module the effective temperature and surface gravity of the starting model, which is used for the initial rescaling (see Sect. C), need to be specified. Thus, a line starting with the keywords ‘starting model’ followed by the ‘teff’ and ‘surface gravity’ values has to be included.
By default the convective flux and overshoot is turn off. To change that, the keyword ‘convection’ sets the convection flux calculations ‘on’ or ‘off’, and a line with ‘overshoot on’ turns on the additional convection flux in the overshoot region (Castelli et al. 1996). The mixinglength value is the number following ‘convection on’, and is usually set between 0.0 and 2.0. Moreover, turbulence can be set on/off, by including the line ‘turbulence’ followed by four numbers. The turbulent velocity, v_{turb}, is computed using the four numbers a, b, c and d after ‘turbulence on’, as follows:
where v_{snd} is the sound speed in centimetres per second, and the term d accounts for the convective turbulence, which is often the only contribution used on latetype stars.
G.1.4. Settings for model and flux calculations
For module II and module III in the case of two standard frequency grids, the keyword ‘frequencies’ determines if the ‘big’ or ‘little’ Kurucz grid is used. Moreover, if for the calculation only a particular wavelength range should be considered, the keyword ‘wave limits’ followed by the starting and end wavelength in nanometre can be specified. If this keyword is not set, the whole available wavelength range of the considered opacity table will be used.
G.1.5. Flux calculation specific settings
Module III has a switch to recalculate the electron number densities in a given atmosphere using the given elemental composition. This switch is in particular important if atmospheres from a 3D MHD cube where n_{e} are not precalculated are used. It can be switched on using ‘recalxne on’. Finally, setting the keyword ‘surface flux’ starts the calculations for the emergent flux at the surface, whereas, the keyword ‘surface intensity’, followed by the number of view angles and their values, results in emergent intensity calculations at the specified view angles.
G.2. Model structure files
For all modules a second input file is read. For the ODF table calculations (module I), this file contains the TP grid on which the ODF table should be calculated, while for modules II and III, it contains the starting atmosphere model and the atmosphere model for which the emergent intensities are calculated, respectively. The structure of the TP grid (odf.tpgrid) is simple, the first two lines are two integers indicating the number of T and P points. They are followed by the temperature values given in each line, and the pressure values.
The atmosphere model for modules II and III has the structure shown in Fig. G.5. The first line gives the number of models, the second line lists the number of depth points. The first line is always 1, except if models form 3D MHD cubes are calculated. Starting in the 3rd line, the model is given for each depth point using seven columns. The first column lists the column mass, followed by the temperature, pressure, electron number density, mean Rosseland opacity, radiation pressure, and turbulent velocity in each depth point. We only consider cases with constant turbulent velocity, while the code allows to use a depth dependent v_{turb}, when several ODF tables for a range of v_{turb} are read in.
Fig. G.5. Example of model atmosphere (model.start or flux.model). 
All Tables
All Figures
Fig. 1. Illustration of ODF generation in one example bin. (a) Detailed highresolution opacity in the bin from 420–422 nm. (b) Sorted opacity without the information of the wavelength, and the corresponding geometric mean values for 12 subbins, which were chosen as in Castelli (2005b). 

In the text 
Fig. 2. Original structure of codes that calculate ODF, model atmospheres, and emergent intensities. Each green bubble indicates a separate code that results in a separate executable. The listed tasks are the main procedures used in each programme. Grey background indicates procedures that are the same for all three calculations. 

In the text 
Fig. 3. Schematic structure of the code modules I–III. 

In the text 
Fig. 4. Input parameters for different code modules together with information on the output. 

In the text 
Fig. 5. Highresolution emergent flux for Vega around the Balmer jump region, in the range of 360 nm to 400 nm. (a) Flux calculated using interpolated opacities from pretabulated highresolution opacities on a TP grid and using opacities calculated for each depth point. (b) Ratio of the fluxes from (a). 

In the text 
Fig. 6. Same as Fig. 5, but for a smaller wavelength interval (395.5 nm to 398.5 nm), corresponding to the Hε line core. 

In the text 
Fig. 7. Ratios of the solar emergent intensity at the disccentre calculated using ODF tables that were generated with different cutoff factors for the line preselection. Here cutoff factor refers to the opacity of a given line relative to the continuum opacity. If this ratio drops below the cutoff factor, this particular line is not included in the computation of the intensity spectrum. The default value of the cutoff factor is 10^{−3}. 

In the text 
Fig. 8. Highresolution solar disccentre intensity in the range 5340 Å–5360 Å computed with MPSATLAS using different line lists together with data from the Hamburg atlas of the solar spectrum (Neckel 1999; Doerr et al. 2016): (a) computed intensity using Kurucz’s line list and the Hamburg atlas data, (b) intensity using the VALD3 line list and the Hamburg atlas data. 

In the text 
Fig. 9. Ratios of the emergent intensity at disccentre calculated using ODF tables with the VALD3 line list to the emergent intensity at disccentre calculated using ODF tables with Kurucz’s original line list. 

In the text 
Fig. 10. Top panel: speedup for ODF table generation as a function of the number of compute cores used (module I from Fig. 3). Bottom panel: speedup for atmosphere model calculation (module II from Fig. 3). Black lines in both panels indicates the speedup averaged over 5 runs. The grey shaded area indicates the spread in the 5 runs, and the red line indicates the ideal speedup. 

In the text 
Fig. 11. Flux comparison between MPSATLAS, KuruczATLAS9, and PHOENIX code for an Atype star with the effective temperature of 8000 K. Flux values are shown (top panel) together with the corresponding flux deviations (bottom panel) in % compared to the original Kurucz calculations. 

In the text 
Fig. 12. Same as in Fig. 11, but for a Ftype star with T_{eff} = 6500 K. 

In the text 
Fig. 13. Flux comparison between MPSATLAS, KuruczATLAS9, and PHOENIX code for a Ktype star with T_{eff} = 4000 K. Flux values are shown (the two top panels) together with the corresponding flux deviations (in the two bottom panels) in % compared to the original Kurucz calculations. 

In the text 
Fig. 14. Irradiance calculated for solar parameters compared to WHI. The light and dark grey shaded area indicates the minimum and maximum measurement error of the SIRS WHI (Solar Irradiance Reference Spectra for the 2008 Whole Heliosphere Interval), respectively (see, Woods et al. 2009, for a detailed description). Solar irradiance (top panels) calculated using MPSATLAS with three different assumptions: no convection, with convection but no overshoot, and with convection and overshoot. The flux deviations between WHI observed irradiance and models in % are shown in the bottom panels. 

In the text 
Fig. 15. Irradiance calculated for solar parameters compared to WHI. Same as in Fig. 14, except that the MPSATLAS model calculations were done without overshoot and using the Anders and Asplund elemental compositions. 

In the text 
Fig. 16. Irradiance calculated for solar parameters compared to WHI. Same as in Fig. 14, but without the Asplund composition, while in addition the solar flux provided on Kurucz’s website and the flux from PHOENIX are shown. 

In the text 
Fig. 17. Comparison of measured and calculated UV flux for Vega. The grey area indicates systematic and statistical measurement uncertainties. For comparison, MPSATLAS flux and PHOENIX flux are plotted. The modelled fluxes are scaled with the same factor and were obtained for the same effective temperature, metallicity, and surface gravity (Table 1). Top panel: absolute flux, with error. Bottom panel: deviations of modelled fluxes from measured flux in %. 

In the text 
Fig. 18. Comparison of measured and calculated flux for Vega. The observed flux consists of data from STIS. The grey area indicates systematic and statistical measurement uncertainties. For comparison, MPSATLAS flux and PHOENIX flux is plotted. The modelled fluxes are scaled with the same factor and were obtained for the same effective temperature, metallicity, and surface gravity. Top panel: absolute flux, with error. Bottom panel: deviations of modelled fluxes from measured flux in %. 

In the text 
Fig. 19. Comparison of STIS flux with MPSATLAS calculation with R = 400. (a) Hydrogen Balmer series. (b) H_{α} line. (c) Paschen series. (d)–(f) the deviation of the calculated flux to the measured flux in %. 

In the text 
Fig. D.1. Highresolution discintegrated flux, F_{ν}, calculated using the iterative solver and the Feautrier solver with 3 view angles. 

In the text 
Fig. D.2. Ratio of the emergent intensities for two different limb positions with the ODF approach calculated using the Feautrier solver to the one using the iterative solver. 

In the text 
Fig. F.1. Ratio between emergent intensities for the Sun obtained including NH photodissociation and without it for two different viewing angle μ. 

In the text 
Fig. G.1. Example of an ‘mpsa.control’ file. 

In the text 
Fig. G.2. Example input for ODF generation (odf.input). 

In the text 
Fig. G.3. Example input file for model calculations (model.input). 

In the text 
Fig. G.4. Example input for specific emergent intensities (flux.input). 

In the text 
Fig. G.5. Example of model atmosphere (model.start or flux.model). 

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.