© ESO, 2017

1. Introduction

In the literature, there are several proposed empirical models specific to sunspot umbrae (Avrett 1981; Maltby et al. 1986; Obridko & Staude 1988; Collados et al. 1994; Severino et al. 1994; Socas-Navarro 2007; Fontenla et al. 2009). Most of these models have been obtained empirically with observational data taken at several wavelengths (UV, visible, etc.). However, although modeling the temperature and vector magnetic magnetic field stratification in the photospheric part of sunspot umbrae is relatively easy, this is not the case for the chromosphere and layers above. In these regions, most observed spectral lines suffer from non-local thermodynamic equilibrium effects (NTLE), which makes it difficult to determine the height variation of several thermodynamic parameters accurately. For this reason, while the different available models look similar to each other in the photosphere these models significantly differ at higher layers, especially at the chromosphere (see Fig. 1). To properly model the thermodynamic parameters at these layers, spectral lines sampling of those heights should be included in the modeling process. This particular problem has been already addressed (see, e.g., review of Solanki 2003). The most recent model we find in the literature is that from Socas-Navarro (2007). These authors inferred a sunspot umbral model from the analysis of several chromospheric Ca ii  and photospheric Fe I  lines. They used an inversion code that is able to synthesize spectral lines in NLTE conditions. Other umbral models do not take into account the chromosphere (Collados et al. 1994).

Fig. 1 Variation of the temperature with atmospheric height corresponding to various atmospheric models available in the literature. In particular, the E model of Avrett (1981), the models by Maltby et al. (1986) and Fontenla et al. (2009) that cover from the photosphere to the chromosphere, and the model of Collados et al. (1994; cool umbra and hot umbra) that only covers the photosphere. The optical depth scale has been calculated using the Maltby et al. (1986) atmospheric model.

Iwai et al. (2016; see also Iwai & Shimojo 2015) have found, using the Nobeyama Radioheliograph (NoRH; Nakajima et al. 1994), that the observed umbral brightness temperature is the same as that of the quiet Sun at 34 GHz (corresponding to an 8.8 mm range). Their findings are inconsistent with current model predictions, which suggest that sunspot umbrae should be brighter than the quiet Sun at the millimeter spectral range. Indeed, none of the umbral models available to the solar community are able to explain these radio observations (Iwai et al. 2016). From the results of Iwai et al. it seems that the temperature should start increasing at a lower geometrical height, so that the minimum temperature region should be located lower in the atmosphere when above sunspot umbrae. As we stated above, constructing sunspot umbrae atmospheric models for the upper photosphere and low chromosphere is not straightforward since it entails the inclusion of NLTE spectral lines in the modeling process. So far, only Socas-Navarro (2007) have included chromospheric lines. In this paper, we present a sunspot observation taken in the Si I 1082.7 nm line ($\mathrm{4}{\mathrm{s}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{2}}^{\mathrm{o}}\mathrm{-}\mathrm{4}{\mathrm{p}}^{\mathrm{3}}{\mathrm{P}}_{\mathrm{2}}$, with Landé factor 1.5, excitation potential 4.95 eV and log  gf = 0.363¹) at high spatial resolution (≈0.̋3) with the GRegor Infrared Spectrograph (GRIS; Collados et al. 2012). Remarkably, the core of the Si I 1082.7 nm line is mainly formed in NLTE conditions (Bard & Carlsson 2008). The novelty of these observations is that, once the stray-light effects, including the telescope, atmospheric seeing, and scattered light, are removed from the observations, the Si I 1082.7 nm line appears in emission at the umbral core. The analysis of these observations, as we detail throughout this paper, suggests that the minimum temperature region should be located at lower heights, at least for this sunspot, with a temperature hump located at about log τ = −2. Should this result be extended to average sunspot umbrae models, it could provide a natural explanation of the controversy between the present empirical models and the radio observations at the millimeter range. Moreover, it could also settle a new paradigm regarding umbral stratification. In Sect. 2 we describe the observations and explain the deconvolution process. In Sect. 3 we detail the Si I 1082.7 nm emission profiles and discuss the results of the deconvolution process, and in Sect. 4 we cover the inversion of the data under NLTE conditions. Finally, in Sect. 5 we summarize the results and highlight the main conclusions and possible impacts.

2. Sunspot observations and PSF correction

Fig. 2 Intensity maps of the observed NOAA 12096 sunspot corresponding to the continuum (left panels) and to the Si I 1082.7 nm line core (center panels). The top panels represent the original data and the bottom panels the reconstructed data. The right-most panels show a zoom over the original (top) and reconstructed (bottom) sunspot umbra in the continuum (see the square in the top left panel). The contours outline the region where the Si I 1082.7 nm continuum is below 25% of the average quiet Sun.

The sunspot observations were taken on 27 June 2014 with the GRIS spectrograph of the 1.5 m German GREGOR telescope (Schmidt et al. 2012) at the Observatorio del Teide (Tenerife, Spain). The sunspot (NOAA 12096) was located close to the disk center (x = −72,5′′, y = 148′′). The GRIS spectograph is equipped with an infrared camera (TIP II; Collados et al. 2007) and a polarization modulation package. Thus, the four Stokes parameters were measured with TIP-II around the near-infrared spectral region containing the He I 1083.0 nm multiplet and the Si I 1082.7 nm photospheric line and an atmospheric water vapour line at 1083.21 nm. The size of the scanned region is 32′′× 63.̋3. The four Stokes parameters were recorded by moving the slit (0.̋135 wide) in 240 step positions with a step size of 0.̋135. The effective pixel size of the data is 0.̋135. The total effective exposure time per slit position was about 1.5 s, so that the full scan lasted 6 min from 10:06 UT until 10:12 UT.

