© ESO 2018

1. Introduction

The boundary between umbra and penumbra of sunspots has long been defined in terms of the continuum intensity I_c. This brightness difference is the consequence of the different magneto-convective processes running in umbrae and penumbrae. We have evaluated magnetic quantities to identify which of them may cause the different behaviour on the two sides of the umbral boundary.

Jurčák (2011) investigated the properties of the magnetic field at umbral boundaries and noted that the vertical magnetic field component |B_ver| changes little along the boundaries of the ten sunspots he analysed and could neither verify nor falsify a dependence of the median value along the boundary on the area of the sunspot. The ten-spot average of the median along the boundary was 1860 G, whereas the mean of the standard deviations along the boundary was given as 190 G for Hinode/SP data.

Jurčák et al. (2015) extended the analysis by investigating a 4.5h time series of a forming sunspot using GFPI/VTT data and noting an increase of |B_ver| at the migrating umbral boundary during penumbra formation and stabilization of this value after completion of the formation. Shortly thereafter, that part of the umbral boundary was observed with Hinode/SP and a |B_ver| value of 1810 G measured. They propose that the umbral mode of magneto-convection prevails in areas with , whereas outside, the penumbral mode takes over.

Following this line of investigation, Jurčák et al. (2017) studied a pore whose |B_ver| remained below this critical value. They found that a developing penumbra completely cannibalized the pore, thus supporting the assertion that in umbral areas with , the penumbral mode of magneto-convection takes over the umbral mode.

Jurčák et al. (2018) extends the analysis of Jurčák (2011) to 88 scans of 79 different active regions again using Hinode/SP and showed that the I_c = 0.5I_qs contours match mostly the |B_ver|=1867 G contours. A Bayesian linear regression showed that a model with constant |B_ver| is more likely to explain the data than a first or second order polynomial with  logarea   as independent variable. Furthermore the most likely |B_ver|=1867 G, with a 99% probability for 1849 G ≤ |B_ver|≤1885 G. A dependence on the solar cycle could not be verified.

These findings have led to the Jurčák criterion, an empirical law stating that the umbral boundary of stable sunspots can be equivalently defined by a continuum intensity I_c or a vertical magnetic field component |B_ver|. In other words, in areas with , only the umbral mode of convection exists, hindering other modes of magneto-convection. A conjecture can also be stated from these findings: umbral areas with are unstable against more vigorous modes of convection, that is, they are prone to vanish.

In this work we have investigated the behaviour of the magnetic field along the umbral boundary in a time series of a single stable sunspot. We used the spot of NOAA AR 11591 during its first disc passage. This allows us to verify whether remains constant over ≈10 days, which would provide support to the Jurčák criterion. Hereby ⟨ ⋅ ⟩_ψ stands for average along the I_c contour.

2. Data and analysis

The used data are retrieved after processing by the Solar Dynamics Observatory’s (SDO) Helioseismic and Magnetic Imager (HMI) Vector Magnetic Field Pipeline (Hoeksema et al. 2014) cutout service for NOAA AR 11591. Using this NOAA AR number on the JSOC website¹ in the im_patch processing option automatically gives the reference coordinates listed in the final three columns of Table 1. t_{< E} and t_{> W} are the first and last time steps processed, where t_{< E} is before the sunspot rotates over the east limb onto the sun and t_{> W} is after the sunspot rotates off the west limb. A cutout size of 500 × 500 pixel was chosen. The data series used are hmi.IcnoLimbDark720s and hmi.B720s. For the full disc passage, there are 1599 time steps.

For the 180°-disambiguation the potential acute solution provided by the pipeline was adopted. This can be done using hmi_disambig with method=0. We note that for all pixels 180° must be added because the azimuth is defined relative to the positive y-axis of the maps in CCD-frame and exportdata’s im_patch option rotates the maps 180° so that solar north is up.

The heliographic Stonyhurst coordinates are calculated using procedures modified from and tested against sswidl’s wcs routines fitshead2wcs, wcs_get_coord, wcs_convert_from_coord and those they call (see Thompson 2006). The canonical value for HMI of R_⊙ = 696 Mm is used. The transformation of the magnetic field vector into the local reference frame was performed with a code modified from and tested against Xudong Sun’s sswidl routine hmi_b2ptr (see Gary & Hagyard 1990; Thompson 2006; Sun 2013).

