Magnetic properties of a long-lived sunspot - Vertical magnetic field at the umbral boundary

Context. In a recent statistical study of sunspots in 79 active regions, the vertical magnetic field component $B_\text{ver}$ averaged along the umbral boundary is found to be independent of sunspot size. The authors of that study conclude that the absolute value of $B_\text{ver}$ at the umbral boundary is the same for all spots. Aims. We investigate the temporal evolution of $B_\text{ver}$ averaged along the umbral boundary of one long-lived sunspot during its stable phase. Methods. We analysed data from the HMI instrument on-board SDO. Contours of continuum intensity at $I_\text{c}=0.5I_\text{qs}$, whereby $I_\text{qs}$ refers to the average over the quiet sun areas, are used to extract the magnetic field along the umbral boundary. Projection effects due to different formation heights of the Fe I 617.3 nm line and continuum are taken into account. To avoid limb artefacts, the spot is only analysed for heliocentric angles smaller than $60^{\circ}$. Results. During the first disc passage, NOAA AR 11591, $B_\text{ver}$ remains constant at 1693 G with a root-mean-square deviation of 15 G, whereas the magnetic field strength varies substantially (mean 2171 G, rms of 48 G) and shows a long term variation. Compensating for formation height has little influence on the mean value along each contour, but reduces the variations along the contour when away from disc centre, yielding a better match between the contours of $B_\text{ver}=1693$ G and $I_\text{c}=0.5I_\text{qs}$. Conclusions. During the disc passage of a stable sunspot, its umbral boundary can equivalently be defined by using the continuum intensity $I_\text{c}$ or the vertical magnetic field component $B_\text{ver}$. Contours of fixed magnetic field strength fail to outline the umbral boundary.


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 |B ver | > B stable ver , whereas outside, the penumbral mode takes over.
videos associated with Fig. 3 are available at http://www.aanda.org 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 |B ver | < B stable ver , the penumbral mode of magneto-convection takes over the umbral mode. Jurčák et al. (2018) extends the analysis of 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 log area 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 |B ver | > B stable ver , 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 |B ver | < B stable ver 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

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 http://jsoc.stanford.edu/ajax/ exportdata.html 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.Ic_noLimbDark_720s & hmi.B_720s. 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 = 696Mm 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 & 0.5 I 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, and2001b), 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 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 1 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 X(t) ∈ B ver ψ (t) , |B| ψ (t), γ lrf ψ (t) . Similarly, standard deviations along the contours σ ψ (t) were calculated. These time series (c.f 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 whereby t is in days and t ∈ N is at noon. X 0 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 & t end have a different time of day. Here we have | X(t) t − X 0 | < 0.5 G for all X(t) in G and < 0.01 • for X(t) = γ lrf ψ (t). While X 1 & X 2 are used internally during the fitting process to guarantee numeric stability, the results are presented with parameters X 0 , X 3 , and X 4 in Table 2. X 3 and X 4 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 0 from the fits to B ver ψ (t), |B| ψ (t), and γ lrf ψ (t) 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 1 We note the subtle difference between the heliocentric angle and the angle between the LOS and the local vertical. The heliocentric angle, θ, is the angle between the centre of the umbra and the observer as measured from the centre of the sun. The angle, α, between the LOS and the local vertical at the umbral centre is given by: α = θ + r, whereby r = x 2 + y 2 . For any position on the solar disc, r is smaller than r ≈ 0.27 • . The angle, α, is used to transform between the LOS and LRF coordinate systems.
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.5I 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, Y(ψ) ∈ {B ver (ψ), |B|(ψ), γ lrf (ψ)} is least square fitted against functions of the form 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 Y 3 (t) t 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 0 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.

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.5 I qs and of Fig. 2 for I c = 0.4 I qs . The lower panels show the residuals after subtracting the fit.

Temporal evolution
The parameters of the sinusoidal fits, offset X 0 , amplitude X 3 , 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.5 I qs , we find for B ver ψ (t) that σ t = 15 G is smaller than the orbital amplitude X 3 = 18 G, with an offset of X 0 = 1693 G. For the contours of I c = 0.4 I qs , σ t = 19 G is also smaller than X 3 = 20 G with an X 0 = 1850 G. For the residuals of B ver no long-term trend is noticeable.
In contrast, the residuals of |B| ψ (t) 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 0 for |B| ψ (t) and γ lrf ψ (t) are 2171 G and 141.4 • respectively at the contours with I c = 0.5 I qs . The fact that the residuals σ t ( B ver ψ (t)) are smaller than σ t ( |B| ψ (t)) 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 B ver ψ (t) offset, X 0 , between the two different intensities amounts to 157 G while that of |B| ψ (t) 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 B ver ψ (t) can be considered constant in time.

Contours
Using the offsets X 0 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 100x100 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.
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 at http://www.aanda.org. 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, d ψ,t (cf. Eq. 3). It is given in the last column of Table 2. d ψ,t 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, d ψ,t = 0.44, is found for I = 0.53I qs with −B ver = 1639 G (X 3 = 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.

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 0 , X 3 , and σ t change only marginally. A plot like in Fig. 1 with ∆h = 0 km looks almost identical (not shown).
However, σ ψ t and d ψ,t increase significantly. For example, for B ver at I = 0.5I qs , σ ψ t and d ψ,t 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 d ψ,t , 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 at http://www.aanda.org

Conclusion
Investigating the physical properties along the umbra-penumbral boundary of a stable sunspot for a time span of approximately ten, 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.5 I 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 |B ver | = 1867 +18 −16 G (99%-error) at I c = 0.5 I qs . The values for |B ver | differ by some 175 G. In general, a difference is expected due to differences in the experimental setup