The standard data reduction software was applied to the GRIS data, which include dark current and flat-field correction, polarimetric calibration, and cross-talk removal (see Collados et al., in prep.). The wavelength calibration was carried out by fitting the observed spectral region to the McMath-Pierce Fourier Transform Spectrometer (FTS) spectrum (Livingston & Wallace 1991; Wallace et al. 1993), taken the Si I 1082.7 nm line of the average quiet Sun and the water vapour line at 1083.21 nm of the atlas as reference wavelengths. The resulting spectral sampling is 54.3 mÅ/pixel, after averaging three pixels in the wavelength dimension. This operation also allowed us to correct the continuum level of the data; that is, we removed the low frequency oscillation of the continuum. The spectral resolution of GRIS at 1083 nm should be around 120 mÅ accordingly to Joshi et al. (2016). More information about these procedures can be found in Franz et al. (2016) and Borrero et al. (2016).

The seeing conditions during the observations were good. The GREGOR Adaptive Optics System (GAOS; Berkefeld et al. 2012) was running continuously. The alto-azimuthal mounting of the telescope introduced an image rotation² smaller than 0.5 [°/ min] at 10 UT (Franz et al. 2016), so that this effect is almost negligible in our data set. Finally, we removed an 8% of spectrograph scattered light from the observations with the procedure described in Franz et al. (2016) and Borrero et al. (2016).

The observations are affected by the telescope optics, the Earth’s atmosphere (seeing), and scattered light. This stray-light PSF(r) can be modeled by the sum of two Gaussian components (Felipe et al. 2016) as follows: $\mathrm{PSF}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{\alpha }{\mathit{G}}_{\mathrm{1}}\mathrm{\left(}\mathit{r,}{\mathit{\sigma }}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}{\mathit{G}}_{\mathrm{2}}\mathrm{\left(}\mathit{r,}{\mathit{\sigma }}_{\mathrm{2}}\mathrm{\right)}\mathit{,}$(1)where G₁(r,σ₁) and G₂(r,σ₂) are two Gaussian functions representing the spatial resolution of the telescope and atmospheric seeing and the light scattered over large angles. The values σ₁ and σ₂ stand for the rms width of both contributions while α is the relative contribution of the Gaussians, and r represents the radial distance r² = x² + y². Both contributions were determined experimentally with observations of the Mercury transit in front of the Sun on 2016 May 9 and GREGOR adaptive optics. Weather conditions were bad during the Mercury transit but the acquired data were good enough to allow the determination of the instrumental stray-light contribution to the GREGOR PSF, which is σ₂ = 5.030 ± 0.001′′and β = (1−α) = 42%. Once the contribution of the second Gaussian is fixed, we used the information provided by the GAOS system to compute the contribution of the telescope and the Earth’s atmosphere to the GREGOR PSF during our observations. In our case, the estimated Fried parameter at 500 nm was r₀ = 15.4 cm with an rms variation of ± 0.6 cm, which corresponds to a σ₁ = 0.24′′, in the near-infrared at 1083 nm. In reality, the PSF varies continuously during the 6 min observation period, mostly due to the changing conditions of the Earth’s atmosphere (turbulence). For the calculations above we took the average Fried parameter over the 6 min period. The rms variation of the Fried parameter was small though, as it was approximately ± 0.6 cm during the observations.

Fig. 3 Variation of the original (blue) and reconstructed Stokes I/IQS and V/IQS profiles as we move from the inner penumbral boundary (first column) to the center of the umbra (last column). The Stokes I profiles in black correspond to the deconvolved profiles taking into account 42% of stray light while for the red profiles we only accounted for 30% of the stray light. For simplicity, we only represent five pixels. The location of the pixels can be found in the right panels of Fig. 2. The first profile corresponds to the lowest point.

We used the principal component analysis (PCA) deconvolution method described by Ruiz Cobo & Asensio Ramos (2013) and Quintero Noda et al. (2015) to remove the stray light. This deconvolution technique has already been applied to GRIS data by Borrero et al. (2016), who describe how the method works in detail. In our case, we limited the information to the first 12 PCA coefficients and performed only 20 iterations within the deconvolution process. We applied the PCA deconvolution method to the spectral region containing only the Si I 1082.7 nm line. It is important to note that the PCA deconvolution method is applied to “solar images” constructed from the slit images (240 in our case). In this sense, the solar scene and the PSF vary along the scan direction, where this issue is one of the main drawbacks of image deconvolution techniques applied to spectroscopic data. However, we find that the most important effect comes from the contribution of the second Gaussian of the PSF (stray light), which we assume to be constant during the observations. The amount of stray-light contribution coming from the second Gaussian depends on the solar scene as well, but we can assume that the evolution of the solar scene along the scanning direction has negligible impact on the stray-light contribution because this contribution comes from wide angles. Finally, the deconvolution process is applied independently to the four Stokes parameters.

Figure 2 shows continuum Stokes I/I_QS and Si I 1082.7 nm line core maps of the observed sunspot before and after applying the PCA deconvolution, where I_QS is the average quiet Sun intensity at the continuum. It can be seen how the contrast and image sharpness has increased considerably after the deconvolution of the data, especially at the line core. The effect in the sunspot umbra is also evident; for instance, the intensity at the umbra has clearly diminished after the deconvolution. However, as we see in detail in the next section, the most conspicuous effect can be found in the spectra recorded at the sunspot umbra. In the original data, the Si I 1082.7 nm line was always in absorption, while it becomes in emission after the deconvolution. The contours in Fig. 2 outline the area where the sunspot umbral intensity is lower than 25% (compared to the quiet Sun). Within that area, the core of the Si I 1082.7 nm line appears in emission.

3. Si I 1082.7 nm profiles in emission