Quiet sun intensity. The limb darkening correction in the HMI pipeline was based on Pierce & Slaughter (1977, Eq. (9)), which does not consider all orbital artefacts introduced into the continuum intensity I_c of SDO/HMI data. Even after limb darkening removal and normalization there is a change over the day in I_c of the order of 1% towards the limb with opposite signs on the western and eastern hemisphere. To compensate for this, the quiet sun intensity I_qs for each time step was chosen such that I_qs is the mean of all the quiet sun pixels within the 500 × 500 cutout, where quiet sun is defined as having I_c >  0.9I_qs.

Contours were taken at I_c = 0.4 and 0.5I_qs, and the positions of the contours are used to interpolate the values of the vertical magnetic field component B_ver, the magnetic field strength |B| and the inclination to the surface normal γ_LRF. Vertical is to be understood in the local reference frame, in other words, it is the direction of the surface normal.

Due to different formation heights of I_c and the Fe I 617.3 nm line, as well as the Wilson depression (Wilson 1774) and differential line-of-sight opacity effects (see e.g. Rimmele 1995; Westendorp Plaza et al. 2001a,b), the magnetic contours are projected towards the limb (i.e. outwards) relative to the intensity contours. To compensate for these shifts and get a better match between I_c and B_ver contours, we transformed the coordinates obtained from I_c contours, (x, y), using

(1)

before retrieving the magnetic field values at coordinates (x′,y′). (x, y),(x′,y′) are helio-projective coordinates in arc-seconds from disc centre and Δh is the formation height difference. Later on, when the contours from magnetic field maps are plotted onto the I_c map (cf. Sect. 3.2 and Figs. 3 and 5), the inverse of Eq. (1) is applied, meaning that the magnetic contours are shifted inward. The value of Δh = 465 km results from a minimization procedure, which is explained on page 4. The effect of neglecting this compensation is discussed in Sect. 3.3.

The limits of the time series we analyse are given as t_start and t_end in Table 1. A total of 1063 time steps in this time range are available. This time range was chosen to select data sets, for which the heliocentric angle² of the centroid of the umbra was smaller than 60°.

Time series fit. For every time step and magnetic quantity, an average was computed along the contours, thereby creating time series of the form , , . Similarly, standard deviations along the contours σ_ψ(t) were calculated. These time series (cf. Sect. 3 and Figs. 1 and 2) show a daily variation of an approximately sinusoidal shape. We believe them to be an artefact of SDO’s geosynchronous orbital motion. For the ranges from t_start to t_end given in Table 1, these time series are least square fitted against functions of the form

Fig. 1. Mean magnetic field strength in black (it’s vertical component in blue) along the Ic = 0.5Iqs contour, with Δh = 465 km accounted for, Sinusodial fits and the residuals for NOAA AR 11591.

(2)

whereby t is in days and t ∈ ℕ is at noon. X₀ is the value we are interested in and will be henceforth called offset. It is used instead of a time average ⟨X(t)⟩_t because it correctly accounts for missing data (most importantly the gap in the afternoon of Oct 17) and that t_start and t_end have a different time of day. Here we have |⟨X(t)⟩_t − X₀| < 0.5 G for all X(t) in G and < 0.01° for . While X₁ and X₂ are used internally during the fitting process to guarantee numeric stability, the results are presented with parameters X₀,  X₃,  and X₄ in Table 2. X₃ and X₄ are the amplitude and phase of the orbital artefacts. Also listed are the standard deviations of the residuals σ_t = σ(X(t)−X_fit(t)) and the means of the standard deviations along the contours over the same range in time ⟨σ_ψ⟩_t.

Levels of magnetic contours. The offsets X₀ from the fits to , and for the 0.5(0.4)I_qs contours are then used as contour level on the B_ver, |B| and γ_LRF maps, respectively. They are discussed in Sect. 3.2 and plotted in Figs. 3 and 5 and the videos.

