© V. Bommier 2020

1. Introduction

This paper is the last in a series aimed to develop a theory capable of partial redistribution (PRD) and line profiles in the atomic density matrix formalism well adapted to describing the polarized atom in nonlocal thermodynamical equilibrium (NLTE) conditions. The aim of this formalism is to write down and solve the statistical equilibrium equations for the atomic density matrix elements, which include sublevel populations but also coherences (or phase relationship) between sublevels, responsible in particular for the Hanle effect. The atom is submitted to an incident radiation and to collisions, which both enter the transition rates between density matrix elements. Once these elements have been computed by solving the statistical equilibrium equations, the radiative transfer equation coefficients may be derived, which are absorption matrix and emissivity for the Stokes parameters (see, e.g., Eq. (6.85) of Landi Degl’Innocenti & Landolfi 2004), and the radiative transfer equation may then be integrated along each line of sight in the medium in order to recompute the radiation incident upon the atom in an iterative scheme.

Partial redistribution and line profiles are introduced by repelling step by step the Markov approximation in the atom–radiation interaction hamiltonian, and by synthesizing a practicable final equation from the contribution of each step. This enables relationships between incoming and outgoing photons as two-step coupled processes and beyond. The principles of the calculation are each described in Bommier (1997a). This latter publication also introduces line profiles in the equations. The profile appears as an infinite but convergent and well-known series as a function of the development, which may then be summed up, leading to the final equations, which are then nonperturbative. Although the calculation principles are more general, the final equations given in Bommier (1997a) are for a two-level atom with an unpolarized lower level. Bommier (1997b) incorporates the magnetic field effect, of arbitrary strength, and derives the partial redistribution matrix for polarized radiation in the presence of an arbitrary magnetic field.

The following paper of the series, Bommier (2016a), presents the multilevel equations able to describe the NLTE problem for a multilevel atom embedded in a magnetic field. The statistical equilibrium equations are given for the multilevel atom (multiterm as well), which have to be resolved. This enables lower-level atomic polarization for example. The advantage conferred by this process is that PRD is included in the basic formalism itself in a self-consistent manner, whereas, usually, complete redistribution (CRD) is first solved in a multilevel scheme, and PRD is then added by considering the levels pair by pair (Uitenbroek 1989, 2001). In Bommier (1997a, 2016a), extension from CRD to PRD is simply accounted for via a new fourth-order term that enters the emissivity and is comprised of a product of three profiles able to include more than two levels. An interesting feature appears in this very general formalism: the atomic density matrix statistical equilibrium has to be resolved for each atomic velocity class, and not for an averaged atom. This additionally and simultaneously enables velocity redistribution to be accounted for.

As an application, the so-called two-term redistribution matrix, which is the redistribution matrix for an atom with fine and/or hyperfine structure responsible for a series of lines connecting a lower and an upper term, as in for instance, the Na I D line pair, or the Mg II h and k lines. This is indeed a multilevel problem and the solution was enabled following Bommier (2016a). This two-term redistribution matrix is given in Bommier (2017), with a corrigendum in Bommier (2018, thanks to Ernest Alsina Ballester).

An application code of the multilevel theory is presented in Bommier (2016b), and is devoted to computing the theoretical linear polarization profile (“second solar spectrum”) of the Na I D lines observed close to the solar limb, where linear polarization is formed by radiative scattering. Observations of this profile are provided, for example, in Bommier & Molodij (2002, Figs. 2 and 3) and in Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1). The D₂ linear polarization profile interestingly displays a line center sensitive to the Hanle effect, and two far wings probably insensitive to the magnetic field and then able to serve as a local reference for the zero-field polarization.

The D₁ line polarization profile is a subject of debate: a net linear polarization is visible in Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1), whereas there is no net linear polarization in Bommier & Molodij (2002, Figs. 2 and 3), or in Trujillo Bueno et al. (2001, Figs. 2 and 3), or in Gandorfer (2000). Observations by Malherbe et al. (2007, Fig. 19) also display a net linear polarization in Na I D₁. However, its shape is different from that of Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1). This net linear polarization is surprising from a theoretical perspective, because the upper term of the D₁ line is 3²P_1/2, with J = 1/2 leading to unpolarizable term and line, even if there is hyperfine structure as in Na I. The hyperfine structure acts as a depolarizing mechanism on an already zero polarization, and therefore the net polarization should remain zero. However, this does not rule out a spectral shape of this polarization, and investigating this shape was one of the motivations of the present theoretical computation.