Figure 3 shows how the shape of Stokes I/I_QS and V/I_QS changes from the inner penumbral boundary to the umbral core for the estimated 42% of stray light and for 30% of stray light for comparison. In the original data, Stokes I/I_QS never goes below 40% of the intensity (relative to the quiet Sun average intensity). Moreover, the profiles are rather regular and show three tiny peaks at the line core, basically as a result of the Zeeman splitting under the presence of a strong magnetic field; the left and right absorption peaks correspond to the σ components while the central absorption peak correspond to the π component. The stray-light contamination contributes to the central absorption peak, since the stray light has, in general, the shape of Stokes I in the absence of a magnetic field. The π component is absent if the vector magnetic field is strictly pointing toward or away from the observer. Similarly, Stokes V/I_QS shows the prototypical anti-symmetrical shape with two lobes of opposite sign. However, the deconvolved profiles show a completely different scenario. In Stokes I, the intensity has dropped considerably, reaching values below 20% at the umbral core where the profile clearly shows the σ components. The line core intensity is much greater than the intensity at the wavelength location of the σ components. As we see later, we ascribe such intensity excess to line emission and refer to it as emission peak. In the profiles located at the umbral core (rightmost profile in Fig. 3), the emission peak is clearly discernible, flanked by two tiny absorptions. Regarding Stokes V/I_QS, their amplitudes with respect to the quiet Sun decrease only slightly after the deconvolution. Interestingly, the Stokes V/I_QS profiles always show a reversal in the zero-crossing wavelength in all cases although it is much more prominent at the umbral core. The intensity reversal is not recovered when the straylight is only reduced to 30%. Regarding Stokes Q/I_QS and Stokes U/I_QS (not shown), their amplitude is very small, i.e., below 1% for the former and negligible (within the noise level) for the latter. This discards the magneto-optical effects as an explanation for the Stokes V/I_QS zero-crossing reversal.

Fig. 4 Variation of umbra continuum intensity with wavelength (right panel) for three different models (left): the FALS, the cool model of Collados et al. (1994), and a super-cool model obtained from the cool model. The vertical line shows 1083 nm wavelength while horizontal lines pinpoint the corresponding continuum intensity values at 1083 nm for the three models.

Overall, the effects of the deconvolution process in the Stokes profiles (especially for Stokes I/I_QS) are rather dramatic. For this reason, we evaluated whether the deconvolution worked properly by looking at possible side effects ascribed to the PCA deconvolution strategy or to an inaccurate determination of the PSF of the telescope. The first thing to check is the recovered continuum intensity levels, which have dropped considerably with respect to the original values. In the darkest part of the umbra, the continuum intensity changes from the 42% in the non-deconvolved data to about 17% after the reconstruction as suggested by the Stokes I/I_QS shown in the rightmost panel of Fig. 3. However, on average, the intensity variation in the darkest area of the umbra (outlined in Fig. 2, right panel) goes from 44% of the continuum in the non-reconstructed data to about 20%. If we take into account the whole umbra, which is defined as the region in the reconstructed data that has a continuum intensity below 40%, the continuum intensity decreases from 50% to about 30% on average. Current models for a sunspot umbra of about 2.5 Mm radius³ suggest about 25% of continuum intensity (following Kiess et al. 2014, see their Fig. 5) for visible wavelengths. The continuum intensity also varies with wavelength, which is typically larger for longer wavelengths. For instance, Maltby et al. (1986) suggest continuum values of about 22% for 1 μm (see their Fig. 2). The continuum intensities recovered in the GREGOR data are slightly lower than the predictions for the darkest part of the umbra. To shed light into this issue we show, in Fig. 4, three different atmospheric models, i.e., the FALS, the cool model of Collados et al. (1994), and a super-cool model, along with the variation of the corresponding intensity values at the continuum with wavelength. The super-cool model was tailored from the cool model to get continuum values of about 17% in the 1083 nm region. To get such continuum intensity values in the infrared, we needed to shift the temperature down about 400 K. Thus, the model would be rather cool. Interestingly, the intensity values that would be recovered with this model for visible wavelengths (around 630 nm) would be about 7%, which is in line with the intensity values reported in Hinode observations with deconvolved data (Ruiz Cobo & Asensio Ramos 2013). Remarkably, the cool model of Collados et al. (1994) already predicts intensity values of about 30%, which is in line with the (averaged) intensity values for the deconvolved sunspot umbra of around 30%.

Fig. 5 Variation of Stokes I/IQS and V/IQS profiles with the amount of scattered light β from 0% to 42%. The profile corresponds to that shown in the fifth panel on Fig. 3.

Therefore, we believe that a β = 42% of scattered light for the GREGOR telescope seams reasonable. Remarkably, if we reduced the stray-light contribution to 30%, we would obtain continuum intensity values significantly higher than those predicted in the models, for the darkest part of the umbra (about 30% on average). In Fig. 5 we represent the variation of the Stokes I/I_QS and V/I_QS shapes with increasing amount of scattered light (β = 0%,10%,20%,30%,and 42%). For β = 30% and 42%, the recovered Stokes I/I_QS can be easily associated with line emission. We only recovered a Stokes V/I_QS line reversal at the zero-crossing point with values above β = 30%. For lower β values, it is not possible (by eye) to guess whether the line is in emission or if it is just the Zeeman splitting of the line. The figure also demonstrates the fact that the contribution of the second Gaussian (scattered light) to the global PSF is the most critical one. Changing the width of the first Gaussian (core of the PSF) basically affects to the spatial sharpness only.

Furthermore, when we recovered the emission feature in Stokes I/I_QS, we simultaneously recover a line reversal around the zero-crossing wavelength in Stokes V/I_QS. Since the PCA deconvolution method is applied separately to the Stokes I/I_QS and Stokes V/I_QS profiles, it could be argued that the emission feature is due to errors within the deconvolution process itself. Thus, it makes sense to check whether the recovered Stokes I/I_QS and V/I_QS are compatible. We discuss this issue in Sect. 4.2.