Distance between contours. To quantify how well two contours match we calculated the average distance between them ⟨d⟩_ψ, which we define as the area of symmetric difference divided by the length of the intensity contour, ℓ(t). The area of symmetric difference, Δa(t), is the area surrounded by either of the contours but not both. When averaging in time we weighed by the contour length, giving

(3)

These average distances between contours are listed in Table 2 in pixel. For ⟨d⟩_ψ ≪ 1pixel only the total ordering should be relied upon due to griding and other computational effects.

Fit along each contour. For every point along a contour, a reference angle ψ = ∡(PCD) is calculated, whereby P is the point on the contour, C is the centroid of the I_c = 0.5 I_qs contour in the CCD frame and D is the centre of the solar disc as observed by SDO. The angles are calculated on the sphere. For every time step and every contour, is least square fitted against functions of the form

(4)

Those fits are plotted in the right panels of the videos (cf. Sect. 3.2 and bottom panels of Figs. 3 and 5). Furthermore the time averages of the fit amplitudes are listed in Table 2.

Optimal height difference.Δh = 465 km was chosen because it minimizes the average distance ⟨d⟩_ψ, t between the I_c = 0.5I_qs contours after transformation with Eq. (1) and the B_ver contours, whereby the contour level X₀ on the B_ver map has been derived with the fit to ⟨B_ver⟩_ψ(t) as described above (Eq. (2)). An optimal height difference of Δh = 465 km means that the intensity contour at the limb is shifted outwards by 465 km ⋅ r_⊙/R_⊙ ≈ 0.65″ ≈ 1.3 pixel. The difference of the formation heights for continuum and Fe I 617.3 nm line core amounts to ≈250 km for a typical umbral model atmosphere (see e.g. Norton et al. 2006, Table 1). The fact that the value for Δh is larger may be explained with the Wilson depression of the umbra, which typically amounts to 800 km. The latter causes the τ = 1 surface to be strongly inclined relative to horizontal. Minimizing the standard deviation of B_ver along the I_c = 0.5I_qs contour (⟨σ_ψ⟩_t column in Table 2) instead would give an optimal Δh = 520 km.

3. Results

Based on the time series of approximately ten days, in which the spot of NOAA AR 11591 has heliocentric angles smaller than 60°, we determine the magnetic properties for two distinct contour levels of the continuum intensity. As intensity levels we use I_c = 0.5 (0.4)I_qs. Along each contour, the azimuthal average of B_ver, |B| and γ_LRF are calculated. The respective values of those averages for B_ver (in blue) and |B| (in black) as well as sinusoidal fit of the orbital variation are displayed in the upper panels of Fig. 1 for I_c = 0.5I_qs and of Fig. 2 for I_c = 0.4I_qs. The lower panels show the residuals after subtracting the fit.

Fig. 2. As Fig. 1, but from contours at Ic = 0.4Iqs.

3.1. Temporal evolution

The parameters of the sinusoidal fits, offset X₀, amplitude X₃, and the rms of the corresponding residuals, σ_t, are given in Table 2 for all considered cases. In addition, they are printed into the plots of Figs. 1 and 2. For the contours at I_c = 0.5I_qs, we find for that σ_t = 15 G is smaller than the orbital amplitude X₃ = 18 G, with an offset of X₀ = 1693 G. For the contours of I_c = 0.4I_qs, σ_t = 19 G is also smaller than X₃ = 20 G with an X₀ = 1850 G. For the residuals of B_ver no long-term trend is noticeable.

In contrast, the residuals of amount to σ_t = 48 G which is larger than the amplitudes of the sinusoidal fit (16 G), and it shows a long-term variation. Since γ_LRF is dependent on B_ver and |B|, it has a long-term variation which compensates for that of |B| (not shown). The offsets X₀ for and are 2171 G and 141.4° respectively at the contours with I_c = 0.5I_qs.

The fact that the residuals are smaller than is remarkable, but is even more remarkable if one considers that the gradient of B_ver perpendicular to the contour is larger than that of |B|. This can be inferred from Table 2: the difference of the offset, X₀, between the two different intensities amounts to 157 G while that of is only 94 G. Hence, a small shift of the contour implies a larger deviation in B_ver than in |B|. Therefore, our result of a smaller deviation in B_ver relative to |B| gives further evidence that can be considered constant in time.