Because of the hyperfine structure, the modeling of the Na I D line polarization is a fully multilevel problem requiring the formalism presented in Bommier (2016a). We developed this application in Bommier (2016b). Nevertheless, the obtained profiles were very far from agreeing with observations (see Fig. 5). As the numerical code and methods are fully described in Bommier (2016b), we do not repeat them here. The present paper is devoted to concluding on this subject by presenting significantly improved results after an approximation described in Sect. 2. The calculation convergence is detailed in Sects. 3 and 4 is a concluding discussion about the remaining gap between the theoretical and observed profiles.

In our calculations, collisional broadening and collisional transitions are fully accounted for, from computations by Roueff (1974) for line broadening and Kerkeni & Bommier (2002) for transitions both due to collisions with neutral hydrogen atoms, and by applying the semi-classical method of Sahal-Brechot (1969a,b) for collisions with electrons and protons. In particular, collisional transitions between the upper J = 1/2 and J = 3/2 are fully accounted for, which remains impossible in the two-term redistribution matrix approach, as discussed in Bommier (2017).

2. Results

The new results we present in this paper were obtained after some error corrections (in summations). Above all, a numerical approximation was carried out which is able to free the results from the unacceptable dispersion shape they show in Fig. 5 of Bommier (2016b). The dispersion shape is due to the imaginary part of the profiles, which enter equations as complex quantities as visible in Eq. (63) of Bommier (1997a). As integrals are all numerical in the code, we attributed these certainly false dispersion profiles to insufficient accuracy of the integrals. As increasing the number of points and weights was impractical due to the long computing time, we forced the convergence by setting the imaginary part of the final profile in the emissivity and absorption coefficient to zero. As this appeared to bring about no effect on the fourth-order contribution responsible for the PRD, we finally did not apply it to this term. This was efficiently applied to the second-order contributions.

The emerging linear polarization is plotted in Fig. 1 together with the line intensity. The limb distance is taken at 4.1 arcsec (μ = 0.092), as in the observation by Bommier & Molodij (2002, Figs. 2 and 3). As discussed in Bommier (2016b), the atmosphere model is limited to the temperature minimum, as in the HOLMUL model (Holweger & Mueller 1974). This avoids a central bump in the line center intensity profile, which results from the temperature rise in the chromosphere, as visible in Fig. 12 of Bommier (2016b). The Na I D lines are formed in the low chromosphere. The line center bump does not exist in observations. With this HOLMUL approximation, already applied in Bommier (2016b), the computed intensity profile is in better agreement with observations.

Fig. 1. Intensity (I) and linear polarization (Q/I) theoretical spectrum of the Na I D lines as observed 4.1 arcsec inside the solar limb. CRD: Assuming complete redistribution of radiation (second-order calculation). PRD: including partial redistribution of radiation (infinite-order calculation).

The computed polarization profile is now in much better agreement with observations. However, the agreement remains qualitative. A small bump is visible in the blue wing of Na I D₂, similar but much smaller than the far wings visible in observations.

The shape of the Na I D₁ polarization profile is very similar to that of observations, at least that of Bommier & Molodij (2002, Figs. 2 and 3). This complex profile is comprised of two parts: central sharp peaks and broad wings. The amplitude of the linear polarization at the sharp peaks is significantly weaker than in the observations. However, the amplitude of the broad wings is very similar to that of the observations. Interestingly, the computations by Belluzzi & Trujillo Bueno (2013, Fig. 4) are of a similar very low polarization level for their sharp peaks. However, the broad wings of their theoretical profile are very different from those of our calculations and those of the observations. The computed Na I D₁ polarization profile seems perfectly antisymmetrical and free from any net linear polarization, which is not surprising because the Na I D₁ line is globally unpolarizable.