Another source of error may be found in the fact that the PCA deconvolution method uses a database to reconstruct the Stokes profiles. Such a database is constructed using the very same (observed) Stokes profiles. Thus, in the case of Stokes I, the database does not contain Stokes I/I_QS profiles in emission. But, as one may erroneously think, this does not limit the PCA method to deliver profiles with emission signatures since they are constructed from linear combinations of all available profiles in the database, i.e., 12 in our case. The exact same situation takes place for Stokes V/I_QS, although Stokes V/I_QS profiles with opposite polarity are available in the database. In any case, and to discard any other effects within to PCA deconvolution process (such as that previously mentioned), we repeated the data deconvolution using a completely different method known as the maximum entropy method (Hollis et al. 1992). This method is based on fast Fourier transforms, so it is rather different from PCA. The results from the deconvolution, shown in Fig. 6, are exactly the same for the sunspot umbra. The only noticeable difference is that the deconvolved profiles using the maximum entropy method are much more affected by the noise inherent to the deconvolution process, and actually an advantage of using PCA is that we can keep the noise at low levels since PCA also acts as a noise filter.

We also varied several parameters within the PCA deconvolution method concerning the number of iterations and number of profiles used to construct the database but none of these parameters affected the deconvolution results in a significant way, i.e., we still get the line core emission (see Appendix A). Interestingly, we see small variations, although at the very limit around the Stokes V/I_QS zero-crossing point. This issue is discussed in detail in Sect. 4.2.

Finally, to explore the deconvolution method further, we applied it to the red side of the observed spectral region containing the Ca I 1083.34 nm (see, e.g., Felipe et al. 2016) and the Na I 1083.48 nm lines both located at the sides of a Ti I 1083.36 nm atmospheric line and the Ca I 1083.9 nm (see, e.g., Joshi et al. 2016, 2017). The results of the deconvolution (see Appendix B for a detailed description) suggest that these lines are not in emission in the Sunspot umbra. This may be simply because these lines are weak and do not reach the upper photosphere or low chromosphere.

Fig. 6 Left panels: umbra Stokes I/IQS and V/IQS profiles reconstructed using the PCA deconvolution and maximum entropy methods. Right panels: umbra Stokes I/IQS and V/IQS profile as seen by GREGOR and the German VTT. The solid line represents the non-reconstructed data while the dotted line stands for the reconstructed data. The black line corresponds to the German VTT (simulated), red to the GREGOR telescope, and blue to the observed VTT profile (see text).

The question remains as to why these emission profiles have not been seen before, for instance, at the German Vacuum Tower Telescope (VTT) of the Observatorio del Teide. The TIP-II instrument has been operating at the VTT for many years and there are many observations of sunspots with that telescope and instrument. The main differences between the VTT and the GREGOR telescopes are the size of the primary mirror (70 cm versus 1.5 m) and that the VTT has, in principle, less scattered light, i.e., about 10% (Beck et al. 2011). We carried out the following test to see how the observed and reconstructed GREGOR umbral profiles would have been seen with the German VTT. We took the reconstructed sunspot data (as observed in GREGOR) and degraded these data with a PSF and a pixel sampling corresponding to the German VTT (assuming same seeing conditions and 10% scattered light). The results are shown in Fig. 6 (right panels). It can be seen that the shape of the simulated Stokes I profile (solid black) resembling the German VTT is very similar to the observed and reconstructed GREGOR profile (dotted red). Moreover, the line core intensity almost reaches that of the continuum level. The large difference between the non-reconstructed GREGOR profile (solid red) and the corresponding simulated VTT profile (solid black) can be directly ascribed to the different amount of scattered light (42% in GREGOR versus 10% in the VTT) and, to a lesser extent, to the spatial resolution. When we compared the deconvolved GREGOR profile with the simulated VTT profile (red and black dotted curves, respectively) the main difference between them is simply due to the smaller spatial sampling (and resolution) of the German VTT (assuming same atmospheric seeing conditions). This simplistic simulation suggests that umbral emission profiles should have been seen earlier at the VTT, without applying any deconvolution method to the data. For this reason, we reviewed previous VTT data, looking for a sunspot at disk center with similar size and good observing conditions. In particular, we made use of a data set taken on 27 September 2011 with a 0.̋32 spatial sampling. The radius of the spot was about 4.8 Mm, which is a bit larger than that presented here and was located at x = 228′′, y = 90′′. A profile taken from the umbral core of the sunspot contained in this data set can be seen in Fig. 6 (solid blue). Remarkably, the continuum level⁴ of the observed Stokes I/I_QS at the VTT is very close to the simulated Stokes I/I_QS with the GREGOR observations, i.e., about 0.3, so a wide-angle scattered light of 10% at the VTT seems to be a good lower limit. The main difference between the simulated and observed intensity profile is that the latter shows less contrast in the wavelength dimension, i.e., it is shallower. This difference could be perfectly ascribed to the spectral resolution and a non-negligible spectral scattered light, which has not been accounted for in the observed VTT profile. Regarding Stokes V/I_QS, although it has a smaller amplitude, its shape is very similar to the simulated Stokes V/I_QS around the zero-crossing wavelengths (solid black). We believe that if we apply a deconvolution on these data, we would recover an emission profile as we do with the GREGOR data. We can argue that, at the limited spatial resolution of the German VTT, the amount of scattered light is critical to the point where if the scattered light is larger than 10%, the emission signals would be significantly suppressed. To summarize, we think that the German VTT telescope was at the limit of the detection of the emission signals in the Si I 1082.7 nm line.

4. NLTE treatment of the Si I 1082.7 nm line