3.2. Contours

Using the offsets X₀ from the fits in Table 2 with I_c = 0.5I_qs and Δh = 465 km, the upper panels of Fig. 3 overplot the contours of intensity I_c = 0.5I_qs (red), |B|=2171 G (green), −B_ver = 1693 G (blue), and γ_LRF = 141.4° (yellow). The background images consist of 100 × 100 pixel cutouts of grey-scale intensity maps with a minimum (maximum) of I_c = 0.1 (1.2)I_qs. A close inspection of the figure shows that the B_ver contour matches best with the intensity contour. The cyan arrow originates in the centroid of the umbra and points towards disc centre. The centroid is determined by the I_c = 0.5I_qs contour and is derived using CCD coordinates.

Fig. 3. NOAA AR 11591, top panels: continuum intensity maps for longitudes −60°, 0° and 30°. The legend in the lower right corner of the top left panel defines the contour levels. Different formation heights are accounted for (Eq. (1), Δh = 465 km). The cyan arrow originates in the centroid of the umbra and points towards disc centre. Bottom three rows: magnetic field parameters retrieved along the Ic = 0.5Iqs contour. The temporal evolution is available online.

The three bottom rows of panels of Fig. 3 show the magnetic field quantities along the I_c = 0.5I_qs contour as well as their sinusoidal fits in black. The azimuth is determined relative to the centroid and the direction towards disc centre, which corresponds to ψ = 0° and runs counter-clockwise.

To quantify the azimuthal variation of the magnetic parameters, Table 2 gives the time average of the standard deviations along the contours, ⟨σ_ψ⟩_t. Again ⟨σ_ψ⟩_t is smaller for B_ver (81 G) than for |B| (111 G). As before, the small value for B_ver is remarkable, since its gradient perpendicular to the contour is larger than for |B|. The lower panels demonstrate that the azimuthal variations are smallest for B_ver. Again, we note that this is remarkable considering the fact that the gradient of B_ver perpendicular to the contour is larger than the gradient of |B|.

A video of the temporal evolution of those contours during the disc passage of the spot is available online. This animation demonstrates that an iso-contour of B_ver = −1693 G coincides nicely with the intensity contour at 0.5I_qs. This animation also demonstrates that contours of |B| and γ_LRF do not coincide.

To quantify the match or mismatch of two contours, we have introduced the average distance between two sets of contours, (cf. Eq. (3)). It is given in the last column of Table 2. is smallest for the B_ver contours with Δh = 465 km (see the first and final row of Table 2).

In Fig. 4 the average distance is plotted for intensities changing from 0.30 to 0.65. The corresponding contour levels for B_ver are calculated as described in Sect. 2 (fit to Eq. (2)). The best match, , is found for I = 0.53I_qs with −B_ver = 1639 G (X₃ = 17 G, σ_t = 15 G, and ⟨σ_ψ⟩_t = 82 G). Distances for |B| and γ_LRF are in all cases larger and not plotted. Hence, by minimizing the distance, −B_ver = 1639 G results as the value that defines the umbral boundary at I = 0.53I_qs. This is additional proof that our chosen value of I = 0.5I_qs is very close to the optimum value.

Fig. 4. Average distance, , between contours of varying intensity and Bver. The contour level of Bver is determined by fitting Eq. (2). has a minimum for Ic = 0.53Iqs corresponding to −Bver = 1639 G.

3.3. Effect of neglecting formation heights compensations

For the results presented so far, we corrected for the projection effects due to different formation heights of continuum and line. As discussed in the end Sect. 2 we assume a height difference of Δh = 465 km. Table 2 also gives the results for the case in which these projection effects are not considered, i.e. Δh = 0 km. As a general trend, it is seen that the values for X₀, X₃, and σ_t change only marginally. A plot like in Fig. 1 with Δh = 0 km looks almost identical (not shown).