In Fig. 1 the linear polarization in Na I D₁ may appear weaker than in observations. However, this may be assigned to a perspective effect due to the strength of the linear polarization peak in Na I D₂, which is 0.4% in our calculations. We note that this strength is rather variable among observations, being 0.5% in Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1), of 0.35% for a different date and place of observations by Stenflo et al. (2000b, Fig. 1 at μ = 0.1), and of 0.3% in the observation by Bommier & Molodij (2002, Fig. 2). The depolarizing effect of the magnetic field may play a role in this line center polarization sensitive to the Hanle effect. Our calculations do not include any magnetic field effect in a first step.

3. Convergence

This section is devoted to demonstrating that the results are essentially converged in terms of polarization degree Q/I. Indeed, the calculation method outlined at the beginning of this paper is of the lambda-iteration type, which is generally considered as poorly convergent. The implementation of acceleration methods was already discussed in Bommier (2016b). Below, we summarize and complement this discussion, and then present the results of the present calculations.

3.1. Convergence acceleration

Ng acceleration was attempted but was found inefficient, as reported in Bommier (2016b). Figure 3 of that paper shows that the emergent polarization is modified–sometimes increased, other times decreased–without any consistent convergence-improving trend. Indeed, this method is nothing more than a linear extrapolation, when the solution of the statistical equilibrium equations is linear in the two-level case only, which is not the present case.

The preconditioning would probably be inefficient for our density matrix comprised of 640 elements, in addition decomposed over 48 velocities. As demonstrated in Sahal-Bréchot et al. (1998) and Bommier (2016a), the correct procedure for summing over atomic velocities is to solve the statistical equilibirum equations for each velocity class, and to integrate the emerging radiation over atomic velocities. The preconditioning would then be applied to a few density matrix elements of 1 velocity over 48. Weak efficiency is then expected. Alternatively, if different velocities become coupled by preconditioning, the full system becomes 640 × 48 = 30720 elements in width, which would be huge. In addition, the complexity of the polarized atom formalism (Racah algebra), which is different for statistical equations (see Eqs. (9)−(13) of Bommier 2016a) and radiative transfer equation coefficients (see Eqs. (14)−(20) of Bommier 2016a) makes the preconditioning implementation such complex task that we renounced. We note that previous works successfully applied preconditioning up to ten levels only, without any velocity splitting up and with angle-averaged PRD for two by two levels (Uitenbroek 2001). Given all these arguments, we decided not to implement any preconditioning, owing to its complexity in the polarized radiation and density matrix formalism, together with little hope of improvement because of the large number of matrix elements expanded over atomic velocities.

We studied the implementation of the Feautrier method in our formalism. The second-order Feautrier difference equations are given in Bommier (2019) in the case of polarized radiation and density matrix formalism. However, in the multilevel case, the solution has to remain iterative, and was found to be overly time- and memory consuming to be applied to the present computation. In simpler trials, ignoring PRD or hyperfine structure, the convergence was found not to be significantly improved with respect to the convergence of the lambda-iteration method we applied here, which includes the short characteristics method renewed by Ibgui et al. (2013) for radiative transfer equation integration.

3.2. Convergence results

We study the outgoing radiation polarization degree Q/I, which is the result of our calculation. The behavior of this quantity as a function of the iteration step is represented in Fig. 3, in the line center and in the line wings (at 5888.4 Å for Na I D₂ and 5894.9 Å for Na I D₁). When polarization quickly converges in the line center, this is not the case in the line wings.

In order to demonstrate that convergence exists in line wings, we compare the behavior of our series with the behavior of the Riemann series:

(1)

When α = 1, this is the harmonic series, which does not converge. However, when α >  1 strictly, this is a Riemann series, which is absolutely convergent.

We therefore compare the step increment of our series to the inverse of the step number 1/N in order to determine whether or not α is larger than unity. The top panel of Fig. 3 shows the ratio of this increment Δ(Q/I) to 1/N, as a function of the iteration step number N. For sufficiently large step numbers (N >  60), Δ(Q/I) can be seen to decrease more quickly than 1/N, which means that Δ(Q/I) behaves as 1/N^α with α >  1, meaning that our series is of the Riemann type, and is therefore absolutely convergent even in the line wings.