In this section we try to infer the physical conditions on the umbra atmosphere of the sunspot that gives rise to the Si I 1082.7 nm line core emission. To this end, we generate synthetic Si I 1082.7 nm Stokes profiles and compare these with the observed profiles, putting special emphasis on the line core emission. This forward modeling process, known as “inversion” of Stokes profiles, allows us to infer the temperature and magnetic field vector stratification in this particular sunspot umbra. The Si I 1082.7 nm line is mainly formed in the photosphere but the line is very deep and its core is dominated by NLTE (e.g., Bard & Carlsson 2008). Therefore, solving the radiative transfer equation in LTE is not sufficient for our purposes (see, e.g., Centeno et al. 2006; Kuckein et al. 2012), since the emission is at the line core. Thus, we need to solve the radiative transfer problem in NLTE. There are codes available for the solar community, such as the MULTI (Carlsson 1986) or the RH (Uitenbroek 2001) codes, but these codes are not able to perform inversions. Here we make use of NICOLE (Socas-Navarro et al. 2000, 2015), which is able to infer the atmospheric parameters through inversions of the Stokes profiles. This code has never been used to invert observations of the Si I 1082.7 nm line so the results presented in this paper concentrate on a single profile alone. A full analysis of the sunspot is presented elsewhere. The first, and more important, ingredient we need to tackle the NLTE problem is a tractable Silicon atom model. The atom model allows us to compute the level populations involved in the Si I 1082.7 nm line transition. There are several atom models available in the literature, for example, those by Bard & Carlsson (2008) or Shchukina et al. (2012). Particularly, in the first case, the authors started from a complex and complete model (as was used in the latter work) and developed a simpler model containing the main characteristics of the silicon atom (including the transition of interest). Thus, they constructed a quintessential silicon model that is computationally tractable for NLTE calculations. The use of simpler atomic models is recommended since they simplify the inversion process, which can take tens of hours in a standard computer. This also speeds up the analysis when one needs to invert a large amount of pixels.

Fig. 7 Spectral atlas of the Si i line (solid). We include the results of a LTE (dashed) and NLTE (dot-dashed) synthesis using the FALC model atmosphere.

4.1. NLTE synthesis

We have checked whether NICOLE produces meaningful results using the silicon atom model of Bard & Carlsson (2008). To this end, we computed the intensity spectrum under LTE and NLTE conditions and compared the resulting profiles (see Fig. 7). In the plot, we included the Si I 1082.7 nm profile taken from the solar atlas of Delbouille et al. (1973). To set the atomic parameters of the transition, we followed the values listed on Table A.1 of Borrero et al. (2003) including collisional broadening effects (Anstee & O’Mara 1995; Barklem & O’Mara 1997). We used the quiet Sun FALC model of Fontenla et al. (1993) to have a fair comparison with the atlas spectrum and also to be able to compare our results with those presented in Bard & Carlsson (2008). The temperature stratification of this model is plotted in Fig. 1. We kept the original microturbulence stratification of the FALC model and added a macroturbulence of 2 km s^-1 to match the typical convective broadening of the lines (Bard & Carlsson 2008, used the same value in their synthesis). Figure 7 shows that the Si I 1082.7 nm line wings are identical in both LTE and NLTE computations, as expected, while there are large differences at the line core wavelengths. In particular, the NLTE profile is much deeper than the LTE profile. A similar result is presented in Fig. 3 of Bard & Carlsson (2008). Regarding the atlas profile, we can see that the NLTE profile is deeper, which could be related to the atomic parameters used here (see the discussion on Bard & Carlsson 2008). We plan to further investigate this particular issue in a future work. These differences were also found on Bard & Carlsson (2008; see their Fig. 18). This simple proof assures us that the atom model has not been perverted during the conversion process.

In order to check if the Si i line appears in emission, we synthesized the line profile using the three semi-empirical umbral models presented in Fig. 1. None of the models provided information about the magnetic field vector so we added a constant magnetic field of 2500 G strength and 170 degrees of inclination to mimic the umbral conditions. We chose the latter value to match the observed Stokes V polarity and the small Stokes Q and U signals detected in the umbra. In addition, to simplify this comparison as much as possible we set the microturbulence constant with height and equal to one km s ^-1 (except in the FALS model where we kept the original microturbulence stratification). Regarding the macroturbulence, we set it to 2 km s ^-1, which is slightly above the 1.41 km s ^-1 of macroturbulence that corresponds to the 120 mÅ spectral resolution of GRIS at 1083 nm (Felipe et al. 2016). This parameter is kept fixed in the inversion since we expect negligible spatial variations of the macroturbulence in the umbra.

The synthetic Stokes I/I_c and V/I_c profiles are shown in Fig. 8. The Stokes I/I_c continuum levels of the Maltby and Avrett umbral models match the observed and deconvolved profile with 42% of scattered light. If we consider only 30% of scattered light, the continuum level is about 0.36 (see second panel in Fig. 3), which is far from the synthetized continuum umbral levels. Regarding the Stokes I/I_c wings, the synthesized profiles show a very similar behavior compared to the observed profile. However, the residual intensity at the line core in the synthesized profiles is much lower than the observed profile, i.e., the line core intensity in the observations is much higher. The two minima in the Stokes I/I_c synthesized profiles correspond to the well-known Zeeman splitting pattern due to a strong magnetic field. In particular these profiles correspond to the σ components, which are located at equidistant wavelength locations form the central wavelength. Their separation depends linearly on the field strength (for instance, see Landi Degl’Innocenti & Landolfi 2004). If the magnetic field is strictly pointing to the observer there is no central unshifted π component. Therefore, the intensity difference between the synthesized profiles and the observed profile at the line core wavelengths ascribed to any other physical mechanisms.

Regarding Stokes V/I_c, none of the models are able to generate a line reversal at the zero-crossing point. We could argue that the Stokes V line reversal appears because Stokes I is in emission at the line core. This would require the temperature to rise at lower atmospheric heights. However, this reasoning is too simplistic. The flanks of the line where the two σ components are located are deeper (in intensity units) than the line core and therefore the σ components would be “formed” higher in the atmosphere than the line core. Thus, if the temperature is rising at lower heights, the σ components in both Stokes I and V would be the first to appear in emission and not the line core. In general, the exact physical mechanisms are much more complex, particularly in NLTE conditions. As explained before, we discarded the magneto-optical effects as the source of the Stokes V zero-crossing reversal since the linear polarization signals are very small.

