© 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 non-LTE (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 non-linear algebraic equations coupled with an integro-partial-differential 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 ALI-scheme. 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 co-moving 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 ALI-scheme has proved to be inefficient. In this paper we present an update of ALO in NESSY and show that the ALI-scheme 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 physics-based 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 NESSY-based models of solar irradiance is given in the reviews by Ermolli et al. (2013) and Thuillier et al. (2014a).

The inefficiency of the old ALI-scheme 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 non-accelerated Λ-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 so-called Λ-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 off-diagonal 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⁽ⁱ⁻¹⁾ = S(n⁽ⁱ⁻¹⁾). The iterations start with an initial approximation for the source function and Λ-operator corresponding to this approximation. At step i the populations n⁽ⁱ⁻¹⁾ come from the previous Λ-iteration and are considered constant. The statistical balance system is solved for n⁽ⁱ⁾ 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⁽ⁱ⁻¹⁾ in the first term of Eq. (5)with S⁽ⁱ⁾, but leaves S⁽ⁱ⁻¹⁾ in the second term. Regrouping terms in Eq. (5)and designating Λ_i−1S⁽ⁱ⁻¹⁾ as J⁽ⁱ⁻¹⁾ 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. Non-linearity of the statistical balance system

In ALI instead of substituting Eq. (5) Eq. (6) one substitutes Eq. (6). This makes the statistical balance system non-linear with respect to populations n⁽ⁱ⁾. This non-linearity 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 non-linear algebraic systems of equations (e.g. Koesterke et al. 1992). In NESSY the non-linearity 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 non-linearity comes from the so-called 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 non-linear statistical balance equations (Koesterke et al. 1992).

These three guidelines led us to choose the so-called local operator (Olson et al. 1986).

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 off-diagonal 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 ray-by-ray 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 ray-by-ray geometry is given in Fig. A.1.

We start with the equation interrelating the vectors of mean-intensity-like and flux-like 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 ray-by-ray form of the stationary spherical radiative transfer equation (Mihalas & Hummer 1974).

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 ALI-scheme vs. old ALI-scheme vs. no-ALI

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 FCHHT-B 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 bound-bound transitions and 84 bound-free 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 non-convergent behavior. The same holds true for He i levels because of the large excitation energy of the second helium level. The temperature in the FCHHT-B 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 non-convergent 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 no-ALI scheme (no acceleration in spectral lines at all) to converge. As discussed below convergence does not guarantee correct solution. In fact no-ALI 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 no-ALI 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 no-ALI scheme stabilizes at a wrong solution. It is typical for the non-accelerated Λ-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₂ 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 non-existent or low expansion velocity field) is capable of slowly releasing Λ-iterations from this “stabilization trap” (Hubeny 2003).

We note, that in the outermost atmospheric layers (>2.05 Mm) the CSA 1280 solution is the same as no-ALI 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 CSA-iterations (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 set-up of the code to reach the correct double-precision 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 ALI-scheme (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 FAL99-C 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 FAL99-C 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 SOLar-STellar 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 semi-empirical model of the atmosphere, abundances, atomic and molecular line-list, among others. In particular, the semi-empirical models have been adjusted to yield a good agreement of the output spectrum employing some specific radiative transfer code (SRPM in the case of FAL99-C 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

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 FCHHT-B 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 code-named 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 Wolf-Rayet 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 non-local, 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 non-local fraction of Lyα and Lyβ photons. The gradually decreasing density leads to increasing mean free path and therefore increasing non-local 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 non-local 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 so-called Hα opacity gap (Rutten & Uitenbroek 2012) in the vicinity of the temperature minimum. The decline of Λ^∗ in the almost isothermal FCHHT-B chromosphere is due to the gradual decrease of density.

5.2. Response to expansion and Doppler broadening velocity change

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

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 non-locality 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₂ and b₃ (as defined in Sect. 3.3) for the FCHHT-B 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₂ 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₂> 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₂< 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 non-Lyman transitions can “leak out”. This leakage creates radiative losses in non-Lyman 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₂> 1.

We attribute the slight decrease of b₂ 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₂ 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 re-enables collisions to drive the second level closer to LTE, thus forming the b₂ 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 ALI-scheme 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 ALI-scheme. The new ALI-scheme is based on the local operator (Olson et al. 1986). We described our derivation of it in CMF.

In the stationary case the old ALI-scheme is able to converge to the correct solution and allows to avoid the pathological behavior of the non-accelerated Λ-iterations. But in doing so it needs an unreasonably large number of iterations. We showed, that in the stationary case the new ALI-scheme 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 200020-153301. 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 Middle-atmosphere 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

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 ray-by-ray 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 mean-intensity-like 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).

We then use expression (7-207) 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 frequency-independent Λ-operator corresponding to mean intensity averaged over the line profile for exactly this mean intensity enters the statistical balance system in the case of bound-bound 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 Newton-like 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 bound-bound 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 off-diagonal 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 multi-variable 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 non-diagonal ones and the lower elements (B.3) normally depend only on the electron concentration.

For the regular Λ-iterations the non-linearity 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 non-linearity 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

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 ALI-scheme 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 free-free emission in the spectral synthesis block of the code;

• implementation of new Gaunt factors for hydrogen free-free 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 ALI-scheme 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 free-free 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.

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 = ν₀V_dop/c is the Doppler width, ν₀ is the central line frequency, c is the speed of light and x_k = (ν_k−ν₀)/Δν_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₂ and b₃ hydrogen departure coefficients: NESSY vs. RH

The coefficients b₂ and b₃ for the FCHHT-B 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₂ and b₃ calculated with NESSY.

Our interpretation of the origin of the b₂ peak at the temperature minimum is the same as that of Rutten & Uitenbroek (2012). However, our b₂ 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 Van-der-Waals 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₂ 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₂ value in the lower chromosphere (where the FCHHT-B temperature profile has a plateau). While b₂ equals to one according to Rutten & Uitenbroek (2012), our calculations yield a significant deviations from LTE due to the photon loss in non-Lyman 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₂ 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₂ 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₂ value to unity.

All Figures