In order to evaluate α in the case of our series, in the bottom panel of Fig. 3 we plot the logarithm of the ratio, and we see that this logarithm decreases linearly for both wings, which determines α = 1.00058 for the D₁ wing and α = 1.0012 for the D₂ wing, which are both larger than unity. As can be seen in Fig. 3, the slope of the linear interpolation is found to be significantly larger than the numerical noise. The series is therefore absolutely convergent. However, the behavior of the series and of its convergence remains unverified for the next iteration steps.

Figure 2 shows that for the 166 iteration steps we performed, the outgoing polarization Q/I is not far from convergence in the line wings, but is not yet completely reached. These computations are rather time consuming, and we spent 26 Mh of computations in the TURING machine at IDRIS¹ to obtain these 166 iteration steps; we therefore estimate the above convergence to be sufficient. Convergence is however not fully reached in the line wings within the allocated computing time. Development of efficient acceleration methods would be desirable for future work.

Fig. 2. Convergence of the linear polarization rate Q/I, for the D2 line (left) and for the D1 line (right) at line center (upper row) and in the line far wing (lower row).

Fig. 3. Convergence of the emerging linear polarization Q/I: evolution of the step increment Δ(Q/I) as a function of the step number N. The aim is to show that the step increment decreases more rapidly than the inverse of the step number 1/N (as for a Riemann series). In this respect, their ratio is the upper plot and the logarithm of their ratio is the lower plot, which also includes linear fits and their equations.

The two terms of the Q/I ratio, namely the Stokes parameters Q and I, each nevertheless converge much less quickly, in such a way that they are not at all converged within the 166 iteration steps. However, the behavior of the ratio Q/I appears to be different, because some simplifications certainly occur when forming the ratio.

4. Concluding discussion

Although the calculation seems reasonably converged in the line wings, the agreement between the observed and computed profiles remains only qualitative. As described below, we investigated the different possible causes of this difference.

We questioned the accuracy of collision rate calculations, in particular those with hydrogen neutral atoms that determine the coherent scattering weight in PRD, and therefore in the far wings. We artificially varied these rates, but the Na I D₂ wing profile was not improved with respect to observations.

We questioned the model electron density, whose determination is very indirect from observations, and in light of the recent investigation by Bommier (2020), who suggests that the electron density at the solar surface would be much higher than in present models, and similar to the neutral hydrogen atom density. We artificially varied the model electron density, but again the Na I D₂ wing profile was not improved with respect to observations.

We questioned the accuracy of our model of Rayleigh scattering on neutral hydrogen atoms, which is responsible for the continuum below spectral lines, which is of more relative importance in their wings. Our model is taken from the MALIP code for computing the Stokes parameter profiles of magnetoactive Fraunhofer lines by Landi Degl’Innocenti (1976), and we did not try to modify it.

We questioned the very rough HOLMUL approximation (Holweger & Mueller 1974) we applied to our atmosphere model. This approximation consists in not considering atmosphere higher than the temperature minimum. The temperature reversal there may be responsible for a central bump in the line profile, following a well-known result for a two-level atom in LTE. As our model considers only the lower and upper levels of the Na I D lines (with all their sublevels and coherences), this is probably the reason why this bump appeared in our calculations, as visible in Fig. 12 of Bommier (2016b), which is not the case in observations. One solution for this problem could lie in considering several upper levels of Na I, as done by Bruls et al. (1992). Leenaarts et al. (2010) suggested modeling the contribution of several higher upper levels by grouping them into a single artificial upper level. We tried to implement this suggestion but it was unsuccessful because we treated the radiative transfer only in the Na I D lines to remain within computing time restrictions.

We question the 1D plane parallel character of the usual atmosphere model above the temperature minimum region. In hydrogen Hα images, this part of the atmosphere contrarily seems highly inhomogeneous and structured into fibrils. Also inspired by the first ionization potential (FIP) effect observed higher in the corona, we tried an increase of the atmosphere model electron density in this part of the atmopshere, as if the electrons were grouped into inhomogeneous structures along the line of sight above the temperature minimum region. This was unsuccessful.