Fig. 8 Synthetic profiles from the semi-empirical models shown in Fig. 1. We also included in black the observed and deconvolved profile with 42% of scattered light shown in the second panel of Fig. 5 for comparison purposes. The observed profile is normalized to IQS, which is the average continuum intensity in quiet Sun areas. Ic stands for the continuum intensity at 500 nm.

4.2. NLTE inversion of the Si I 1082.7 nm line

We used NICOLE code to invert several selected pixels from the observed map shown in Fig. 2. We started the inversion using the FALS atmosphere as a guess model with a constant magnetic field of 2500 G, 170 degrees of inclination, and 90 degrees of azimuth. In order to invert each profile we performed three iterative cycles. In each one, we increased the number of positions in the atmosphere where the RF are computed (nodes) following the values shown in Table 1.

Fig. 9 Left panels: observed Stokes Iand V profiles (open circles connected with black solid line) and their corresponding fits (color) for two of the inversion runs. Right panels: temperature and magnetic field stratification, corresponding to the initial FALS model, and final inferred model by NICOLE, corresponding to cycles 2 and 3, are shown. The difference between the two cycles is in the number of nodes used for the vector magnetic field, which was a constant for one of these cycles (red) and three nodes for the other (blue).

This is the first time NLTE inversions of the Si i have been attempted and we must say that there is additional work that needs to be performed in the future. The most important part is reducing the complexity of the atomic model, if possible. The time required to perform a single inversion cycle ranges from 1~8 h using an intel i7-6700 CPU whose frequency ranges from 3.7–4.0 Ghz, i.e., it is included among the fastest frequencies that can be achieved with present CPUs without overclocking them. This long time per cycle means that for computing a full inversion, between 2–3 cycles, almost one day is required. This is more than 100 times slower than the NLTE inversion of the Ca ii 8542 Å line, which can be modeled with a much more simple atom (e.g., Shine & Linsky 1974). Thus, it would be desirable to work with simpler atom models, such as that proposed by Shchukina et al. (2017), to speed up the initial trial and error analysis before analyzing maps. In any case, NICOLE can work in parallel, i.e., we can speed up the inversion linearly with the number of cores available. In this section we discuss the results of the best-fitted profile. We would like to emphasize that it was extremely difficult to fit profiles such as those presented in the right-most column of Fig. 5, where the line core intensity is above the continuum level.

The inversion results corresponding to the second and third cycle are shown in Fig. 9. The NICOLE code fits to Stokes I are rather satisfactory overall. When we look into details, we can see that the line core is not very well fit when we use a model with constant magnetic field (about 3.5 kG and corresponding to cycle 2). The fit is more accurate when we increase the degrees of freedom in the magnetic field (cycle 3). In the latter case, the inferred magnetic field decreases linearly from 3.5 kG at log τ = 0 to 2.5 kG at log τ = −2 and then it increases again reaching 3.5 kG at log τ = −5. Regarding the temperature stratification, the inversion, in both cases, tries to introduce a hump around log τ ≈ −2 and an overall increase of the temperature beyond log τ = −2. Interestingly, none of the inversions are “apparently” able to fit the reversal in Stokes V at the zero-crossing point. As we mentioned in Sect. 3, the PCA deconvolution method is applied separately to Stokes I and V. So, there is a possibility of having two Stokes profiles that are not compatible at the central wavelength. For instance, the Stokes V profile, at the zero-crossing wavelength, is usually more affected by noise during the deconvolution process. Indeed, we have seen small variations around the Stokes V zero-crossing point when changing the number of iterations (see Appendix A). Another reason could be that the amount of stray light is too large; the reversal is recovered only when we take 42% of stray light, as opposed to the emission in Stokes I, which already appears with 30% levels. However, as we showed before, we do not recover any Stokes V reversal at the zero-crossing point either in the Ca I 1083.34 nm or Na I 1083.48 nm and Ca I 1083.9 nm lines, after applying the same deconvolution procedure and with noisier Stokes V profiles. Finally, the linear polarization profiles were also successfully fit in the inversion process. The inclination of the field is, at all heights, almost vertical (within 170 and 180 degrees).

5. Summary and discussion

In this paper we have shown that, at high spatial resolution and in the absence of scattered light, the Si I 1082.7 nm line is in emission in sunspot umbrae. To this end, we took observations in the near-infrared spectral range containing the Si I 1082.7 nm line with the GRIS instrument of the German GREGOR telescope at a spatial resolution of ≈0.̋3. Using a Mercury transit and the information provided by the adaptive optics system of GREGOR, we modeled the GREGOR PSF during the observations and corrected the data using a novel PCA deconvolution method. The results suggest there is a possibility that the Si I 1082.7 nm line appears in emission in sunspot umbrae with a line reversal in Stokes V at the zero-crossing point. We have shown how a large amount of stray light coming from wide angles can easily hide the emission peak in the Si I 1082.7 nm line. Only removing a 30% of the GREGOR stray light would suffice to start detecting the silicon line in emission. Remarkably, the applied amount of stray light has been determined by a Mercury transit, about 42%. With this amount of stray light the umbral continua of the recovered profiles match both: the sunspot umbral continuum levels predicted with the current model for a sunspot umbra of about 2.5 Mm radius (Kiess et al. 2014; Maltby et al. 1986) and the continuum level obtained by synthesizing the Si I 1082.7 nm spectral line with the umbral models of Avrett (1981) and Maltby et al. (1986). However, the latter models are not able to reproduce the line core intensity after taking into account NLTE effects. Remarkably, neither of the lines nearby the Si I 1082.7 nm line, i.e., the Ca I 1083.34 nm and Na I 1083.48 nm lines, which are both located at the sides of a Ti I 1083.36 nm atmospheric line, are in emission after the deconvolution.

