Issue 
A&A
Volume 603, July 2017



Article Number  A27  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201628574  
Published online  03 July 2017 
NESSY: NLTE spectral synthesis code for solar and stellar atmospheres
^{1} Astrophysics Group, Imperial College London, London SW7 2AZ, UK
email: tagirovrinat@gmail.com
^{2} Institute for Astronomy, ETH Zürich, 8093 Zurich, Switzerland
^{3} Physical Meteorological Observatory Davos, 7260 Davos Dorf, Switzerland
^{4} Max Planck Institute for Solar System Research, 37077 Göttingen, Germany
Received: 22 March 2016
Accepted: 20 February 2017
Context. Physicsbased models of solar and stellar magneticallydriven variability are based on the calculation of synthetic spectra for various surface magnetic features as well as quiet regions, which are a function of their position on the solar or stellar disc. Such calculations are performed with radiative transfer codes tailored for modeling broad spectral intervals.
Aims. We aim to present the NLTE Spectral SYnthesis code (NESSY), which can be used for modeling of the entire (UVvisibleIR and radio) spectra of solar and stellar magnetic features and quiet regions.
Methods. NESSY is a further development of the COde for Solar Irradiance (COSI), in which we have implemented an accelerated Λiteration (ALI) scheme for comoving frame (CMF) line radiation transfer based on a new estimate of the local approximate Λoperator.
Results. We show that the new version of the code performs substantially faster than the previous one and yields a reliable calculation of the entire solar spectrum. This calculation is in a good agreement with the available observations.
Key words: methods: numerical / line: formation / opacity / radiative transfer / scattering / Sun: atmosphere
© ESO, 2017
1. Introduction
Stellar spectra are the primary instrument used to study the atmospheres of stars. Various chemical elements in the stellar atmospheres manifest themselves in the millions of spectral lines observed in the stellar spectra. Calculations of such spectra require the application of radiative transfer codes.
Among the simplifications that can be made in such codes the central one is that of local thermodynamic equilibrium (LTE). LTE implies that the energy level populations are not put out of thermal equilibrium by the radiation field. LTE modeling has significantly furthered the understanding of stellar spectra, but a nonLTE (NLTE) approach has to be used to attain a higher degree of physical realism (see e.g. the Historical Précis in Hubeny & Mihalas 2014).
However, NLTE approach renders the problem much more computationally demanding. In NLTE the populations of energy levels and the intensity of radiation affect each other and the equations governing the populations have to be solved simultaneously with the radiative transfer equation. Mathematically, in the general case, the problem is a search for the solution of a system of nonlinear algebraic equations coupled with an integropartialdifferential equation. To solve the equations a technique known as accelerated Λiterations (ALI) has proved to be very successful (Hubeny 2003). The central element of ALI is the notion of approximate Λoperator (ALO; see e.g. Cannon 1973; Scharmer 1981).
In this paper we present an update of the COde for Solar Irradiance (COSI). In particular, we describe the update of COSI ALIscheme. The acronym COSI has been used both for the NLTE radiative transfer code (see Haberreiter et al. 2008; and Shapiro et al. 2010, further referenced here as S2010) as well as to model solar irradiance variability (Shapiro et al. 2011; Thuillier et al. 2014b). Spectral synthesis and irradiance variability modeling are two adjacent, but separate problems. Therefore, we gave the code the new name NESSY (NLTE Spectrum SYnthesis) to make the purpose of the code unambiguous.
NESSY is a combination of two codes. The first is HMINUS (Hamann & Schmutz 1987; Schmutz et al. 1989), which solves the population equations simultaneously with the radiative transfer equation for a set of energy levels of chemical elements from hydrogen to zinc. The second code is FIOSS, which was adapted from SYNSPEC (Hubeny 1981) and first used in its merged form by Schmutz et al. (1991). FIOSS performs the final spectral synthesis taking into account a much wider set of energy levels than HMINUS and, consequently, millions of spectral lines. The roots of NESSY go back to a code which was originally designed to model the conditions in expanding atmospheres of massive hot stars (see e.g. Hamann & Schmutz 1987). The code has been subsequently applied for solar calculations.
The initial purpose of the code predisposed it to solve the radiative transfer equation in spherical symmetry and the comoving frame of reference (CMF, Mihalas et al. 1975). The part of the code’s ALI scheme for the radiative transfer in spectral lines was originally based on the core saturation ALO (Scharmer 1981) modified by Hamann (1985, 1986) for the case of an expanding atmosphere. However, for solar conditions, where the expansion velocities are much smaller than characteristic chaotic and convection velocities in the plasma, the original ALIscheme has proved to be inefficient. In this paper we present an update of ALO in NESSY and show that the ALIscheme based on the new ALO can handle both expanding and stationary case.
While the 1D radiative transfer calculations are gradually being replaced by a 3D approach (see e.g. Koesterke et al. 2008; Uitenbroek & Criscuoli 2011, for the intercomparison of 1D and 3D calculations), the 1D approach still remains at the forefront of the solar irradiance variability modeling. The main reason for this is that the 3D models suitable for the irradiance calculations are still not available (see, however, Criscuoli & Uitenbroek 2014, for the first attempt to base the irradiance modeling on the 3D radiative transfer). Consequently, the main application of NESSY has been calculating the spectra of the quiet Sun and magnetic features for the physicsbased modeling of solar irradiance variability. The code has been used in several reconstructions of the historic solar irradiance (Shapiro et al. 2011; Thuillier et al. 2012) as well as for modeling of the solar irradiance during the satellite era (Shapiro et al. 2012; Thuillier et al. 2014b). Recently, Shapiro et al. (2015) incorporated NESSY calculations into the SATIRE model (Spectral And Total Irradiance REconstruction, see Fligge et al. 2000; Krivova et al. 2003) of solar irradiance variability. A detailed discussion of the NESSYbased models of solar irradiance is given in the reviews by Ermolli et al. (2013) and Thuillier et al. (2014a).
The inefficiency of the old ALIscheme mentioned above made the calculation of strong spectral lines as well as far UV and radio spectral regions in the solar case rather cumbersome. Here we present an improved version of NESSY which reaches a converged solution of the NLTE problem significantly faster than the version published by Haberreiter et al. (2008) and S2010 and thus is better suited for irradiance modeling from UV to radio. This improvement not only advances the use of NESSY in 1D modeling, but also makes possible to employ it in the 1.5D regime (i.e. combining the intensities along parallel rays which sample the whole multidimensional model) for calculating the spectra emergent from 3D MHD models. The 1.5D approach has been demonstrated to give reliable results for spatially averaged spectra (Holzreuter & Solanki 2012), which justifies its application for modeling of the solar irradiance variability. In this paper we present the modification of the code, while its applications for the irradiance modeling will be discussed in forthcoming publications.
In Sect. 2 we give a brief description of the NLTE radiative transfer problem in general. In Sect. 3 we discuss the derivation of the new ALO in the framework of NESSY and compare its performance with the old ALO and nonaccelerated Λiterations. In Sect. 4 we compare the spectra synthesized with the old and new versions of the code to the available observations. Section 5 is devoted to the behavior and performance of the new ALO in stationary and expanding atmospheres. We also show the impact of different values of chaotic velocities. In Sect. 6 we analyze the departure coefficients for hydrogen levels calculated with NESSY. Finally, in Sect. 7 we summarize our results.
2. Problem description
2.1. General formulation
An NLTE radiative transfer problem in a medium with 1D geometry, and known temperature and density structures is described by the system of the radiative transfer and statistical balance equations.
For a given atmospheric density and temperature distribution on a discrete geometrical depth grid the formal solution of the radiative transfer equation for a given spectral line in the operator form is: (1)Here, S is the vector of line source functions at each depth. The source function depends on the elements of n – the vector of populations of all atomic levels. J refers to the vector of line intensities averaged over angle and spectral profile and Λ is the socalled Λoperator. Equation (1)can be written for the continuum intensity as well, in which case S, J and Λ would correspond to a certain continuum frequency.
The statistical balance system is a set of rate equations depending on the mean intensity J and augmented by a set of number density conservation equations. The entire system can be written as (2)Here, P(J) is the statistical balance matrix and an element of vector b is equal either to zero (for a rate equation) or to the abundance of the corresponding chemical element (for a conservation equation).
For rate equations in Eq. (2)the offdiagonal element of matrix P(J) describes the transition between levels i and j and it is equal to the sum of the corresponding radiative and collisional rates per second per particle: (3)where describes the collisions with species s.
The conservation equations in Eq. (2)express the mere fact that the abundance A_{k} of the chemical element k remains the same: (4)Here, n_{l} is the population of level l, l_{k} and u_{k} are the indices of the lowest and highest levels of chemical element k, is its number density and n^{tot} the total particle number density.
2.2. Accelerated Λiterations
The radiative transfer equation (Eq. (1)) and the statistical balance system (Eq. (2)) are coupled to each other via the vector of atomic level populations n and the mean intensity J. In order to consistently solve Eqs. (1)and (2)one can iterate between them. Such iterations are called Λiterations (Rutten 2003, p. 122; Hubeny & Mihalas 2014, p. 381). At step i they are described as follows (Koesterke et al. 1992): where S^{(i−1)} = S(n^{(i−1)}). The iterations start with an initial approximation for the source function and Λoperator corresponding to this approximation. At step i the populations n^{(i−1)} come from the previous Λiteration and are considered constant. The statistical balance system is solved for n^{(i)} and thus, in the case of regular Λiterations, it is linear with respect to populations.
In optically thick case Λiterations “exhibit a pathological behavior in that the solution stabilizes (i.e., relative changes of the source function become extremely small) long before the correct solution is reached” (Hubeny 2003). A technique known as ALI (Rutten 2003, p. 125; Hubeny & Mihalas 2014, p. 421) can be implemented in order to circumvent this behavior. According to ALI approach the Λoperator is split in two (Cannon 1973): (5)The operator is called ALO. It describes the radiative transfer in the optically thick regime, where the regular Λiterations have troubles converging. To obtain the iterative relation one replaces S^{(i−1)} in the first term of Eq. (5)with S^{(i)}, but leaves S^{(i−1)} in the second term. Regrouping terms in Eq. (5)and designating Λ_{i−1}S^{(i−1)} as J^{(i−1)} we arrive at the ALI update formula: (6)For the method to properly work may take any form as long as it is easily invertible and adequately describes the optically thick radiative transfer (Rutten 2003, p. 126).
2.3. Nonlinearity of the statistical balance system
In ALI instead of substituting Eq. (5) Eq. (6) one substitutes Eq. (6). This makes the statistical balance system nonlinear with respect to populations n^{(i)}. This nonlinearity has to be dealt with either by means of linearization of the statistical balance system (e.g. Rybicki & Hummer 1991, 1992) or by means of an iterative method for the solution of nonlinear algebraic systems of equations (e.g. Koesterke et al. 1992). In NESSY the nonlinearity is increased due to electron concentration entering the statistical balance system (and getting updated after each ALI) via the collision coefficient (Eq. (3)) and the charge conservation equation (7)where Z_{l} is the charge number corresponding to level l. Further nonlinearity comes from the socalled radiative bracket (Rutten 2003, p. 49; Hubeny & Mihalas 2014, p. 473) implemented in the code.
3. Accelerating Λiterations in NESSY
3.1. The choice of the approximate Λoperator
Up until now the part of NESSY ALI scheme treating spectral lines has been based on the core saturation operator (Scharmer 1981) optimized by Hamann (1985, 1986) for expanding atmospheres. Runs of the previous code version for the solar case showed that the ALI scheme operating on the core saturation operator basis exhibits a very slow convergence. Moreover, it contains the escape probability parameter γ (see Hamann 1985), which has to be adjusted during the course of Λiterations, which makes the functioning of the code less autonomous.
Therefore, we eliminated these drawbacks by replacing the core saturation operator. In doing so we tried to find an ALO:

that would work autonomously (i.e., no free parameter adjustment would be needed);

that would still work adequately in the case of an expanding atmosphere;

whose construction would be relatively simple in the framework of NESSY ALI scheme based on the Broyden method for the solution of nonlinear statistical balance equations (Koesterke et al. 1992).
In a stationary case the local operator is defined as the diagonal of the original Λoperator (see Eq. (1)): (8)The local operator has been successfully implemented in the literature dealing with radiative transfer problems in stationary media (Rybicki & Hummer 1991, 1992) and it has also been implemented in COSI for the radiative transfer in continuum.
According to the definition of Λoperator (see Eq. (1)) for a given depth point its diagonal element Λ^{∗} gives the fraction of radiation contributed by the local source function to the total radiation field at this depth point. Therefore 0 < Λ^{∗} < 1, and in an optically thin case Λ^{∗} ≈ 0 reducing Eq. (6) to (5) which describes the regular Λiteration (which does not show poor convergence in optically thin case). In optically thick case Λ^{∗} ≈ 1 and the offdiagonal elements are close to zero. Therefore, the requirements for the choice of approximate Λoperator mentioned in Sect. 2.2 are fulfilled.
3.2. Derivation of the local operator in NESSY
An element of the local operator describes the contribution of the source function at a given depth point to the radiation field at the same depth point. In the following subsection we discuss one of the possible ways to obtain an estimate for this contribution in CMF.
The NESSY code follows the raybyray formulation in Mihalas et al. (1975; further referenced here as MKH) to obtain the solution of the CMF radiative transfer equation (see Eq. (2.1) of MKH). The depiction of the raybyray geometry is given in Fig. A.1.
We start with the equation interrelating the vectors of meanintensitylike and fluxlike variables and at two consecutive frequencies k and k−1 along ray j. In matrix form (for all depths at once) together with the boundary conditions this equation reads (see Eq. (2.31) of MKH) (9)where is the source function. The elements of matrices , , are given in Appendix C of MKH. All matrices contain the information about the expansion velocity field. The presence of terms with differing frequency indices originates from the term containing the derivative of intensity with respect to frequency in the analytical CMF equation (see MKH). This derivative describes the Doppler shift of photons as they move from layer to layer. In the case of no velocity field the second and third term on the left hand side of Eq. (9)disappear recovering the raybyray form of the stationary spherical radiative transfer equation (Mihalas & Hummer 1974).
Fig. 1 Function C_{H}(i) (Eq. (11) for hydrogen levels only), corresponding to: the new ALI scheme for spectral lines based on the local operator (LOC), the old ALI scheme for spectral lines based on the core saturation operator (CSA), no acceleration in spectral lines (noALI). LSP denotes the acceleration with local operator in case when the spectral lines are treated in the statistical balance system with single precision instead of double. The small wiggles of LOC curve after the convergence limit has been reached is due to the update of the Eddington factors (see e.g. Hubeny & Mihalas 2014, p. 418) each 10th iteration. The same holds true for the smaller wiggles along CSA and noALI curves. The dashed lines denote the convergence limits used for the calculations with single and double precision. 
To calculate the elements of the local operator in CMF we follow the formalism of stationary spherical radiative transfer equation while keeping some information about the expansion. To this purpose we omit the second and third terms on the left side of Eq. (9). Thereby at the statistical balance stage of ALI we erase the information about the Doppler frequency shift and a part of the information about the expansion. Equation (9)is then reduced to (10)where the remaining matrix still depends on the velocity field (see Appendix C of MKH or Appendix D of this paper for details). The information about the Doppler shift of photons and the full information about the expansion still enter the ALI scheme at the stage of solution of the radiative transfer equation and the proposed approach proves to be working reasonably well in the case of an expanding atmosphere. Neglect of the information about photons coming from the neighboring frequencies at the stage of statistical balance solution gets iterated out during the ALI process, but the ultimate number of iterations required for the ALI scheme to converge depends on the expansion velocity field (see Sect. 5).
There has been a previous implementation of the local operator together with the CMF formalism in PHOENIX code. The corresponding derivation is given in Hauschildt & Baron (1999). It differs from our derivation significantly and therefore in Appendix A we provide our detailed procedure for obtaining the elements of the local operator. Moreover, the solution of the statistical balance equations in PHOENIX is treated with the linearization method (Rybicki & Hummer 1991, 1992), whereas in NESSY the Broyden method is used (Koesterke et al. 1992; Nicolas et al. 2011). Thereby in Appendix B we show how to take the local operator acceleration into account when calculating the Jacobian of the statistical balance system.
3.3. New ALIscheme vs. old ALIscheme vs. noALI
Fig. 2 Second hydrogen level departure coefficient in the FCHHTB solar atmosphere (Fontenla et al. 2009) for noALI and CSA cases after the indicated number of iterations. The converged LOC and LSP cases are shown for comparison. 
In this subsection we illustrate the convergence of the new and old ALI schemes used in the code for the radiative transfer in spectral lines. The calculations were performed using the exemplary FCHHTB solar atmosphere model (Fontenla et al. 2009) assuming no expansion. We employed the default NLTE atomic setup of the code consisting of the first 30 chemical elements of the periodic table. For each chemical element its neutral and first ionized state are considered. For each ionized state only the ground level is included. All in all this setup amounts to 114 energy levels (84 atomic levels + 30 ionized states), 217 boundbound transitions and 84 boundfree transitions.
We note here that the main purpose of NESSY is modeling the overall energy distribution (as opposed to modeling of detailed line profiles) in the solar and stellar spectra from UV to IR. It was shown in S2010 that the NLTE atomic setup described above is sufficient for this task in the solar case. In this work we confirm this conclusion (see Sect. 4).
The convergence at iteration i can be quantified by the function (11)Here, indices l and d run over all atomic levels and depths eligible for the calculation of C(i). Once C(i) reaches a preset convergence limit the NLTE problem described in Sect. 2.1 is considered to be solved.
In our calculations this limit was made equal to 10^{4} because once C(i) corresponding to the new ALI scheme (LOC), based on the local operator, reaches this limit the convergence rate significantly decreases. It can be seen in Fig. 1, where we show the function C(i) calculated, for the sake of simplicity, only for hydrogen levels included in our default atomic model (H^{−}, H ii and the first ten levels of H i). From here on we will denote this function as C_{H}(i). We attribute the aforementioned abrupt decrease to the limitations of the machine precision, and in order to support this hypothesis we also plotted the convergence curve for the LOC case, but having treated the line transitions in the statistical balance system with single precision instead of double (LSP). One can see that the behavior similar to the flattening of the LOC curve takes place at much higher values of C_{H}(i) close to 10^{2}. The fact that the LSP case is limited at the square root of the LOC case indicates that the limiting factor of the LOC convergence at 10^{4} is the numerical precision.
The convergence of almost all other chemical elements beyond hydrogen (not shown) reveals a similar behavior and approximately the same rate. However, because of the relatively large ionization potentials of helium and neon with respect to the other chemical elements included in our NLTE atomic model, He ii and Ne ii states show a nonconvergent behavior. The same holds true for He i levels because of the large excitation energy of the second helium level. The temperature in the FCHHTB model below the transition region is not high enough to yield the sufficient collision rates for the corresponding ionization and excitation states. Hence, the populations of He i levels as well as He ii and Ne ii states turn out to be much lower compared to the rest of the chemical elements included. This causes numerical problems for these levels and states and their subsequent nonconvergent behavior.
Figure 1 also shows C_{H}(i) for the old ALI scheme (CSA), based on the core saturation operator. This scheme requires 6459 Λiterations to reach the preset convergence limit. That is 224 iterations more than it takes the noALI scheme (no acceleration in spectral lines at all) to converge. As discussed below convergence does not guarantee correct solution. In fact noALI converges to a wrong solution. It should be noted here that for CSA case we used the same value of γ (see Sect. 3.1) as in S2010 which, however, could probably be more optimized for the case of no expansion.
To give the reader a better idea about the convergence rate and behavior of CSA and noALI schemes Fig. 2 shows the second hydrogen level NLTE departure coefficient for both cases after a certain number of iterations have been completed. The departure coefficient is calculated according to the Zwaan definition as , where n_{l} is the population of a given atomic level in NLTE and is the LTE population calculated with the NLTE electron density.
One can see that the CSA scheme converges to the LOC solution, while the noALI scheme stabilizes at a wrong solution. It is typical for the nonaccelerated Λiteration failure that the stabilized solution hardly differs from the starting value. Such behavior is seen in the left panel of Fig. 2 in the photosphere and chromosphere where the converged b_{2} hardly differs from its LTE value. It is interesting to learn from the right panel of Fig. 2 that even an inefficient acceleration (as it is in the CSA case with nonexistent or low expansion velocity field) is capable of slowly releasing Λiterations from this “stabilization trap” (Hubeny 2003).
Fig. 3 Comparison of COSI (S2010) and NESSY spectra to the ATLAS3 (Thuillier et al. 2004) and SIRS WHI (Woods et al. 2009) observations. COSI spectrum is available only from 90 nm onward. The gray areas show the uncertainties of the observed spectra. Longward of about 115 nm the uncertainties become indiscernible on the logarithmic scale. The dashed line in both panels designates the wavelength below which the COSI spectrum in S2010 was not shown because the authors have been aware of the nonconverged hydrogen level populations in the outer atmospheric layers. 
Fig. 4 Comparison of NESSY and SIRS WHI spectra. The dark and light gray areas designate the minimum and maximum uncertainty ranges of SIRS WHI (Woods et al. 2009). 
We note, that in the outermost atmospheric layers (>2.05 Mm) the CSA 1280 solution is the same as noALI solution after 1280 iterations and differs by about a factor of two from the LOC solution. This is important because the Lyman lines are formed in those layers. Consequently, the S2010 solar spectrum and Shapiro et al. (2011) solar irradiance reconstruction have been provided only for the wavelengths longward 125 nm, since in those studies the code was run for about 1000 CSAiterations (see next section).
Finally, Fig. 2 shows the LSP solution right after its 45th Λiteration, whereupon it stabilizes (see Fig. 1). As can be seen from Fig. 2, this solution is identical to the converged LOC solution. Comparing Figs. 1 and 2, we conclude that fewer iterations may be needed for the current NLTE atomic setup of the code to reach the correct doubleprecision solution.
4. Comparison with observations
In this section we present the comparison of two synthetic solar spectra to two observed solar reference spectra available today. The first synthetic spectrum (COSI) is calculated with the version of the code used in S2010. The second synthetic spectrum (NESSY) is calculated with the new version of the code. The reference spectra are ATLAS3 (Thuillier et al. 2004) and the Solar Irradiance Reference Spectrum (SIRS) for the 2008 Whole Heliosphere Interval (WHI) (Woods et al. 2009). We note that since S2010 publication the code has undergone a number of changes besides the upgrade of the ALIscheme (Sect. 3). In this section we present the NESSY spectrum that accounts for all of those changes. The changes are described in more detail in Appendix C. The synthetic spectra were computed with the FAL99C solar atmosphere model (Fontenla et al. 1999) and the set of abundances used in S2010.
The COSI spectrum in S2010 was calculated taking into account the line blanketing in UV and visible. To this purpose three iterations of ODFs (Haberreiter et al. 2008) were performed. Each ODF iteration requires a cycle of ALI. One more ALI cycle is needed to initiate the ODF iterations. In each ALI cycle the code was run up to 320 iterations. In this work, the NESSY spectrum was calculated using the same ODF iteration procedure. In each ALI cycle NESSY requires 74, 48, 24 and 3 Λiterations, respectively, to converge.
Figure 3 shows the comparison to the ATLAS3 and SIRS WHI spectra in 70−170 nm spectral range. In the top panel the high resolution (0.005 Å) spectra of COSI and NESSY spectra are convolved with the Gaussian profile corresponding to the spectral resolution of ATLAS3 (G. Thuillier, priv. comm.). In the bottom panel of Fig. 3 the 0.1 nm running means of COSI and NESSY spectra are shown. The running means are calculated to match the spectral resolution of SIRS WHI in the considered wavelength range.
Figure 3 shows, that NESSY better reproduces the continuum level and the core of the Lyα line. However, as can be seen (especially in the bottom panel of Fig. 3), most of the spectral lines in the considered wavelength range are either missing or appear in absorption instead of emission. This is because the level populations for these lines are calculated in LTE with respect to the ground level of a given chemical element (calculated in NLTE). To avoid an erroneously large emission the source function for these lines is made equal to the Planck function at the temperature minimum of the FAL99C model. In order to properly reproduce the missing or reversed lines the corresponding level populations must be treated in full NLTE.
Figure 4 shows the comparison of NESSY spectrum with SIRS WHI in 170−2400 nm range. In 170−310 nm wavelength range SIRS WHI was measured by SOLarSTellar Irradiance Comparison Experiment (SOLSTICE, McClintock et al. 2005) onboard the SOlar Radiation and Climate Experiment (SORCE, Rottman 2005) and has the resolution of 0.1 nm. In the 310−2400 nm range it was measured by the Spectral Irradiance Monitor (SIM, Harder et al. 2005) with a resolution decreasing from about 1 nm near 310 nm to about 34 nm near 1200 nm. Consequently, we convolve high resolution NESSY spectrum with the spectral window of SIM/SORCE (J. Harder, priv. comm.). Shortward 310 nm NESSY and SIRS WHI spectra were averaged over 1.5 nm bins to match the spectral resolution near the transition wavelength from SOLSTICE/SORCE to SIM/SORCE measurements. As can be seen from Fig. 4, NESSY spectrum is in a reasonable agreement with observations except for the two features around 320 nm and 345 nm, which are attributed to the inaccuracy of the molecular line lists used in this study (see S2010 for details).
We note that in S2010 the comparison of COSI and SIRS WHI spectrum in 400−700 nm range yielded an agreement within about 2%. However, the slightly smaller error in S2010 does not imply, that the previous version of the code was more accurate than the one presented in this paper. The output spectrum depends not only on the code, but also on the utilised semiempirical model of the atmosphere, abundances, atomic and molecular linelist, among others. In particular, the semiempirical models have been adjusted to yield a good agreement of the output spectrum employing some specific radiative transfer code (SRPM in the case of FAL99C model). In the forthcoming publications we will test the performance of NESSY in 1.5D radiative transfer based on 3D cubes produced by realistic MHD simulations.
5. Behavior and performance of the local operator in CMF
Fig. 5 Local operator element Λ^{∗} as a function of height for SA and FEA velocity profiles. The left panel shows 1−Λ^{∗} for Lyα and Lyβ lines. The right panel shows Λ^{∗} for Hα line. The figure is plotted for the default value of the Doppler broadening velocity km s^{1}. The black curve, corresponding to the right yaxis in both panels, shows the temperature profile of the FCHHTB solar atmosphere model (Fontenla et al. 2009). 
In this section we verify our derivation and implementation of local operator in CMF by means of investigating the behavior of its element as a function of height, expansion velocity and Doppler broadening velocity.
In Sect. 5.1 we show the dependence of local operator element on height in the exemplary FCHHTB solar atmosphere (Fontenla et al. 2009) and in Sect. 5.2 examine its response to a change of the expansion and Doppler broadening velocity as well as the corresponding response of the convergence rate.
We present the results of running NESSY for two cases codenamed as “stationary atmosphere” (SA) and “fast expanding atmosphere” (FEA). The FEA case implies the expansion law that is a monotonically increasing function outward (see the detailed description in Hamann & Schmutz 1987) with lower boundary value equal to 240 km s^{1} and the outermost expansion velocity equal to 2400 km s^{1} which represent the typical velocities in WolfRayet winds (Schmutz 1997). The SA case implies no expansion.
5.1. Local operator element as a function of height
The local operator element Λ^{∗} can be interpreted as the fraction of photons that originate at a certain depth point and remain at this point (see Sect. 3.1). Accordingly, its complement to one 1−Λ^{∗} yields the fraction of photons that penetrate the neighboring layers. Henceforth we will call the corresponding fractions of photons local and nonlocal, respectively.
Figure 5 shows 1−Λ^{∗} for Lyα and Lyβ lines and Λ^{∗} for Hα as functions of height for the SA and FEA velocity profiles. Because of high Lyα excitation energy most of the hydrogen atoms reside in their ground state everywhere in the atmosphere up to the transition region. It results in large opacities in Lyman series and, consequently, small nonlocal fraction of Lyα and Lyβ photons. The gradually decreasing density leads to increasing mean free path and therefore increasing nonlocal fraction of Lyα and Lyβ radiation. The local maximum of 1−Λ^{∗} at the temperature minimum in the SA case is caused by the finer height grid of the atmosphere model in this region. The closer to each other the grid points are the easier it is for the photons at each point to come from an adjacent layer. Around 2000 km, in the transition region, the temperature sharply rises and starts playing a dominant role in determining the behavior of the local operator element. The hydrogen ionization fraction there drastically increases making all hydrogen lines optically thin. Therefore, an increasingly larger portion of the photons in the corresponding layers becomes nonlocal and propagates through them unimpeded, eventually leaving the atmosphere. One can also see in the left panel of Fig. 5 the increase of 1−Λ^{∗} if instead of Lyα line Lyβ is considered. This behavior is consistent with the physical meaning of the local operator as Lyβ is less optically thick than Lyα due to the smaller spontaneous Lyβ transition coefficient.
The right panel of Fig. 5 shows Λ^{∗} for the Hα line. Because of high temperature and density in the lower photosphere Hα is optically thick there and Λ^{∗} = 1. Toward the temperature minimum Λ^{∗} decreases as the second hydrogen level becomes less populated at lower temperatures. Higher up the chromospheric temperature rise repopulates the second level thus making Hα more optically thick. The corresponding increase of Λ^{∗} together with the photospheric drop forms the socalled Hα opacity gap (Rutten & Uitenbroek 2012) in the vicinity of the temperature minimum. The decline of Λ^{∗} in the almost isothermal FCHHTB chromosphere is due to the gradual decrease of density.
5.2. Response to expansion and Doppler broadening velocity change
Fig. 6 Response of to different values of V_{dop} expressed in terms of the ratio to corresponding to km s^{1}. The figure is plotted for the SA velocity field. The linear scaling of with V_{dop} up to 2000 km is explained in detail in Appendix D. 
Fig. 7 Left panel: hydrogen convergence (see Sect. 3.3) in the SA case for the three indicated values of V_{dop}. Right panel: hydrogen convergence for the FEA and SA expansion velocity fields and the default Doppler broadening velocity. The dashed line designates the convergence limit. 
Figure 5 shows the decrease of Λ^{∗} in the FEA case. Owing to the gradient of the FEA velocity field the redshift of photons, seen from the consecutive atmospheric layers, transfers them away from the core of a line to the wing where the mean free path is longer. This results in the decrease of the local source function contribution into the radiation field at a given depth and therefore diminishes the local operator element. The disappearance of the 1−Λ^{∗} temperature minimum peak in the SA case mentioned in the previous section is therefore readily explained. The grid refinement no longer plays a significant role in facilitating the photon transfer since it is already facilitated by the atmosphere expansion.
In Fig. 6 we show the response of 1−Λ^{∗} for the Lyα line to a change of the Doppler broadening velocity. The Doppler broadening velocity is given by: (12)where is the average thermal velocity of gas particles (T− temperature, − average mass of gas particles) and V_{turb} is the average microturbulent velocity, that is, the average velocity of separate parcels of gas that have a size smaller than the photon mean free path. In the NLTE block of NESSY V_{dop} is treated as a free parameter and has the default nominal value of km s^{1}.
It is clearly seen that smaller Doppler broadening velocities cause the decrease of , while higher V_{dop} values produce the reverse effect. This is understood as relative thermal velocities of atoms along with turbulent motions of separate parcels of gas broaden line profiles and prolong the photons’ mean free path, hence delocalizing the radiative transfer.
Thereby, Figs. 5 and 6 clearly show how increasing Doppler broadening velocity, as well as increasing the absolute values and gradient of the expansion velocity field, works toward making the atmosphere more transparent.
Fig. 8 Left panel: NLTE departure coefficients (see Sect. 3.3) for the second and third hydrogen levels calculated with NESSY for the FCHHTB solar atmosphere model (Fontenla et al. 2009). Right panel: second hydrogen level population in LTE and NLTE cases expressed in terms of the ratio to the total number density. The black curve, corresponding to the left yaxis in both panels, shows the FCHHTB temperature profile. 
Fig. 9 Main processes (according to NESSY calculations) contributing to the statistical balance of hydrogen levels in photosphere, temperature minimum and chromosphere. Each arrow corresponds to a net rate. Wavy red arrows are transitions dominated by radiation, straight blue arrows are collisionally dominated transitions. Mixed arrows indicate a significant contribution both from collisions and radiation either everywhere in the region or in a part of it. 
In a more transparent atmosphere the interconnection between its layers is stronger and, as illustrated by Fig. 7, the local operator becomes less efficient. The figure shows the function C_{H}(i) (see Sect. 3.3) corresponding to some of the cases considered in this section. The slower convergence arising from the change of expansion or Doppler broadening velocity fields is readily seen.
6. Analysis of the hydrogen departure coefficients
In this section we discuss the effects arising from the nonlocality of radiation field (as calculated with NESSY) and leading to the departures from LTE. The left panel of Fig. 8 shows the second and third hydrogen level NLTE departure coefficients b_{2} and b_{3} (as defined in Sect. 3.3) for the FCHHTB atmosphere model (Fontenla et al. 2009). In this section we provide the interpretation of our result, while the comparison with the result of Rutten & Uitenbroek (2012) is given in Appendix E.
6.1. Photosphere (125−800 km) and temperature minimum (800−1050 km)
In the innermost layers (up to 125 km) the density and temperature are high enough to drive the populations of all levels to LTE by collisions. Upwards of 125 km the second hydrogen level becomes progressively underpopulated due to Balmer continuum penetrating from the hotter layers below and giving rise to 2 → p^{+} radiative net rate (see the left panel of Fig. 9). This net rate can be compensated only by Lyα collisions since the photosphere is too dense for Lyα photons to penetrate from below and counter the underpopulating Balmer continuum. Due to the decreasing temperature and density the Lyα collisions become less efficient, which leads to the drop of b_{2} in the upper photosphere.
Around the temperature minimum the density becomes low enough for photons to penetrate from the neighboring warmer layers and make the Lyα transition radiatively dominated causing b_{2}> 1. We note here that this effect may be small, but, since the second level population is also small around the temperature minimum, even a tiny amount of Lyα photons coming from the hotter surroundings can cause a significant second level departure from LTE.
The right panel of Fig. 8 shows the increase of absolute NLTE and LTE populations of the second level along the rise to the chromospheric temperature (900−1100 km). The rising temperature there enables the collisions to repopulate the second level and drive it toward LTE. Eventually, the LTE curve goes above the NLTE one which corresponds to b_{2}< 1 around 980 km. This slight underpopulation is caused by the loss of Lyα photons that contribute to the overpopulation at the temperature minimum.
6.2. Chromosphere (1050−2000 km), transition region (2000−2050 km) and the outermost layers (2050−2140 km)
The density in the chromosphere is not as high as that in the photosphere. Nevertheless since almost all hydrogen resides in its ground state it makes the chromosphere opaque for the lines in Lyman series.
The high opacity in Lyman lines together with the low chromospheric temperature gradient lead to the radiative balance in Lyman series whereas the photons corresponding to nonLyman transitions can “leak out”. This leakage creates radiative losses in nonLyman series and the underpopulation of n> 2 levels (shown only for n = 3 in the left panel of Fig. 8). The radiative losses result in the radiative net rate to the second level (see the right panel of Fig. 9, the downward wavy arrows) which, analogously to the photospheric case, can not be compensated by Lyα photons because of the radiative balance. Therefore, this net rate can only be balanced by Lyα collisions. If the collisions are not efficient enough it leads to b_{2}> 1.
We attribute the slight decrease of b_{2} from 1100 km to 1700 km to the slow increase of temperature changing by about 400 K in this region. The subsequent rise of b_{2} between 1800 km and 2025 km is due to the increased loss of Hα photons because of the surface effect. Around 2025 km the increasing temperature reenables collisions to drive the second level closer to LTE, thus forming the b_{2} peak in the transition region. In the outermost layers (2050−2140 km) all photons escape the atmosphere, but the loss of Balmer photons overweights the Lyα loss diverting the second level population away from its LTE value.
7. Summary
We have presented an advancement of the COde for Solar Irradiance (S2010), which we now call NLTE Spectrum SYnthesis (NESSY) code. Our work can be summarized as follows.
The ALIscheme previously employed in the code for the radiative transfer in spectral lines was not efficient. Recent solar applications of the code required to upgrade the ALIscheme. The new ALIscheme is based on the local operator (Olson et al. 1986). We described our derivation of it in CMF.
In the stationary case the old ALIscheme is able to converge to the correct solution and allows to avoid the pathological behavior of the nonaccelerated Λiterations. But in doing so it needs an unreasonably large number of iterations. We showed, that in the stationary case the new ALIscheme reaches the correct solution of the NLTE radiative transfer problem much faster and it is able to handle the case of an expanding atmosphere as well.
In the 70−170 nm we compared the old (COSI) and new (NESSY) synthesized spectra to the available reference spectra ATLAS3 (Thuillier et al. 2004) and SIRS WHI (Woods et al. 2009). NESSY can properly reproduce the core of Lyα line and the continuum level in this range. In the 170−2400 nm range we compared NESSY and SIRS WHI spectra. NESSY yields an agreement within the upper limit of the uncertainty of the measurements except for the two features around 320 nm and 345 nm where the discrepancy reaches 20%. We attribute it to the inaccuracies in the molecular line lists used.
The analysis of the behavior of the local operator element as a function of height has been performed along with the analysis of the local operator response to the change of the expansion and Doppler broadening velocities. The convergence rate in the cases of increased expansion or Doppler broadening velocities becomes lower, as expected.
The output spectra of the code have been already employed for modeling of solar brightness variations (Shapiro et al. 2015). In the future we plan to employ NESSY in 1.5D modeling of solar and stellar brightness variability. This has become possible only following the advancement of NESSY described in Sect. 3 was accomplished.
Acknowledgments
The research leading to this paper has received funding from SNF grant 200020153301. A.I.S. acknowledges funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007−2013) under REA grant agreement No. 624817. It also got a financial support from the Germany science ministry (BMBF, Bunderministerium fuer Bildung und Forschung) as part of the ROMIC program (Role Of Middleatmosphere In Climate). The authors thank Robert Rutten for inspiring discussions and providing his calculations in digital form and the referee S. A. Korotin for the helpful review.
References
 Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Criscuoli, S., & Uitenbroek, H. 2014, ApJ, 788, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Ermolli, I., Matthes, K., Dudok de Wit, T., et al. 2013, Atmos. Chem. Phys., 13, 3945 [Google Scholar]
 Fligge, M., Solanki, S. K., & Unruh, Y. C. 2000, A&A, 353, 380 [NASA ADS] [Google Scholar]
 Fontenla, J., White, O. R., Fox, P. A., Avrett, E. H., & Kurucz, R. L. 1999, ApJ, 518, 480 [NASA ADS] [CrossRef] [Google Scholar]
 Fontenla, J. M., Curdt, W., Haberreiter, M., Harder, J., & Tian, H. 2009, ApJ, 707, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Haberreiter, M., Schmutz, W., & Hubeny, I. 2008, A&A, 492, 833 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R. 1985, A&A, 148, 364 [NASA ADS] [Google Scholar]
 Hamann, W.R. 1986, A&A, 160, 347 [NASA ADS] [Google Scholar]
 Hamann, W.R., & Schmutz, W. 1987, A&A, 174, 173 [NASA ADS] [Google Scholar]
 Harder, J. W., Fontenla, J., Lawrence, G., Woods, T., & Rottman, G. 2005, Sol. Phys., 230, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 1999, J. Comput. Appl. Math., 109, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Holzreuter, R., & Solanki, S. K. 2012, A&A, 547, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubeny, I. 1981, A&A, 98, 96 [NASA ADS] [Google Scholar]
 Hubeny, I. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 17 [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton, NJ: Princeton University Press) [Google Scholar]
 Janicki, C. 1990, Comput. Phys. Comm., 60, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Jefferies, J. T. 1968, Spectral line formation (Walthan, Mass: Blaisdell) [Google Scholar]
 Johnson, L. C. 1972, ApJ, 174, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Koesterke, L., Hamann, W.R., & Kosmol, P. 1992, A&A, 255, 490 [NASA ADS] [Google Scholar]
 Koesterke, L., Allende Prieto, C., & Lambert, D. L. 2008, ApJ, 680, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Krivova, N. A., Solanki, S. K., Fligge, M., & Unruh, Y. C. 2003, A&A, 399, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McClintock, W. E., Snow, M., & Woods, T. N. 2005, Sol. Phys., 230, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar atmospheres, 2nd edn. (San Francisco: W.H. Freeman and Co.) [Google Scholar]
 Mihalas, D., & Hummer, D. G. 1974, ApJS, 28, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1975, ApJ, 202, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Nicolas, S., Bigarré, L., & Paletou, F. 2011, A&A, 527, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Rottman, G. 2005, Sol. Phys., 230, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres (The Nertherlands: Ultrecht University) [Google Scholar]
 Rutten, R. J., & Uitenbroek, H. 2012, A&A, 540, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1992, A&A, 262, 209 [NASA ADS] [Google Scholar]
 Scharmer, G. B. 1981, ApJ, 249, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Schmutz, W. 1997, A&A, 321, 268 [NASA ADS] [Google Scholar]
 Schmutz, W., Hamann, W.R., & Wessolowski, U. 1989, A&A, 210, 236 [NASA ADS] [Google Scholar]
 Schmutz, W., Leitherer, C., Hubeny, I., et al. 1991, ApJ, 372, 664 [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]
 Shapiro, A. I., Schmutz, W., Rozanov, E., et al. 2011, A&A, 529, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shapiro, A. V., Rozanov, E., Shapiro, A. I., et al. 2012, Atmos. Chem. Phys., 12, 3181 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, A. I., Solanki, S. K., Krivova, N. A., Tagirov, R. V., & Schmutz, W. K. 2015, A&A, 581, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thuillier, G., Floyd, L., Woods, T. N., et al. 2004, in Solar Variability and its Effects on Climate. Geophysical Monograph 141, eds. J. M. Pap, P. Fox, C. Frohlich, et al. (Washington DC American Geophysical Union), 171 [Google Scholar]
 Thuillier, G., Deland, M., Shapiro, A., et al. 2012, Sol. Phys., 277, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Thuillier, G., Melo, S. M. L., Lean, J., et al. 2014a, Sol. Phys., 289, 1115 [NASA ADS] [CrossRef] [Google Scholar]
 Thuillier, G., Schmidtke, G., Erhardt, C., et al. 2014b, Sol. Phys., 289, 4433 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [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: Derivation of the local operator for spectral lines
In this appendix we provide the details of the local operator element derivation in the framework of NESSY. Once the frequency shift terms are discarded in Eq. (9)it is formally reduced to the stationary radiative transfer equation in spherical symmetry. The solution of radiative transfer equation in spherical symmetry is done by means of the raybyray method. The geometry of this method is depicted in Fig. A.1.
Our derivation shows that the estimate of the Λoperator elements (for the radiative transfer in spectral lines) can be obtained in the same way as the value of line profile averaged mean intensity can be obtained from the meanintensitylike function (see Eq. (2.19) of MKH).
Appendix A.1: Angular averaging
We started by inverting Eq. (10)and using the designation . The solution reads (A.1)Here, d refers to a radial depth point, NI_{j} is the number of intersections of ray j with the spherical shells (see Fig. A.1), and (A.2)where , are the opacity and emissivity at the center of a line, φ_{k} are the normalized line profile values (we assume a Gaussian profile with complete frequency redistribution), and , are the background continuum opacity and emissivity at the corresponding line frequencies. We omitted the ray index from the designation of the source function as the latter is dependent only on radial depth in our calculations. The indices j and d have to be consistent in Eq. (A.1)because is defined at the intersections of the rays with the spherical shells (see Fig. A.1).
Fig. A.1 Geometry of the raybyray method for the solution of the radiative transfer equation in spherical symmetry. ND is the total number of radial depth points (or spherical shells) and NR_{ND} is the total number of rays. NR_{ND} is always greater than ND by the number of rays intersecting the vertical axis below the radius of the innermost spherical shell. 
We then use expression (7207) in (Mihalas 1978, p. 254) to express mean intensity J_{dk} at frequency k in terms of : (A.3)Here, NR_{d} designates the number of rays intersecting the spherical shell with index d and are the appropriate quadrature weights corresponding to depth d.
Expanding the summation in (A.1) up to ND (the total number of radial depth points) by putting for n>NI_{j} we substitute (A.1) into (A.3) arriving at (A.4)Using the definition of Λoperator (Eq. 1), we can express J_{dk} in terms of S_{nk}: (A.5)Changing the order of summation and taking the source function out of the sum over j in Eq. (A.4)we compare it to Eq. (A.5)and conclude that (A.6)
Appendix A.2: Line profile averaging
Expression (A.6)gives the matrix values of J_{dk}producing Λoperator which depends on frequency. However, our ultimate goal is to obtain an expression for the diagonal elements of the frequencyindependent Λoperator corresponding to mean intensity averaged over the line profile for exactly this mean intensity enters the statistical balance system in the case of boundbound transitions.
To proceed with the derivation we have to neglect the values of and in expression (A.2)which make the source function frequency dependent. The omission of continuum opacities and emissivities is an additional simplification in our ALI scheme that gets iterated out along with the simplification of Eq. (10). Expression (A.2)therefore becomes reduced to the line source function S_{n} = η_{n}/χ_{n} (index l is now omitted).
To find the elements of J_{dk}producing Λoperator we average Eq. (A.5)over the line profile: (A.7)where NF is the number of frequency points in the line profile and : (A.8)Using (A.6)we arrive at the final result: (A.9)In the line radiation transfer part of our ALI scheme we take the values of Λ^{∗} to be the diagonal elements of matrix (A.9): (A.10)We note, that in practice we do not invert matrix in order to derive (A.10). Instead we use the algorithm described in Appendix B of Rybicki & Hummer (1991) to directly obtain the diagonal of from .
Appendix B: The derivative of line radiative bracket with the local operator
Because of the numerical precision limitations the rate equations in the statistical balance system have to be written in terms of the radiative bracket (Rutten 2003, p. 49; Hubeny & Mihalas 2014, p. 473). If any Newtonlike method is used to solve the statistical balance system (as is the case with NESSY, see Sect. 3.1) it requires the calculation of its Jacobian. This calculation is reduced to the evaluation of the radiative bracket derivative.
It is crucial to take into account the influence of acceleration on this derivative. Otherwise the initial approximation to the solution of the statistical balance system can be too far away from the correct solution. As our numerical experiments showed it leads to numerical instabilities in the ALI scheme preventing it from convergence. In this appendix we show how the local operator element enters the derivative of the radiative bracket for boundbound transitions.
The radiative bracket Z_{ul} is used to calculate the radiative net rate in a numerically advantageous way. In the case of scattering with complete frequency redistribution (adopted in the NLTE block of NESSY) we can write the radiative net rate as follows: (B.1)Here, n_{u} and n_{l} are populations of the upper and lower line levels, R_{lu} and R_{ul} are upward and downward radiative rates per second per particle, A_{ul} is the Einstein coefficient for spontaneous emission, S_{ul} is the monochromatic line source function and is the mean line intensity averaged over the spectral profile. The offdiagonal elements of the rate matrix in case of line transitions therefore take the form: where C_{lu} and C_{ul} are the upward and downward collisional rates.
Consider the ith equation of the statistical balance system (see Eq. (2)): (B.2)An element of the Jacobian of system (2)is by definition a partial derivative of the multivariable function F_{i} with respect to jth level population: (B.3)Hence the elements of Jacobian corresponding to line transitions are comprised of derivatives of the upper rate matrix elements (B.2) since the diagonal elements are expressed in terms of the nondiagonal ones and the lower elements (B.3) normally depend only on the electron concentration.
For the regular Λiterations the nonlinearity of the radiative bracket comes only from the denominator of . In this case († denotes the quantities taken from previous Λiteration) while S_{ul} is expressed in terms of current level populations, thereby: (B.4)When ALI is applied the nonlinearity comes from the numerator as well and the radiative bracket derivative can be conveniently expressed in terms of a correction factor for Eq. (B.4): (B.5)where is the element of the local operator at a given depth for the transition in question (Eq. (A.10)).
Expression (B.5)remains valid for any diagonal approximate Λoperator and under the assumption of partial frequency redistribution. In the latter case the spontaneous emission profile enters expression (B.1)as an additional multiplier which does not change the ultimate result.
When the optical depth is low, the solution is not affected by the ALI correction factor as . When the optical depth is high and and therefore the correction factor becomes significant. In the lower photosphere, however, the solution of the statistical balance system readily converges to LTE populations and although the correction factor there approaches zero it eventually does not have any influence on the convergence process.
Appendix C: NESSY spectrum vs. COSI spectrum
Fig. C.1 Relative deviation of the NESSY spectrum from the one calculated with COSI. COSI spectrum is available only from 90 nm onward. (NESSY − COSI)/COSI ratio is shown for the three consecutive changes in the code that had the largest impact on the resulting spectrum: the update of ALIscheme (ALI), implementation of the new hydrogen collision coefficients (ALI + COL) according to Johnson (1972) and the correction of an error in the treatment of the reaction (C.1)(ALI + COL + C.1). All spectra have been calculated using the FAL99C model of solar atmosphere (Fontenla et al. 1999) and smoothed with the ATLAS3 Gaussian (see Sect. 4) on a 0.5 Å resolution grid. The shaded part of the spectrum was not published by S2010. Note the logarithmic scale of the xaxis in the bottom panel. 
Since its last publication in S2010 the code has undergone a number of changes. We will refer to the version of the code published in S2010 by using its old name COSI (see Sect. 1). The changes can be summarized as follows:

implementation of the ALIscheme based on the local operator(Sect. 3);

implementation of hydrogen collision coefficients according to Johnson (1972); up until now the hydrogen collisions have been calculated using the formulas given in Jefferies (1968);

correction of an error in the treatment of one of the reactions governing the population of H^{−}: (C.1)

correction of an error in the treatment of the stimulated freefree emission in the spectral synthesis block of the code;

implementation of new Gaunt factors for hydrogen freefree emission (Janicki 1990).
The last change in the list was made for more accurate calculations in the radio wavelength range and it has a negligible effect on the other ranges. The effects of the other four changes are shown in Fig. C.1, where the relative deviation of NESSY spectrum from the one calculated with COSI is shown.
The deviation introduced by the implementation of the new ALIscheme is indicated in Fig. C.1, as are the added effects of the newly implemented hydrogen collision coefficients and the correction of an error in the treatment of the reaction (C.1). The plateau between 400 nm and 1000 nm reflects the change in the concentration of H^{−}. The 0.25% deviation of all three curves longward of approximately 1600 nm is caused by the correction of an error in the treatment of the stimulated freefree emission in the spectral synthesis block of NESSY (the fourth item on the list).
Appendix D: The linear scaling of 1 – with V_{dop}
The behavior of the curves up to 2000 km in Fig. 6 can be understood as follows. According to expressions (C3) and (C4) of MKH in optically thick case matrix in Eq. (9)is diagonal and its element reads: (D.1)where d stands for the depth index and is a dimensionless quantity given by expression (2.30) of MKH. It follows from Eqs. (2.28) and (2.3) of MKH that scales as: (D.2)with V_{exp}(d) being the expansion velocity field at depth point d and χ_{dk} – the opacity at depth point d and frequency k.
Fig. E.1 NLTE departure coefficients b_{2} and b_{3} (see Sect. 3.3) calculated with RH (Rutten & Uitenbroek 2012, Fig. 10) and NESSY for the FCHHTB solar atmosphere model (Fontenla et al. 2009). NOLBKG curve is the default calculation with the NLTE block of NESSY (see Fig. 8, left panel) without accounting for the effect of line blanketing. LBKG curve is the output of NESSY NLTE block taking into account the effect of blanketing by all lines within 100−1000 nm wavelength range. LOPA is the “low opacity” curve not taking into account the effect of line blanketing and corresponding to V_{dop} = 2000 km s^{1} when values are 400 times larger than those of for km s^{1}. 
Up to 2000 km Lyα opacity is high enough (see Fig. 5) for in (A.9)to be approximated with its Taylor expansion up to the first order terms. Doing so we obtain: (D.3)where we used the designation (A.10). In the NLTE block of NESSY the line opacities have the Doppler profile: (D.4)where Δν_{D} = ν_{0}V_{dop}/c is the Doppler width, ν_{0} is the central line frequency, c is the speed of light and x_{k} = (ν_{k}−ν_{0})/Δν_{D} is the relative frequency shift from the line center. Therefore, from (D.2)−(D.4)we conclude that (D.5)in an optically thick regime.
Appendix E: Comparison of b_{2} and b_{3} hydrogen departure coefficients: NESSY vs. RH
The coefficients b_{2} and b_{3} for the FCHHTB solar atmosphere model (Fontenla et al. 2009) were also calculated by Rutten & Uitenbroek (2012; see their Fig. 10) with RH code (Uitenbroek 2001). We show their result in Fig. E.1 along with b_{2} and b_{3} calculated with NESSY.
Our interpretation of the origin of the b_{2} peak at the temperature minimum is the same as that of Rutten & Uitenbroek (2012). However, our b_{2} value there is lower by a factor of 40. Such a difference can probably be attributed to a different treatment of the Lyα profile in the ALI schemes of NESSY and RH. In NESSY we assume complete frequency redistribution and Doppler profile for the Lyα line (the VanderWaals and Stark broadening is, nevertheless, taken into account in the final spectral synthesis calculated after ALI has converged). A broader profile of Lyα would cause photons to penetrate from deeper and warmer levels and contribute to larger values of b_{2} at the temperature minimum.
To verify this hypothesis we artificially broadened the Lyα profile by increasing the V_{dop} value (see Sect. 5.2) in the NLTE block of NESSY. It follows from Eq. (D.5), that for atmospheric regions where Lyα is optically thick the V_{dop} = 2000 km s^{1} corresponds to the values of that are 400 times larger than values for km s^{1}. We certainly understand, that V_{dop} = 2000 km s^{1} is unrealistic and we consider it here with the only purpose to illustrate the effect of value on the population of the second hydrogen level. We note, that value strongly depends on the assumed approximation for the line profile (see Eq. (A.9)) and especially on the treatment of the line wings. The resulting departure coefficients are plotted in Fig. E.1 (LOPA curves) and one can see that larger value leads to a much better agreement between our results and that of Rutten & Uitenbroek (2012) at the temperature minimum. The agreement is better in the middle and upper photosphere as well. The latter happens due to the Lyα photons being able to penetrate from lower photosphere and counter the underpopulation of the second level by Balmer continuum (see Sect. 6.1).
Another significant difference between our calculations and that of Rutten & Uitenbroek (2012) is the b_{2} value in the lower chromosphere (where the FCHHTB temperature profile has a plateau). While b_{2} equals to one according to Rutten & Uitenbroek (2012), our calculations yield a significant deviations from LTE due to the photon loss in nonLyman series (see Sect. 6.2). Such an effect may arise due to the different treatment of Balmer ionization. The more efficient Balmer ionization is, the more photons are lost in the subsequent cascade downward transitions (see Fig. 9, right panel). This effect works in both directions, since Balmer ionization also counteracts the second level overpopulation. However, the net effect of the interplay between the Balmer ionization and the downward cascade transitions eventually depends on the relative efficiency of the transitions involved.
To illustrate this point we show in Fig. E.1 the departure coefficients calculated taking into account line blanketing (LBKG curves) as discussed in Haberreiter et al. (2008) and S2010. One can see that line blanketing leads to larger b_{2} values in the lower chromosphere. As has been mentioned above, the illumination of the chromosphere by the Balmer continuum contributes to counteracting the second hydrogen level overpopulation. Consequently, the damping of Balmer continuum by blanketing decreases this contribution increasing b_{2} in the chromosphere.
Another important mechanism which has a potential to strongly affect the statistical equilibrium calculations of hydrogen, but is not considered in the present study, is the photon loss due to blends of metal lines and the Lyα line (Schmutz 1997). In particular, such a photon loss may lead to a compensation of the net rate to the second level of hydrogen in the lower chromosphere and push the b_{2} value to unity.
All Figures
Fig. 1 Function C_{H}(i) (Eq. (11) for hydrogen levels only), corresponding to: the new ALI scheme for spectral lines based on the local operator (LOC), the old ALI scheme for spectral lines based on the core saturation operator (CSA), no acceleration in spectral lines (noALI). LSP denotes the acceleration with local operator in case when the spectral lines are treated in the statistical balance system with single precision instead of double. The small wiggles of LOC curve after the convergence limit has been reached is due to the update of the Eddington factors (see e.g. Hubeny & Mihalas 2014, p. 418) each 10th iteration. The same holds true for the smaller wiggles along CSA and noALI curves. The dashed lines denote the convergence limits used for the calculations with single and double precision. 

In the text 
Fig. 2 Second hydrogen level departure coefficient in the FCHHTB solar atmosphere (Fontenla et al. 2009) for noALI and CSA cases after the indicated number of iterations. The converged LOC and LSP cases are shown for comparison. 

In the text 
Fig. 3 Comparison of COSI (S2010) and NESSY spectra to the ATLAS3 (Thuillier et al. 2004) and SIRS WHI (Woods et al. 2009) observations. COSI spectrum is available only from 90 nm onward. The gray areas show the uncertainties of the observed spectra. Longward of about 115 nm the uncertainties become indiscernible on the logarithmic scale. The dashed line in both panels designates the wavelength below which the COSI spectrum in S2010 was not shown because the authors have been aware of the nonconverged hydrogen level populations in the outer atmospheric layers. 

In the text 
Fig. 4 Comparison of NESSY and SIRS WHI spectra. The dark and light gray areas designate the minimum and maximum uncertainty ranges of SIRS WHI (Woods et al. 2009). 

In the text 
Fig. 5 Local operator element Λ^{∗} as a function of height for SA and FEA velocity profiles. The left panel shows 1−Λ^{∗} for Lyα and Lyβ lines. The right panel shows Λ^{∗} for Hα line. The figure is plotted for the default value of the Doppler broadening velocity km s^{1}. The black curve, corresponding to the right yaxis in both panels, shows the temperature profile of the FCHHTB solar atmosphere model (Fontenla et al. 2009). 

In the text 
Fig. 6 Response of to different values of V_{dop} expressed in terms of the ratio to corresponding to km s^{1}. The figure is plotted for the SA velocity field. The linear scaling of with V_{dop} up to 2000 km is explained in detail in Appendix D. 

In the text 
Fig. 7 Left panel: hydrogen convergence (see Sect. 3.3) in the SA case for the three indicated values of V_{dop}. Right panel: hydrogen convergence for the FEA and SA expansion velocity fields and the default Doppler broadening velocity. The dashed line designates the convergence limit. 

In the text 
Fig. 8 Left panel: NLTE departure coefficients (see Sect. 3.3) for the second and third hydrogen levels calculated with NESSY for the FCHHTB solar atmosphere model (Fontenla et al. 2009). Right panel: second hydrogen level population in LTE and NLTE cases expressed in terms of the ratio to the total number density. The black curve, corresponding to the left yaxis in both panels, shows the FCHHTB temperature profile. 

In the text 
Fig. 9 Main processes (according to NESSY calculations) contributing to the statistical balance of hydrogen levels in photosphere, temperature minimum and chromosphere. Each arrow corresponds to a net rate. Wavy red arrows are transitions dominated by radiation, straight blue arrows are collisionally dominated transitions. Mixed arrows indicate a significant contribution both from collisions and radiation either everywhere in the region or in a part of it. 

In the text 
Fig. A.1 Geometry of the raybyray method for the solution of the radiative transfer equation in spherical symmetry. ND is the total number of radial depth points (or spherical shells) and NR_{ND} is the total number of rays. NR_{ND} is always greater than ND by the number of rays intersecting the vertical axis below the radius of the innermost spherical shell. 

In the text 
Fig. C.1 Relative deviation of the NESSY spectrum from the one calculated with COSI. COSI spectrum is available only from 90 nm onward. (NESSY − COSI)/COSI ratio is shown for the three consecutive changes in the code that had the largest impact on the resulting spectrum: the update of ALIscheme (ALI), implementation of the new hydrogen collision coefficients (ALI + COL) according to Johnson (1972) and the correction of an error in the treatment of the reaction (C.1)(ALI + COL + C.1). All spectra have been calculated using the FAL99C model of solar atmosphere (Fontenla et al. 1999) and smoothed with the ATLAS3 Gaussian (see Sect. 4) on a 0.5 Å resolution grid. The shaded part of the spectrum was not published by S2010. Note the logarithmic scale of the xaxis in the bottom panel. 

In the text 
Fig. E.1 NLTE departure coefficients b_{2} and b_{3} (see Sect. 3.3) calculated with RH (Rutten & Uitenbroek 2012, Fig. 10) and NESSY for the FCHHTB solar atmosphere model (Fontenla et al. 2009). NOLBKG curve is the default calculation with the NLTE block of NESSY (see Fig. 8, left panel) without accounting for the effect of line blanketing. LBKG curve is the output of NESSY NLTE block taking into account the effect of blanketing by all lines within 100−1000 nm wavelength range. LOPA is the “low opacity” curve not taking into account the effect of line blanketing and corresponding to V_{dop} = 2000 km s^{1} when values are 400 times larger than those of for km s^{1}. 

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.