We finally note that the Na I D line intensity profiles are broader in our computation results than in all the observations. The central peak of the Na I D₂ line polarization is accordingly broader. This may be the reason why the Na I D₂ wing bump is much less visible in our results than in observations. In relation, we question the fact that we neglected the inhomogeneous upper part of the atmosphere. In other words, we assign this discrepancy to the inadequacy of the plane parallel atmosphere model to describe the chromosphere, where the Na I D lines are formed. We think that this is the main limitation of our computation. As all integrals over frequencies and velocities, and their directions, are numerical, the calculation could probably be adapted to more realistic and higher dimensional model atmospheres.

Nevertheless, our approach is the first self-consistent treatment of redistribution in NLTE conditions. The lower level atomic polarization is fully accounted for, although found to be negligible, unlike what was assumed by Landi degl’Innocenti (1998); it is destroyed by collisions that are fully accounted for. However, in the case of a line with hyperfine structure, as in the Na I D lines, the calculation is time consuming.

This work was motivated by the existence of a net linear polarization in Na I D₁, as observed by Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1), as well as by Stenflo et al. (2000b, Fig. 1). On the contrary, Bommier & Molodij (2002, Figs. 2 and 3), Trujillo Bueno et al. (2001, Figs. 2 and 3), and Gandorfer (2000) did not observe any such net linear polarization in Na I D₁. As expected from the Na I D₁ upper level kinetic momentum quantum number J = 1/2, our computation does not reveal any net linear polarization in Na I D₁. We note that Bommier & Molodij (2002) and Trujillo Bueno et al. (2001) both observed with the THÉMIS² telescope, which was polarization-free at that time, and that both works did not apply the same data-reduction code. Gandorfer (2000) observed at IRSOL³, where the telescope is polarization-free around equinox, whereas Stenflo & Keller (1997) and Stenflo et al. (2000a,b) observed with the Zürich Imaging Polarimeter ZIMPOL mounted on the McMath-Pierce telescope at Kitt Peak, which has large and varying instrumental polarization (Stenflo et al. 2000b, p. 782). Instrumental polarization correction was installed with a tilting glass plate. The compensation level is about 0.1% (Stenflo et al. 2000b, p. 782), which is also the order of magnitude of the Na I D₁ linear polarization. Therefore, the compensation may be insufficient given the low linear polarization of the Na I D lines, and the net linear polarization observed in Na I D₁ with this instrument only, lacking confirmation by theoretical calculations, could be an artefact of instrumental polarization.

In addition, as already stated at the end of Sect. 2, the linear polarization observed by Stenflo & Keller (1997, Fig. 2) and Stenflo et al. (2000a, Fig. 1) in the Na I D₂ line center is also higher by 0.1% than the present calculation result. This cannot be assigned to a magnetic field effect not accounted for in our calculations because the Hanle effect acts as a depolarizing effect in the second solar spectrum. The difference in polarization degree between our calculation result and these observations is about 0.1%, which is comparable to the instrumental compensation level. The linear polarization in Na I D₁ similarly shows a narrow peak of 0.1% linear polarization at line center, which is not present in the calculation result. The 0.1% excess polarization observed in both line centers could therefore be an artefact of instrumental polarization.

Acknowledgments

For set up, this work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche (France). For exploitation, this work was granted access to the HPC resources of IDRIS under the allocation 2017-A0020407205 made by GENCI (France).

References

All Figures

	Fig. 1. Intensity (I) and linear polarization (Q/I) theoretical spectrum of the Na I D lines as observed 4.1 arcsec inside the solar limb. CRD: Assuming complete redistribution of radiation (second-order calculation). PRD: including partial redistribution of radiation (infinite-order calculation).
In the text

	Fig. 2. Convergence of the linear polarization rate Q/I, for the D2 line (left) and for the D1 line (right) at line center (upper row) and in the line far wing (lower row).
In the text

Fig. 3. Convergence of the emerging linear polarization Q/I: evolution of the step increment Δ(Q/I) as a function of the step number N. The aim is to show that the step increment decreases more rapidly than the inverse of the step number 1/N (as for a Riemann series). In this respect, their ratio is the upper plot and the logarithm of their ratio is the lower plot, which also includes linear fits and their equations.

In the text