We investigated whether the emission is a byproduct of the PCA deconvolution technique; but since we obtained the same results using an alternative deconvolution technique (based in the maximum entropy method), and because of the similarities of the observed GREGOR Stokes profiles with those coming from the German VTT, we think the recovered Stokes profiles are possibly real. Still, the presence of noise in the Stokes V profiles and the number of iterations can slightly alter the results around the zero-crossing point in Stokes V.

To analyze the Si I 1082.7 nm spectral line, and because it suffers from NLTE effects around the line core, we adapted the silicon atom model of Bard & Carlsson (2008) to NICOLE and used the latter code to invert the observed Stokes profiles. The NICOLE inversion shows that to get the emission in the Si I 1082.7 nm line, the temperature stratification should have a hump located at about log τ = −2 and then should start rising at lower heights when going into the transition region. This result may explain the controversy between the present empirical models and the radio observations at the millimetre range (Iwai et al. 2016; Iwai & Shimojo 2015). The fit is better when we increase the degrees of freedom in the magnetic field strength although the recovered field strength stratification, which first decreases linearly from 3.5 kG at log τ = 0 to 2.5 kG at log τ = −2 and then increases again afterward, seems unrealistic; in magnetostatic equilibrium and in a stratified atmosphere (with gas pressure and density decreasing with height) the magnetic field should decrease with height. Remarkably, NICOLE is unable to reproduce the Stokes V line reversal. The reason may be that the reversal is not real either because the noise levels have increased after the deconvolution (see Appendix A) or the amount of stray light is too large; the reversal is recovered only when we take 42% of stray light, as opposed to the emission in Stokes I, which already appears with 30% levels. We do not detect the Stokes V line reversal when deconvolving the Ca I 1083.34 nm and Na I 1083.48 nm lines, and in this case, the noise in polarization is larger since the signals are much smaller than for Si I 1082.7 nm.

In summary, the current results point to a temperature stratification in sunspot umbrae that is much more complex that that provided by current umbral models. The main relevance of these results is that the temperature minimum region may be located at a lower height above sunspot umbrae than in the photosphere. This particular issue has been already discussed by Severino et al. (1994), where the authors point out that the temperature gradient at the high photosphere should be steeper than in current models.

Our finding provides insights to understand why the sunspot umbra emission in the millimeter spectral range is less than that predicted by current empirical umbral models. However, we need to carry out a more systematic and conscientious investigation of the formation and sensitivities of the Si I 1082.7 nm line, using numerical simulations and additional observations, to fully understand the formation of this line in strong magnetic field concentrations, such as network, pores, and sunspots. Other effects, such as the influence of the Zeeman splitting in the determination of the level populations or abundance variations, has not been taken into account in the inversion process. A deep analysis of the sensitivity of the Si I line, including the magnetic field in the statistical equilibrium equations, is a work in progress. The next step will be to analyze the full sunspot, including the Ca I 1083.34 nm and Na I 1083.48 nm lines. Thus, the results presented here represent the first step toward the inversion of the Si I 1082.7 nm line including NLTE effects.

Acknowledgments

Special thanks to Mats Carlsson (who provided the atom model) and Han Uitenbroek for their support during the NLTE analysis of this work. The 1.5 m GREGOR solar telescope was built by a German consortium under the leadership of the Kiepenheuer-Institut für Sonnenphysik in Freiburg with the Leibniz-Institut für Astrophysik Potsdam, the Institut für Astrophysik Göttingen, and the Max-Planck-Institut für Sonnensystemforschung in Göttingen as partners, and with contributions by the Instituto de Astrofísica de Canarias and the Astronomical Institute of the Academy of Sciences of the Czech Republic. This work was partly supported by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea. This study is supported by the European Commissions FP7 Capacities Programme under the Grant Agreement number 312495. The GRIS instrument was developed thanks to the support by the Spanish Ministry of Economy and Competitiveness through the project AYA2010-18029 (Solar Magnetism and Astrophysical Spectropolarimetry) and AYA2014-60476-P. This work has also been supported by Spanish Ministry of Economy and Competitiveness through projects ESP2014-56169-C6-1-R and ESP-2016-77548-C5-1-R.

References

Appendix A: Variations of deconvolution parameters

The PCA deconvolution method has two fundamental tuneable parameters: (i) the number of profiles used to construct the database and (ii) the number of iterations used within the Richardson-Lucy deconvolution internal algorithm. The determination of the database size is not trivial. It is usually carried out manually by looking at the profiles. Usually, the user cuts the database at a level where the profiles look like pure Poisson noise. In our case, the level is about n = 12. If we decrease the number of profiles in the database (n < 12) the PCA deconvolution algorithm can miss-reproduce the observed profiles. The number of iterations used for the deconvolution can also alter the final result. Too many iterations can amplify the noise in the final deconvolved profiles. Figure A.1 shows how the Stokes I/I_QS (left panels) and Stokes V/I_QS (right panels) change for different n (top panels) and number of iterations (bottom panels). It can be seen that Stokes I/I_QS remains about the same if n> 8. For small n values (small database) the recovered Stokes I/I_QS profile is completely different. In case of Stokes V/I_QS, it is clear that the use of a small database can produce the line reversal at the zero-crossing wavelength point. Regarding the number of iterations, there is no apparent change in Stokes I/I_QS when the number of iterations varies from 5 to 15, except a slight change in the continuum. In Stokes V/I_QS there are tiny variations around the zero-crossing points. We also show, in the top panels, the result using a smaller field of view when deconvolving the maps. This is important because the map contains a spot and although the amount of scattered light is constant for every pixel, it is more critical when the image contrast is large (as in a sunspot).

Fig. A.1 Variations of a reconstructed umbra Stokes I/IQS (left panels) and Stokes V/IQS (right panels) profiles when changing the number n of the profiles within the PCA database and the number of iterations.

Appendix B: Deconvolution of nearby lines