However, ⟨σ_ψ⟩_t and increase significantly. For example, for B_ver at I = 0.5I_qs, ⟨σ_ψ⟩_t and increase by more than 30% from 81 to 113 G, and from 0.45 to 0.59 pixel, respectively. This is illustrated in Fig. 5, which shows the same snapshot as in the left column of Fig. 3, with the only difference that Δh = 0 km. In this case, the heliocentric angle is 60°. It is seen that the magnetic contours are shifted relative to the intensity, which results in an increase of , and the variation of B_ver along the contour (bottom panels) are larger for Δh = 0 km. This can also be seen in the corresponding video of the disc passage of the spot, which is available online.

Fig. 5. As Fig. 3, left column, but without compensation for different formation heights (Δh = 0). The temporal evolution is available online.

4. Conclusion

Investigating the physical properties along the umbra-penumbral boundary of a stable sunspot for a time span of approximately ten days, we find three main results:

1. B_ver averaged along the I = 0.5I_qs contour is nearly constant in time.

2. Contours of intensity and of B_ver match at the umbral boundary. The best match is obtained for I = 0.53I_qs and |B_ver|=1639 G.

3. Projection effects due to different formations height of the spectral line and continuum need to be considered. If not, variation of B_ver along the contour increases significantly.

These results are obtained by analysing 1063 consecutive SDO/HMI data sets (with a time step of 12 min) of the first disc passage of NOAA AR 11591.

Using I_c = 0.5I_qs to define the umbral boundary, we obtain |B_ver|  =  1693 G  ±  15 (1σ_t-error). Jurčák et al. (2018) used Hinode/SP data to find (99%-error) at I_c  =  0.5I_qs.

The values for |B_ver| differ by some 175 G. In general, a difference is expected due to differences in the experimental setup and analysis methods. Sainz Dalda (2017) investigates the differences between HMI and SP vector magnetograms and obtained comparable differences. He concludes that the filling factor followed by spatial and spectral resolution are the main source. At the umbral boundary the filling factor is 1, and causes therefore no differences. The other effects are particularly strong at the sharp boundary between umbra and penumbra, where the intensity gradient is large.

Hence, these investigations provide evidence that |B_ver| is constant for a statistical sample of sunspots as well as during the evolution of one stable spot, thereby supporting the Jurčák criterion.

Acknowledgments

We wish to thank Jan Jurčák, Juan Manuel Borrero and the anonymous reviewer for valuable discussions, Xudong Sun for making available the vector transformation routine as well as Hanna Strecker and various members of the IDL user group (http://www.idlcoyote.com/comp.lang.idl-pvwave/) for help with IDL and tex.sx (https://tex.stackexchange.com/) users for help with LaTeX. The data used is courtesy of NASA/SDO and the HMI science team (see e.g. Metcalf 1994; Leka et al. 2009; Borrero et al. 2011; Pesnell et al. 2012; Schou et al. 2012; Hoeksema et al. 2014). This research has made use of NASA’s Astrophysics Data System (http://adsabs.harvard.edu/abstract_service.html).

References

Movies

Movie of Fig. 3 (Original mp4) (Access here)

Movie of Fig. 3 (Original mov) (Access here)

Movie of Fig. 5 (Original mp4) (Access here)

Movie of Fig. 5 (Original mov) (Access here)

All Tables

All Figures

	Fig. 1. Mean magnetic field strength in black (it’s vertical component in blue) along the Ic = 0.5Iqs contour, with Δh = 465 km accounted for, Sinusodial fits and the residuals for NOAA AR 11591.
In the text

	Fig. 2. As Fig. 1, but from contours at Ic = 0.4Iqs.
In the text

Fig. 3. NOAA AR 11591, top panels: continuum intensity maps for longitudes −60°, 0° and 30°. The legend in the lower right corner of the top left panel defines the contour levels. Different formation heights are accounted for (Eq. (1), Δh = 465 km). The cyan arrow originates in the centroid of the umbra and points towards disc centre. Bottom three rows: magnetic field parameters retrieved along the Ic = 0.5Iqs contour. The temporal evolution is available online.

In the text

	Fig. 4. Average distance, , between contours of varying intensity and Bver. The contour level of Bver is determined by fitting Eq. (2). has a minimum for Ic = 0.53Iqs corresponding to −Bver = 1639 G.
In the text

	Fig. 5. As Fig. 3, left column, but without compensation for different formation heights (Δh = 0). The temporal evolution is available online.
In the text