Figure B.1 shows the original (blue color) and deconvolved (black color) Stokes I/I_QS and V/I_QS profiles corresponding to the Ca I 1083.34 nm and Na I 1083.48 nm lines, both located at the sides of a Ti I 1083.36 nm atmospheric line (top panels) and the Ca I  line at 1083.80 nm located further in the red (bottom panels). The very first important point to highlight is that none of the Stokes V/I_QS get a reversal at the zero-crossing wavelengths. In the Ca I  line at 1083.80 nm we do see some variations in the profile at the umbral core (right most bottom panel) but these are due to noise. In the case of Stokes I/I_QS it is interesting to see how the core of the Na I 1083.48 nm line does not reach local continuum values. This is not the case for the Ca I  lines, especially for the Ca I  line at 1083.80 nm (bottom panels). In this case, the splitting is clearly visible in the line and the residual intensity at the line core of the darker profiles (right most panel) is larger than the continuum. The same happens in the Ca I 1083.34 nm line, although it is less noticeable. The fact that the Na i line does not behave similarly suggests that the higher continuum intensity values of the line core in the Ca I  lines may be due to the noise. Both Ca I  lines are very shallow. However, in practice we cannot distinguish this effect from an excess of scattered light in the deconvolution.

Fig. B.1 Variation of the original (blue color) and reconstructed (black color) Stokes I/IQS and V/IQS profiles corresponding to Ca I 1083.34 nm, Na I 1083.48 nm, and Ca I  1083.80 nm as we move from the inner penumbral boundary (first column) to the center of the umbra (last column). The pixels correspond to those shown in Fig. 3. The dotted profile in the first panel corresponds to a quiet Sun profile (multiplied by ×0.5 to fit within the Y-axis range). The horizontal dotted lines are drawn for reference.

All Tables

All Figures

Fig. 1 Variation of the temperature with atmospheric height corresponding to various atmospheric models available in the literature. In particular, the E model of Avrett (1981), the models by Maltby et al. (1986) and Fontenla et al. (2009) that cover from the photosphere to the chromosphere, and the model of Collados et al. (1994; cool umbra and hot umbra) that only covers the photosphere. The optical depth scale has been calculated using the Maltby et al. (1986) atmospheric model.

In the text

Fig. 2 Intensity maps of the observed NOAA 12096 sunspot corresponding to the continuum (left panels) and to the Si I 1082.7 nm line core (center panels). The top panels represent the original data and the bottom panels the reconstructed data. The right-most panels show a zoom over the original (top) and reconstructed (bottom) sunspot umbra in the continuum (see the square in the top left panel). The contours outline the region where the Si I 1082.7 nm continuum is below 25% of the average quiet Sun.

In the text

Fig. 3 Variation of the original (blue) and reconstructed Stokes I/IQS and V/IQS profiles as we move from the inner penumbral boundary (first column) to the center of the umbra (last column). The Stokes I profiles in black correspond to the deconvolved profiles taking into account 42% of stray light while for the red profiles we only accounted for 30% of the stray light. For simplicity, we only represent five pixels. The location of the pixels can be found in the right panels of Fig. 2. The first profile corresponds to the lowest point.

In the text

	Fig. 4 Variation of umbra continuum intensity with wavelength (right panel) for three different models (left): the FALS, the cool model of Collados et al. (1994), and a super-cool model obtained from the cool model. The vertical line shows 1083 nm wavelength while horizontal lines pinpoint the corresponding continuum intensity values at 1083 nm for the three models.
In the text

	Fig. 5 Variation of Stokes I/IQS and V/IQS profiles with the amount of scattered light β from 0% to 42%. The profile corresponds to that shown in the fifth panel on Fig. 3.
In the text

Fig. 6 Left panels: umbra Stokes I/IQS and V/IQS profiles reconstructed using the PCA deconvolution and maximum entropy methods. Right panels: umbra Stokes I/IQS and V/IQS profile as seen by GREGOR and the German VTT. The solid line represents the non-reconstructed data while the dotted line stands for the reconstructed data. The black line corresponds to the German VTT (simulated), red to the GREGOR telescope, and blue to the observed VTT profile (see text).

In the text

	Fig. 7 Spectral atlas of the Si i line (solid). We include the results of a LTE (dashed) and NLTE (dot-dashed) synthesis using the FALC model atmosphere.
In the text

	Fig. 8 Synthetic profiles from the semi-empirical models shown in Fig. 1. We also included in black the observed and deconvolved profile with 42% of scattered light shown in the second panel of Fig. 5 for comparison purposes. The observed profile is normalized to IQS, which is the average continuum intensity in quiet Sun areas. Ic stands for the continuum intensity at 500 nm.
In the text

Fig. 9 Left panels: observed Stokes Iand V profiles (open circles connected with black solid line) and their corresponding fits (color) for two of the inversion runs. Right panels: temperature and magnetic field stratification, corresponding to the initial FALS model, and final inferred model by NICOLE, corresponding to cycles 2 and 3, are shown. The difference between the two cycles is in the number of nodes used for the vector magnetic field, which was a constant for one of these cycles (red) and three nodes for the other (blue).

In the text

	Fig. A.1 Variations of a reconstructed umbra Stokes I/IQS (left panels) and Stokes V/IQS (right panels) profiles when changing the number n of the profiles within the PCA database and the number of iterations.
In the text

Fig. B.1 Variation of the original (blue color) and reconstructed (black color) Stokes I/IQS and V/IQS profiles corresponding to Ca I 1083.34 nm, Na I 1083.48 nm, and Ca I  1083.80 nm as we move from the inner penumbral boundary (first column) to the center of the umbra (last column). The pixels correspond to those shown in Fig. 3. The dotted profile in the first panel corresponds to a quiet Sun profile (multiplied by ×0.5 to fit within the Y-axis range). The horizontal dotted lines are drawn for reference.

In the text