© C. Baumgartner et al. 2022

1. Introduction

The magnetic field in solar active regions is often modeled by coherent bundles of magnetic field lines, so-called flux tubes. The magnetic helicity H = ∫_VA ⋅ B dV, where A is the magnetic vector potential and B is the magnetic field, can be used to describe the topological structure of flux tubes fully contained in a volume V (Berger & Field 1984). The magnetic helicity of a single flux tube has two components: writhe, which measures the deformation of its axis, and twist. If one imagines the magnetic field as a straight ribbon with its two ends rotated in opposite directions, the twist T measures how often the ribbon turns around its straight axis

$$\begin{array}{c}\hfill T=qL,\end{array}$$(1)

where L is the length of the ribbon and q is the twist density, which counts how often the ribbon fully turns per unit length.

Measurements of the twist of magnetic field play an important role in many different areas of solar physics: The twist distribution in the photosphere constrains models of the solar dynamo and magnetic flux emergence (e.g., Gilman & Charbonneau 1999; Brandenburg 2005; Pipin et al. 2013). For example, the hemispheric helicity sign rule describes an observed latitudinal dependence of the twist with predominantly negative (counter-clockwise) or positive (clockwise) twist in the northern or southern hemisphere (Seehafer 1990; Pevtsov et al. 1995; Longcope et al. 1998; Nandy 2006), which is a key ingredient that solar dynamo models should be able to reproduce (Charbonneau 2020). The twist in the magnetic field plays an essential part in the dynamics of the solar atmosphere; for example a highly twisted flux tube can become susceptible to kink instability, which leads to a deformation of the axis of the flux tube in exchange for its twist. This is a possible trigger mechanism for solar flares and CMEs (e.g., Török & Kliem 2003; Török et al. 2004; Leka et al. 2005; Fan 2005). Furthermore, the observed twist of photospheric magnetic field is used as an input to inject twist into coronal magnetic field extrapolations (e.g., Yeates et al. 2008; Wiegelmann & Sakurai 2012).

Various methods have been developed to measure the twist density of the magnetic field in active regions directly from individual photospheric observations. These methods either use the force-free parameter, α, as a proxy for the twist density or try to fit the twist density directly.

Woltjer (1958) shows that the force-free parameter α in a closed system corresponds to the helicity content of a linear force-free field structure in its lowest attainable energy state. Therefore, α is used in observations as a proxy for the helicity of the magnetic field (Pevtsov et al. 2014).

As we have observations at only one height in the photosphere from instruments like HMI, we cannot measure α directly. We are limited to calculating the vertical current density J_z and consequently α_z = J_z/B_z at one height, where B_z is the vertical field strength. Burnette et al. (2004) studied 34 active regions and show that a two-dimensional spatial average of α_z over an active region is correlated with the α value corresponding to the three-dimensional linear-force-free extrapolation with the best least-squares fit to the observed field. Therefore, spatial averages of α_z are often used to characterize the helicity and twist of the magnetic field in active regions or individual sunspots (e.g., Longcope et al. 1998; Hagino & Sakurai 2004). Leka et al. (2005) chose a single peak value of α_z close to center the of a sunspot (α_Peak) to characterize the twist in the magnetic field. This is because in simple models α_z only relates directly to the twist density at the axis of a flux tube (α_Peak = 2q).

Nandy et al. (2008) suggested a method to infer the twist density of the magnetic field q that avoids the force-free assumption. These authors assume that the magnetic field in a sunspot can be approximated by a monolithic vertical flux tube with a constant twist density. They then use a least-squares fitting approach to fit the observations to this reference model to obtain the twist density.

Various tests of these methods have been conducted. Leka & Skumanich (1999) and Leka (1999) used observations to compare methods of using moments of the distribution of α_z, a global α from force-free extrapolations, and a fitting approach of the function J_z = αB_z. These authors find quantitative agreement between these methods. They also assessed the influence of instrumental effects like spatial resolution and a limited field-of-view, and considered noise by restricting the α_z measurements to certain noise thresholds. Leka et al. (2005) successfully retrieved the twist density using their α_Peak method on a model by Fan & Gibson (2004) in the absence of errors. Crouch (2012) evaluated different least-squares fitting methods of the twist density from a model flux tube, and find that the inferred twist density can be significantly different depending on the model assumptions used for fitting, also in the absence of noise. Tiwari et al. (2009a) used a linear force-free magnetic field model to test the effect of random polarimetric noise on estimates of the global α value of the synthetic field structure. They find that noise does not influence the sign of α and the global twist can be measured accurately.

To interpret twist density measurements from a single observation, we need to understand how these measurements are affected by temporal variations of the magnetic field. HMI observes the magnetic field in the photosphere where the force-free assumption is thought to be violated (Gary 2001). Temporal fluctuations of the magnetic field may arise in such an environment; for example a twisted magnetic field structure can be distorted by its surrounding plasma flows. We need to model the fluctuations of the magnetic field in a sunspot from SDO/HMI observations to characterize the sensitivity of twist measurements to these fluctuations.

We tested the robustness of existing methods to infer the twist of the magnetic field under temporal variations of the magnetic field from SDO/HMI observations. We modeled the well-established leading sunspot of active region NOAA 11072 (observed by SDO/HMI at 2010.05.25 03:00:00 TAI) with the semi-empirical sunspot model by Cameron et al. (2011) with added uniform twist. We studied the spatial covariance of the temporal fluctuations of the magnetic field and created a model based on our findings. We tested the robustness of methods to measure twist using Monte-Carlo simulations of the sunspot and fluctuation model.

In Sect. 2 we present vector magnetic field observations of the reference sunspot in active region NOAA 11072. Sections 3 and 4 describe the fluctuation and sunspot model, respectively. In Sect. 5 we present a summary of the twist measurement methods that we test in this paper as well as their implementation. We then qualitatively compare our sunspot model to the SDO/HMI observations of the reference sunspot in Sect. 6. In Sect. 7 we present Monte-Carlo simulations to test how the twist measurement methods fare under the influence of temporally fluctuating magnetic field.

2. SDO/HMI vector magnetogram observations of our reference sunspot in active region NOAA 11072

In this section, we present a sunspot observed by SDO/HMI that we selected as a reference for our sunspot model. The reference sunspot should closely resemble the assumption that its magnetic field structure could be described as a monolithic uniformly twisted flux tube. Therefore, we looked for sunspots that are well established, roughly circular, and under little influence from other strong magnetic field in its vicinity. We model the sunspot with uniform twist to test various methods to measure the twist under temporal fluctuations of the magnetic field.

We chose the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI), which was located about 30° away from disk center at the Stonyhurst heliographic coordinates 27° east and 13° south. Figure 1 shows the Postel-projected sunspot transformed to local cylindrical coordinates from SDO/HMI vector magnetogram observations (hmi.b_720s, Hoeksema et al. 2014). B_z is the component normal to the surface. B_r and B_θ are located in a plane parallel to the surface. B_r points radially away from the center of the spot and B_θ is always perpendicular to B_r. See Appendix A for a detailed description of the coordinate systems and transformations used.

Fig. 1. SDO/HMI vector magnetogram of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI) in cylindrical coordinates. Br (top), Bθ (middle), and Bz (bottom) show the radial, azimuthal, and vertical component of the magnetic field.

The sunspots we considered have a dominant radial B_r and weak azimuthal B_θ component, as shown for the example of the reference sunspot in Fig. 2. This is a known characteristic of sunspots (Borrero & Ichimoto 2011). We noticed that the azimuthal component – which carries information about the handedness of the twist in a uniformly twisted flux tube – shows that the sunspot has regions with opposite sign of twist.

Fig. 2. SDO/HMI vector magnetogram of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI) with the horizontal field Bhor plotted as arrows on top of the vertical vector component Bz.

3. Estimating the spatial covariance of magnetic field fluctuations from the observations

We aimed to derive a model for the temporal fluctuations of the magnetic field in SDO/HMI vector magnetogram observations. To do so, we used an approximately seven-hour time-series of our reference spot (2010.05.25, 03:00:00−9:48:00 TAI) to look at each local Cartesian vector component (B_x, B_y, B_z) in Postel-projected maps individually. The observations have a cadence of 12 min. B_z is the vector component normal to the surface, B_x and B_y point from solar east to west and south to north, respectively. The time-frame was chosen so that the spot is stable.

We tracked the proper motion of the sunspot by first calculating the flux-weighted centroid of |B_z| within the sunspot in each observation. We defined the area for this calculation based on pixels with a value below 0.85 in normalized continuum maps of the sunspot. We found that the sunspot moved approximately three pixels over this time-period almost linearly. We fit a line to the location of the centroid in the x- and y- directions, and used the fit to shift the centroids of each image to the same location.

We then detrended the time-series of each pixel by fitting a third-order polynomial to the data and keeping the residuals (sketched in Fig. 3) to remove any long-term trends. We calculated the spatial correlation of these detrended time-series using the Pearson correlation coefficient. Figure 4 shows for each vector magnetic field component the average correlation within the sunspot of a pixel relative to its neighbors. We find that, on average, pixels in the observations are correlated with their neighbors up to 3−4 pixels away.

Fig. 3. Sketch of the detrending process. A series of consecutive vector magnetograms (e.g., the Bz component) labeled with different time-steps (t1, t2, t3, ...), is shown on the left. The black line (data) in the top right panel represents the temporal evolution of one pixel in this time-series (marked in red in the left panel). The blue line is a third-order polynomial fit to the data. Bottom right panel: detrended time-series, which shows the residuals of the data with respect to the fit.

Fig. 4. Average spatial correlation of the detrended time-series of a pixel within the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00−2010.05.25 9:38:00 TAI) relative to its neighbors for the magnetic field vector components Bx, By, and Bz.

In Appendix B, we consider whether the observed correlations were due to the Postel projection, the detrending method, the 12 and 24 h periodicity of HMI caused by the satellites changing radial velocity relative to the Sun (Hoeksema et al. 2014), or HMI’s point spread function (PSF). We concluded that these correlations are caused by the PSF and the dynamic changes of the magnetic field over time.

We incorporated the information about these time variations of the magnetic field into our model. Knowing that adjacent pixels are correlated but the strength of correlation depends on their position within the spot, we calculated the covariance matrix for all detrended time-series and for each vector component individually. We used the Cholesky decomposition (Haddad 2009) of these covariance matrices to create random correlated maps that reflect spatially correlated temporal fluctuations of the magnetic field. Figure 5 shows example realizations of such maps for each vector component.

Fig. 5. Each panel represents one realization of magnetic field fluctuations for each vector component. All panels are plotted on the same scale in units of Gauss displayed by the color bar at the top. The black solid line represents the observed sunspot boundary.

4. A model for the reference sunspot

In this section we present the semi-empirical sunspot model developed by (Cameron et al. 2011) modified to represent a uniformly twisted field structure. Cameron et al. (2011) describes a three-dimensional magnetic field model of an axisymmetric sunspot with a radial B_r and vertical B_z component. It lacks the azimuthal component B_θ, which is essential for creating twist. We added a B_θ component that is only dependent on the radius without violating the requirement of ∇ ⋅ B = 0 in 3D space. We were only interested in the magnetic field structure on the photospheric level (z = 0) to model a sunspot observed by SDO/HMI. Specifically, the magnetic field at the photosphere is

$$\begin{array}{cc}& B_z(r)=B_{0}exp[-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2],\hfill \end{array}$$(2)

$$\begin{array}{cc}& B_{\mathrm{r}}(r)=\frac{r{B}_z(r)}{8a_0\sqrt{1+b^2}},\hfill \end{array}$$(3)

$$\begin{array}{cc}& B_{\theta }(r)=b{B}_{\mathrm{r}}(r),\hfill \end{array}$$(4)

where B₀ is the magnetic field strength at the center of the spot, r is the distance from that center, and h₀ defines the radius of the umbra–penumbra boundary. The parameter a₀ controls the inclination of the field and b governs the amount of twist in the model. The model of Cameron et al. (2011) is designed so that the inclination of the magnetic field at the umbra–penumbra boundary is 45°. Due to the additional azimuthal component B_θ and to keep the same inclination profile we adjusted the parameter controlling the inclination a₀ in accordance with the injected twist by multiplying it with $\sqrt{1+b^2}$.

The pixel scale of our model is the same as HMI’s pixel scale of 0.5″, which corresponds to approximately 0.35 Mm at disk center. We fit the four free parameters (B₀, h₀, a₀, b) to the reference spot. We use a least-squares fitting approach to best match the azimuthal averages of B_z, B_r, and B_θ around the flux-weighted center of the reference sunspots |B_z|.

The information about the twist of the spot is stored in B_θ. As shown in Figs. 1 and 2, B_θ does not show an symmetric behavior about the center of the spot and azimuthal averages do not represent the local twist present in the observation. We find that fitting only the positive or negative values of B_θ yields values of 0.125 and −0.175 for the parameter b, respectively. As the magnetic field twist of the reference sunspot has a preference to be negative, we chose b = −0.15.

The parameters that we find to best describe the reference sunspot with uniform twist are B₀ = 2370 G, h₀ = 5.2 Mm, a₀ = 0.77 Mm, and b = −0.15.

5. Summary and implementation of twist measurement methods

Now we present various methods proposed in the literature to estimate the twist density directly from a single photospheric observation and describe their numerical implementation.

5.1. The twist proxy α

The force-free parameter α can serve as a helicity proxy (Pevtsov et al. 2014) to estimate the twist in uniformly twisted flux tubes (see Appendix F for an interpretation of α in terms of twist). We lack information in photospheric observations of how the B_x and B_y components of the magnetic field change as a function of height (z-direction in a local Cartesian coordinate system) and we can only compute the vertical current density J_z (e.g., Pevtsov et al. 1994; Longcope et al. 1998).

This vertical current density can be calculated from the force-free equation,

$$\begin{array}{c}\hfill \mathbf{\nabla }\times \mathit{B}=\mathit{J}=\alpha \mathit{B},\end{array}$$(5)

with

$$\begin{array}{c}\hfill J_z=\frac{\partial B_y}{\partial x}-\frac{\partial B_x}{\partial y}·\end{array}$$(6)

Consequently we can compute the local twist proxy α_z at a specific location with

$$\begin{array}{c}\hfill {\alpha }_z=\frac{J_z}{B_z}·\end{array}$$(7)

We note that α is a pseudo-scalar and the subscript z denotes that it was derived only from the vertical field component and vertical current density. Positive (negative) values of α_z correspond to right-handed (left-handed) magnetic field twist, respectively.

Numerically, we calculated the derivatives using the Savitzky-Golay filter (Savitzky & Golay 1964, Appendix C) of cubic and quartic order and a stencil size of five pixels. We note that J_z can also be calculated in integral form using Stokes’ theorem (see Appendix D). The stencil size was chosen based on our own tests on how the stencil size impacts spacial averages of α_z (see Appendix E) and the results by Fursyak (2018).

5.2. Average twist

Pevtsov et al. (1995) used a single best fit value of α from linear force-free field extrapolations to characterize the twist for whole active regions. Longcope et al. (1998) used an average ⟨α_z⟩ for entire active regions, which can be calculated from photospheric observations without the need of any extrapolations:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{av}}=⟨{\alpha }_z⟩=\langle \frac{J_z}{B_z}\rangle .\end{array}$$(8)

Hagino & Sakurai (2004) proposed two weighted averages of α_z over a whole active region or spot to determine its twist:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{av}}^{\mathrm{abs}}=\frac{\sum J_z\mathrm{sign}\left[B_z\right]}{\sum \left|B_z\right|}\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{av}}^{\mathrm{sqr}}=\frac{\sum J_z{B}_z}{\sum B_z^2}·\end{array}$$(10)

${\alpha }_{\mathrm{av}}^{\mathrm{abs}}$ and ${\alpha }_{\mathrm{av}}^{\mathrm{sqr}}$ are weighted by absolute and squared B_z, respectively, which is represented by the superscripts “abs” and “sqr”. These latter authors argue that these weighted averages have the advantage of putting less weight on weak field, especially close to the polarity inversion line, where singularities of α_z are more likely to occur.

The area over which α_z is averaged depends on the area of interest that is studied. In Sect. 7 we investigate averages over a small central area of the spot, the umbra, and the whole spot.

5.3. Peak twist

Under the assumption that a spot can be approximated by a monolithic uniformly twisted magnetic flux tube, Leka et al. (2005) show that the α_z profile of this field structure has a peak directly at the center of the spot (flux tube axis), which they name α_Peak. Based on the flux tube model by Gold & Hoyle (1960), Leka et al. (2005) demonstrate that α_Peak directly relates to the constant twist density q of the field’s structure (α_Peak = 2q).

In order to calculate α_Peak for a single sunspot, we estimated the location of the flux tube axis by computing the flux-weighted center |B_z| of the sunspot. We calculated a map of α_z values (Eq. (7)) for each pixel within the spot and boxcar-smoothed this map to 2″ as suggested by Leka et al. (2005). We used the absolute values of this smoothed map to detect the α_z-peak closest to the estimated flux tube axis. α_Peak is then the signed and smoothed α_z value at the location of the peak.

5.4. Twist density

Nandy et al. (2008) proposed to fit the twist density q to quantify the magnetic twist in a single spot. Again, under the assumption that the magnetic field in a sunspot resembles a uniformly twisted flux tube, q can be measured by fitting the slope of the equation

$$\begin{array}{c}\hfill \frac{B_{\theta }}{B_z}=qr+d,\end{array}$$(11)

where r is the distance from the flux tube axis, and B_θ and B_z are the magnetic field in azimuthal direction and along the axis of the tube, respectively. The axis of the tube is estimated by calculating the flux-weighted center of |B_z| of the sunspot. We note that Nandy et al. (2008) allowed a nonzero intercept d to occur in their Fig. 2. This violates the assumption of a vertical uniformly twisted flux tube, where the fitted function is expected to go through zero at r = 0, which is equivalent to a vanishing B_θ at the axis of the flux tube. A physical interpretation of this intercept is not clear to us.

This method was carried out in a cylindrical coordinate system (B_r, B_θ, B_z, see Appendix A). We fit the ratio B_θ over B_z for each pixel as a function of the distance between the pixel and the estimated flux tube axis r. The resulting slope corresponds to the twist density q. We tested this method by both allowing an intercept and by forcing the fit through the origin (d = 0).

6. Example sunspot model with correlated magnetic field fluctuations compared to HMI observations

In this section we qualitatively compare the magnetic field and its twist between our sunspot model with one realization of magnetic field fluctuations and the SDO/HMI observations of the reference sunspot in active region NOAA 11072.

6.1. The vector magnetic field

Figure 6 shows that visually the model spot with correlated magnetic field fluctuations resembles the HMI observation of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI). It even exhibits a filamentary structure especially noticeable in B_z, which is not present in the model without fluctuations of the magnetic field. The fluctuation maps (Fig. 5) have smoother and weaker fluctuations within the sunspot, and stronger, more chaotic fluctuations in the quiet Sun, as one would expect from the observations.

Fig. 6. Comparison of the vector components (Bx, By, Bz) of the magnetic field. First row: HMI observation of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI), which serves as the reference for our sunspot model. Second and third rows: model with and without temporal fluctuations, respectively. Every column is plotted on the same scale in units of Gauss displayed by the color bar at the top of the column.

6.2. The magnetic field’s twist

Figure 7 compares α_z-profiles of the original SDO/HMI observation of the reference sunspot and our model with and without one realization of temporal fluctuations of the magnetic field. Our sunspot model without fluctuations describes a uniformly left-handedly twisted magnetic field structure. The α_z-profile is azimuthally symmetric and α_z increases with distance from the centers of the spots. The sign of α_z changes in the penumbra, which signals the presence of return currents (see Appendix F).

Fig. 7. Comparison of αz maps measured from the original HMI observation (top), and from the model without (middle) and with (bottom) fluctuations of the magnetic field. The color bar at the top shows the αz values in Mm−1.

The reference sunspot has a more complicated structure than the model without temporal fluctuations. Even in the center of the spot, areas of opposite sign of α_z exist. Towards the penumbra, positive values of α_z become more frequent and one could assume a ring of return currents similar to the model. After applying magnetic field fluctuations to the model, a similar structure of the α_z pattern compared to the observations develops.

Our definition of fluctuations includes dynamic variations of the magnetic field. Correlated changes in the direction of magnetic field in adjacent pixels can produce spatially coherent changes in twist and its sign in our model. Sunspot observations typically show a strong radial field component and exhibit only weak twist (i.e. a weak azimuthal field component, B_θ ≪ B_r). A source of fluctuations of the magnetic field is the forced environment of the photosphere, where the magnetic field can be buffeted by plasma flows, which can cause sign changes of the real twist. Such an effect is expected to be stronger where the magnetic field strength is weaker – for example in the penumbral parts of the sunspot – where we see the strongest variations of α_z within sunspots. Also, interactions with magnetic flux surrounding a sunspot could cause deviations from a uniformly twisted field structure. Complex patterns of α_z and magnetic field twist even within the umbra of sunspots have been described in the literature (e.g., Pevtsov et al. 1994; Socas-Navarro 2005; Su et al. 2009).

Figures 8 and 9 show the temporal evolution of α_z and its sign for the leading sunspot of active region NOAA 11072, respectively. In Fig. 8 we find that in most parts of the umbra, α_z is on the order of 10⁻² Mm⁻¹, while penumbral α_z values are typically at least one order of magnitude larger, which is consistent with other α_z-measurements in sunspots (e.g., Tiwari et al. 2009b; Wang et al. 2021).

Fig. 8. Temporal evolution of αz in the leading sunspot of active region NOAA 11072. The time relative to the first observation (2010.05.25 03:00:00 TAI) is shown in hours in the top left of each panel. The black solid line represents the umbra–penumbra boundary.

Fig. 9. Temporal evolution of the sign of αz in the leading sunspot of active region NOAA 11072. The time relative to the first observation (2010.05.25 03:00:00 TAI) is shown in hours in the top left of each panel. The black solid line represents the umbra–penumbra boundary.

We find patches of α_z with opposite signs throughout the reference sunspot, which is consistent with previous studies of sunspots (e.g., Pevtsov et al. 1994; Socas-Navarro 2005; Su et al. 2009). The shape of these patches are in agreement with findings by Su et al. (2009), who describe a mesh-like pattern in the umbra and a thread-like pattern in the penumbra. Figure 9 shows that many of these patches persist over timescales of hours.

We find a similar distribution of α_z values and patterns in the sunspot model with a realization of magnetic field fluctuations. The uncertainty on the mean α_z (standard deviation of α_z divided by the number of pixels considered) within the model sunspot with fluctuations is of the order of 10⁻³ Mm⁻¹. This is of the same order as the uncertainty Leka & Skumanich (1999) and Leka (1999) measured in active regions observed with the Image Vector Magnetograph at Mees Solar Observatory.

Features in the pattern structure appear to be on a smaller scale in our model. Pevtsov et al. (1994) studied the helical structure of the magnetic field of three active regions and estimated that the lifetime of such patches can exceed a day. Our model only describes the spatial correlation of magnetic field fluctuations but does not address their temporal correlation. Therefore, these patches appear uncorrelated from one realization to another. Whether or not the temporal variations of the magnetic field are responsible for the twist and current patterns that we can observe in real sunspots cannot be deciphered yet from our model. To further investigate this question, one also has to consider the temporal correlation of the magnetic fluctuations.

7. Sensitivity of twist measurements to correlated temporal fluctuations of the magnetic field

We used Monte-Carlo simulations to test the sensitivity of twist measurement methods described in Sect. 5 to fluctuations of the magnetic field. We used our magnetic field fluctuation model (described in Sect. 3) in 10 000 realisations to create different fluctuation maps and superimposed them on the sunspot model (Sect. 4). We evaluated in each iteration the different twist proxies in the umbra (up to r = h₀). We also tested the robustness of the twist measurement methods based on the area over which they are averaged. We evaluated α_av in a small umbral area with a radius of 1.75 Mm from the center of the spot (${\alpha }_{\mathrm{av}}^{\mathrm{center}}$, white circle in Fig. 10) and up to the penumbra–quiet Sun boundary (${\alpha }_{\mathrm{av}}^{\mathrm{spot}}$, black dotted circle in Fig. 10, up to r = 2h₀). Fluctuations of the magnetic field in our model can result in vertical field B_z close to zero in the penumbra, which can create singularities when calculating α_z = J_z/B_z. Therefore, we only measure α_z for pixels that are above a B_z threshold of 50 Gauss. The analytically calculated reference values of α_z and q that represent the uniform twist of our model best (see Appendix F) are ${\alpha }_{\mathrm{Peak}}^{\mathrm{ref}}=2q\approx -0.048$ Mm⁻¹.

Fig. 10. HMI continuum image of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI), which was used as a reference for the model presented in this work. The black solid line outlines the umbral area that was considered for most twist calculation methods. The white solid line and the black dotted line correspond to the areas that were used to get spatial averages of αz close to the center of the spot (${\alpha }_{\mathrm{av}}^{\mathrm{center}}$) and over the whole spot (${\alpha }_{\mathrm{av}}^{\mathrm{spot}}$), respectively.

Figures 11 and 12 show the distribution of the calculated twist proxies from the Monte-Carlo simulations for each method. Table 1 compares the mean result and standard deviation from the Monte-Carlo simulations against the expected value from the model without fluctuations. It is important to note that the “Model” values in Table 1 are derived when a method is applied in the fluctuation-free model. The errors given in Table 1 show the standard deviation of the Monte-Carlo simulations. We find that the expectation values of the averaging methods and the twist density fits are not biased by magnetic field fluctuations.

Fig. 11. Monte-Carlo simulation results for methods described in Sect. 5. The q values have been multiplied by two in order to make them directly comparable to αz (see Appendix F). The black line indicates the reference value of the model, which is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the Monte-Carlo mean.

Fig. 12. Monte-Carlo simulation results when αz is spatially averaged over the whole sunspot. The solid black line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the mean.

The averaging methods of α_z show a large spread in their results, but have robust measurements under magnetic field fluctuations. These different spatial averages of α_z can still be related to the twist density based on the azimuthal symmetric behavior of α_z in our simple model. We derive ${\alpha }_z(r)=2q[1-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2]$ in Appendix F. We can calculate the average value of α_z in a circular area with radius r around the center of the spot:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {⟨{\alpha }_z⟩}_{\mathrm{r}}=\frac{1}{\pi r^2}{\displaystyle {\int }_0^r{\int }_0^{2\pi }{\alpha }_z(\stackrel{\sim }{r})\stackrel{\sim }{r}\mathrm{d}\stackrel{\sim }{r}\mathrm{d}\theta =[2-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2]q.}\end{array}\end{array}$$(12)

Consequently, we find the following relation between the twist density q and ⟨α_z⟩_r:

$$\begin{array}{c}\hfill q=\frac{{⟨{\alpha }_z⟩}_{\mathrm{r}}}{2-\left({log}_{\mathrm{e}}2\right){(r/h_0)}^2}·\end{array}$$(13)

This equation is consistent with the relation of α_Peak = 2q at the center of the spot. It also describes the different spatial averages of α_z that we measure, when the averaging radius was changed (see Figs. 11 and 12). We measure a larger ⟨α_z⟩ when averaged over a small area at the center of the spot (${\alpha }_{\mathrm{av}}^{\mathrm{center}}$) compared to the average over the umbra (α_av). When the averaging radius becomes sufficiently large (e.g., ${\alpha }_{\mathrm{av}}^{\mathrm{spot}}$), ⟨α_z⟩ retrieves the opposite sign compared to α_z in the center of the spot.

The expectation value of the α_Peak method is noticeably biased, causing an overestimation of the twist density in the presence of magnetic field fluctuations. As this method characterizes the twist density with a single peak value closest to the center of the spot, it is likely to pick up any enhanced signal caused by the fluctuations of the magnetic field.

The twist density fit (Nandy et al. 2008) retrieves the analytical twist density value. We find a larger spread in the resulting values when no zero-intercept is forced. These results are expected, because the equations of our sunspot model can be exactly reduced to the fitting equation by Nandy et al. (2008). Crouch (2012) shows that fitting techniques can be sensitive to small discrepancies between the fitting and reference model. Observations suggest that sunspots are more complicated and can have nonuniformly twisted field structures (e.g., Socas-Navarro 2005; Su et al. 2009). Another complication for real sunspots is that this method requires the exact location of the axis of a flux tube as a center for the coordinate transformation to cylindrical coordinates. This is a simple task in our model without fluctuations, because the central axis of the model spot and its flux-weighted center fall into the same place by definition; even with magnetic field fluctuations, they are always located close to each other. In real sunspots, the axis of the spot does not have to be in the same place as its flux-weighted center. Also, a sunspot may not have an underlying uniformly twisted vertical field structure.

In Sect. 4 we chose the model parameter b = −0.15, which governs the twist in our model. We tested various values of b, ranging from untwisted field (b = 0) to highly twisted field (b = 10) and found that the amount of twist in the model does not influence the findings about the robustness of the twist measurements described in this section.

8. Summary and conclusion

We derived a model for the spatially correlated fluctuations of the magnetic field in a sunspot based on HMI observations of the leading sunspot of the active region NOAA 11072. We superposed realizations of the fluctuations on the magnetic field of the semi-empirical sunspot model described in Cameron et al. (2011) with added uniform twist. We carried out Monte-Carlo simulations to test the robustness of the different measures of the twist to the fluctuations.

We considered measurement methods based on estimating either the force-free parameter α_z or the twist density of the magnetic field q. In the absence of fluctuations, and for the sunspot model used in this paper, the value of α_z at the center of the sunspot is twice the twist parameter q (see Leka et al. 2005). For our chosen twist profile, α_z is not uniform and changes sign in the penumbra of the model spot.

Including the spatially correlated temporal fluctuations of the magnetic field qualitatively reproduces features seen in vertical current density J_z and α_z observations. Patches with opposite sign of α_z appear throughout the sunspot model at random locations from one realization to another. Although we do not consider temporal correlations of the fluctuations, we note that such features can persist for hours in observations (Pevtsov et al. 1994).

All measures except α_Peak have expectation values consistent with that of the model without fluctuations. The fluctuations do not introduce a bias. The measures based on spatial averages of α_z have different expectation values, because they average over different portions of the nonuniform α_z profile. Due to the sign change of α_z in the penumbra of our model, any spatial averages of α_z that reach too far out from the center of the spot can have the opposite sign from α_z near the center of the spot. The expectation value of α_Peak is biased with respect to the model without fluctuations. The magnitude of the bias is related to the level of the magnetic field fluctuations.

For most methods, the spread of results from the Monte-Carlo simulations is less than the spread in expectation values of the individual methods. Therefore, the choice of method is more significant than the impact of the magnetic field fluctuations on the measurements.

Our results are for the particular sunspot model given by Eqs. (2)–(4). This is a particularly simple axisymmetric sunspot model. For sunspots with a more complex structure, for example nonuniformly twisted sunspots, the applicability and meaning of the different measures needs to be carefully considered (e.g., Crouch 2012). In this regard, we note that the SDO/HMI observations of the leading spot of active region NOAA 11072 had fine structure in α_z which persisted for longer than 7 h. This persistent fine structure suggests a more complicated underlying field structure than that of our uniformly twisted model. The long-lived fine structure found in active region NOAA 11072 is consistent with previous studies of other sunspots (e.g., Pevtsov et al. 1994; Su et al. 2009).

As was previously noted by Leka (1999), a single parameter will not in general characterize the twist of a sunspot. In this paper, we show that a range of different complementary parameters exists, all of which describe somewhat different aspects of the magnetic field line twist in sunspots. We show that most of these measures are robust to fluctuations of the field.

Acknowledgments

C.B. is a member of the International Max Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. C.B. conducted the analysis, contributed to the interpretation of the results, and wrote the manuscript. We thank Jesper Schou for helpful discussions. We thank Graham Barnes for useful comments on the manuscript. The HMI data used are courtesy of NASA/SDO and the HMI Science Team. A.C.B., R.H.C. and L.G. acknowledge partial support from the European Research Council Synergy Grant WHOLE SUN #810218. The data were processed at the German Data Center for SDO, funded by the German Aerospace Center under grant DLR50OL1701. This research made use of the Astropy (Astropy Collaboration 2013, 2018), Matplotlib (Hunter 2007), NumPy (Van Der Walt et al. 2011) and SciPy (Virtanen et al. 2020) Python packages.

References

Appendix A: Vector transformation

A.1. Local Cartesian coordinates

The hmi.b_720s (Hoeksema et al. 2014) provides the vector magnetic field in spherical coordinates aligned with the line of sight (LoS); it includes the absolute field strength B, the inclination angle “inc” with respect to the LoS, and the azimuth angle “azi” (whose ambiguity has to be resolved) measured in a plane perpendicular to the LoS. We study the magnetic field in a local Cartesian and cylindrical coordinate system, where B_z is pointing radially outwards from the Sun. In order to transform the spherical coordinate system to a local Cartesian coordinate system, we followed the transformations described in Gary & Hagyard (1990).

First, the spherical vector components were transformed to a Cartesian system, where B_ζ is aligned with the LoS, while B_ξ and B_η lay in the plane perpendicular to the LoS:

$$\begin{array}{cc}\hfill B_{\xi }& =-Bsin(\mathrm{inc})sin(\mathrm{azi})\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill B_{\eta }& =Bsin(\mathrm{inc})cos(\mathrm{azi})\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill B_{\zeta }& =Bcos(\mathrm{inc}).\hfill \end{array}$$(A.3)

We then rotated the coordinate system (Eq. A.4) so that B_ζ points radially outward from the Sun and becomes B_z. Here, B_x and B_y point from solar east to west and south to north, respectively. The rotation matrix is

$$\begin{array}{c}\hfill \left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]=\left[\begin{array}{c}{a}_{11}{a}_{12}{a}_{13}\\ a_{21}{a}_{22}{a}_{23}\\ a_{31}{a}_{32}{a}_{33}\end{array}\right]\left[\begin{array}{c}{B}_{\xi }\\ B_{\eta }\\ B_{\zeta }\end{array}\right],\end{array}$$(A.4)

with the transformation matrix coefficients a_ij:

$$\begin{array}{cc}\hfill a_{11}=& -sin{B}_{0}sinPsin(L-L_0)+cosPcos(L-L_0)\hfill \\ \hfill a_{12}=& +sin{B}_{0}cosPsin(L-L_0)+sinPcos(L-L_0)\hfill \\ \hfill a_{13}=& -cos{B}_{0}sin(L-L_0)\hfill \\ \hfill a_{21}=& -sinB[sin{B}_{0}sinPcos(L-L_0)+cosPsin(L-L_0)]\hfill \\ \hfill & -cosB[cos{B}_{0}sinP]\hfill \\ \hfill a_{22}=& +sinB[sin{B}_{0}cosPcos(L-L_0)-sinPsin(L-L_0)]\hfill \\ \hfill a_{23}=& -cos{B}_{0}sinBcos(L-L_0)+sin{B}_{0}cosB\hfill \\ \hfill a_{31}=& +cosB[sin{B}_{0}sinPcos(L-L_0)+cosPsin(L-L_0)]\hfill \\ \hfill & -sinB[cos{B}_{0}sinP]\hfill \\ \hfill & +cosB[cos{B}_{0}cosP]\hfill \\ \hfill a_{32}=& -cosB[sin{B}_{0}cosPcos(L-L_0)+sinPsin(L-L_0)]\hfill \\ \hfill & +sinB[cos{B}_{0}cosP]\hfill \\ \hfill a_{33}=& +cosBcos{B}_{0}cos(L-L_0)+sinBsin{B}_0.\hfill \end{array}$$

Here, L and B describe the heliographic longitude and latitude of the individual pixels, while L₀ and B₀ are the longitude and latitude of the solar disk center, respectively, and P is the solar position angle.

A.2. Cylindrical coordinates

We transformed from the local Cartesian coordinate system to a local cylindrical coordinate system by calculating

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left[\begin{array}{c}{B}_r\\ B_{\theta }\\ B_z\end{array}\right]=\left[\begin{array}{cccc}& cos\theta & sin\theta & 0\\ & -sin\theta & cos\theta & 0\\ & 0& 0& 1\end{array}\right]\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right],\end{array}\end{array}$$(A.5)

with θ = arctan2(y, x). The flux-weighted center of a spot is defined as the origin (x = 0, y = 0) for this transformation. In this coordinate system, B_z is the component normal to the surface; B_r and B_θ are located in a plane parallel to the surface; B_r points radially away from the center of the spot; and B_θ is always perpendicular to B_r and B_z.

Appendix B: Possible causes for the measured spatial correlation of magnetic field fluctuations

We tested various effects that could introduce the spatial correlation of temporal fluctuations of neighboring pixels. To assess the effect of the Postel projections on these correlations, we created artificial full-disk HMI maps. We filled pixels with random Gaussian white noise with a mean of zero and standard deviation of one. We created Postel-projected time-series of submaps at different locations on the solar disk based on the center to limb angle (CTL) of the center of the submap. We followed the same detrending procedure as described in section 3, but find no correlation of adjacent pixels either at disk center (CTL = 0°) or close to the limb (CTL = 60°).

We compared different ways of detrending the data: different order of polynomials for fitting (order 3, 4, and 5), differences to previous data points, and running averages. All different detrending methods show similar correlations of adjacent pixels. We find no relation of the detrended signals to the known 12- and 24-hour period systematic errors of HMI that are caused by the satellite’s orbit and change in radial velocity relative to the Sun.

The PSF of HMI definitely contributes to the correlations of adjacent pixels. Figure B.1 shows an estimate of the PSF of HMI by Yeo et al. (2014). Figure B.2 shows cuts along the x- and y-axis through the center of the normalized PSF, and the average correlations of adjacent pixels are shown in Fig. 4. We find that the fluctuations of the vector magnetic field components are correlated up to spatial scales that are 30% larger than that expected from the HMI PSF. We suggest that these correlations are caused by the PSF and the dynamic changes of the magnetic field over time with respect to the underlying global field structure of the sunspot.

Fig. B.1. Normalized PSF for HMI estimated by Yeo et al. (2014).

Fig. B.2. Cuts through the center along the x-axis (top panel) and the y-axis (bottom panel) of the average spatial correlation of neighboring pixels (Fig. 4) and the normalized estimated PSF by Yeo et al. (2014) (Fig. B.1). The errorbars represent the standard error of the spatially averaged correlation.

Appendix C: Derivatives with different stencil sizes

We can evaluate derivatives with different stencil sizes using Savitzky-Golay filter (Savitzky & Golay 1964) of cubic/quartic order for stencil sizes of 5 and 7 pixels or take a central differences approach for a stencil of 3 pixels:

$$\begin{array}{c}\hfill \frac{\partial f(x)}{\partial x}=\frac{1}{\mathit{Nh}}{\displaystyle \sum _{i=-(S-1)/2}^{(S-1)/2}}g(i)f(x+i),\end{array}$$(C.1)

where x describes the location of a pixel counted as integer, for which we want to calculate the derivative. S is the stencil size, N is a normalization factor, h is the pixel scale, and g is a weighting factor. The weighting and normalization factor for the different stencil sizes are listed in Table C.1.

Appendix D: Using Stokes’ theorem for calculating J_z

Instead of using derivatives for calculating J_z (Eq. 6), one can also use Stokes’ theorem:

$$\begin{array}{c}\hfill J_z^{\mathrm{Int}.}=\frac{1}{{\mathrm{A}}_{\mathrm{L}}}{\displaystyle {\oint }_L}{\mathbf{B}}_{\mathrm{hor}}·d\mathbf{L},\end{array}$$(D.1)

where the vertical current density is calculated for a pixel in the center of an area A_L outlined by a contour L. We calculated the integral by using the composite Simpson’s rule over a square, where L is the edge of the square with a side length of 3, 5, or 7 pixels, in accordance to the stencil size when using derivatives.

Using Stokes’ theorem for calculating α_z leads on average to slightly lower values compared to the differential form. This effect can be attributed to the integral form using more pixels for the calculations and averaging over an area. The results of Monte-Carlo simulations are shown in Fig. D.1 and Fig. D.2 as well as Table D.1.

Fig. D.1. Monte-Carlo simulation results for methods based on αz when Stokes’ theorem is used for calculating Jz. The black solid line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The dark/light gray shaded areas represent the range of 1, 2, or 3σ around the Monte-Carlo mean.

Fig. D.2. Monte-Carlo simulation results when αz is calculated with Stokes’ theorem and spatially averaged over the whole sunspot. The black solid line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the mean.

Appendix E: Comparison of different stencil sizes for calculating J_z

We tested stencil sizes of 3, 5, and 7 pixels in Monte-Carlos simulations as described in section 7 to determine how many pixels should be considered for calculating the vertical current density J_z and subsequently α_z. The average twist proxy values and their standard deviation from these simulations are listed in Tables E.1 and E.2. The different stencil sizes do not noticeably impact the end result when we used the Savitzky-Golay filter. In contrast, we find that with bigger areas the measured twist goes down when Stokes’ theorem (see Appendix D) was applied. Similar to the Monte-Carlo simulations described in section 7 only the expectation value of the α_Peak method is biased in all test cases.

Fursyak (2018) tested the effect of differently sized areas for calculating J_z in differential and integral forms using observations from Hinode and HMI. They conclude that integrating over side lengths of five pixels provides the best compromise between smoothing noise but still preserving significant features. As the results from our tests do not show a clear favorite, we chose to calculate J_z with stencil sizes of 5 pixels in accordance to Fursyak (2018).

Appendix F: Interpretation of α_z

It is often difficult to interpret α_z measurements. Ideally we would expect to measure the same α_z or twist density q value at each location within a uniformly twisted flux tube (which is the case in thin flux tube models). While this is true for the twist density q in our sunspot model without magnetic field fluctuations, α_z varies with the distance from the center of the spot, where the α_z profile peaks.

Leka et al. (2005) used the Gold-Hoyle flux tube model (Gold & Hoyle 1960) to show the connection between α_z and the twist density. This model consists of an axial and azimuthal magnetic field component. Leka et al. (2005) demonstrate that only at the center of the flux tube, where the radial distance from its axis r equals zero, the twist density can be retrieved from α_z measurements directly. They describe the following relationship between α_z, the twist density q and the distance from the flux tube axis as:

$$\begin{array}{c}\hfill {\alpha }_z(r)=\frac{2q}{1+q^2{r}^2}·\end{array}$$(F.1)

We note that in this configuration α_z goes towards zero when r increases to infinity. This finding led them to propose the α_Peak method. We can find a similar relationship between α_Peak and the twist density q for our empirical sunspot model. Taking the model’s equation for B_θ (Eq. 4) with $q=\frac{b}{8a_0\sqrt{1+b^2}}$, we get the twist density equation from Nandy et al. (2008):

$$\begin{array}{c}\hfill B_{\theta }(r)=qr{B}_z.\end{array}$$(F.2)

Together with the model’s B_z component,

$$\begin{array}{c}\hfill B_z(r)=B_{0}exp[-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2],\end{array}$$(F.3)

we can use derivatives in cylindrical coordinates to calculate α_z. We note that B_r in Eq. 3 does not depend on the azimuthal angle θ. Then the radial profile of α_z is:

$$\begin{array}{cc}\hfill {\alpha }_z& =\frac{1}{r{B}_z}[\frac{\partial (r{B}_{\theta })}{\partial r}-\frac{\partial B_r}{\partial \theta }]\hfill \\ \hfill & =\frac{1}{r{B}_z}\frac{\partial }{\partial r}\{q{B}_0{r}^{2}exp[-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2]\}\hfill \\ \hfill & =2q[1-\left({log}_{\mathrm{e}}2\right){\left(\frac{r}{h_0}\right)}^2].\hfill \end{array}$$(F.4)

We get a maximum value of α_z (α_Peak) at the center of the flux tube, i.e. r = 0. Therefore, we find the same relationship of α_Peak = 2q as Leka et al. (2005).

In our empirical sunspot model α_z changes sign at a certain distance from the suncenter of the spot. We can determine where alpha changes sign by setting α_z in Eq. F.4 equal to zero to get

$$\begin{array}{c}\hfill r_{{\alpha }_z=0}=\frac{h_0}{\sqrt{{log}_{\mathrm{e}}2}}\approx 1.2h_0.\end{array}$$(F.5)

The parameter h₀ describes the distance of the umbra–penumbra boundary measured from the model center of the spot. Disregarding fluctuations of the magnetic field, within the umbra the sign of the twist and α_z are the same, according to our model.

Since α_z = J_z/B_z, such a sign reversal can only happen in our sunspot model with positive polarity, if the sign of the vertical current density J_z changes. Such rings of return currents at the umbra–penumbra boundary of observed sunspots are reported in literature (e.g., Tiwari et al. 2009b). We note that positive and negative vertical currents are balanced in our model.

All Tables

All Figures

	Fig. 1. SDO/HMI vector magnetogram of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI) in cylindrical coordinates. Br (top), Bθ (middle), and Bz (bottom) show the radial, azimuthal, and vertical component of the magnetic field.
In the text

	Fig. 2. SDO/HMI vector magnetogram of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI) with the horizontal field Bhor plotted as arrows on top of the vertical vector component Bz.
In the text

Fig. 3. Sketch of the detrending process. A series of consecutive vector magnetograms (e.g., the Bz component) labeled with different time-steps (t1, t2, t3, ...), is shown on the left. The black line (data) in the top right panel represents the temporal evolution of one pixel in this time-series (marked in red in the left panel). The blue line is a third-order polynomial fit to the data. Bottom right panel: detrended time-series, which shows the residuals of the data with respect to the fit.

In the text

	Fig. 4. Average spatial correlation of the detrended time-series of a pixel within the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00−2010.05.25 9:38:00 TAI) relative to its neighbors for the magnetic field vector components Bx, By, and Bz.
In the text

	Fig. 5. Each panel represents one realization of magnetic field fluctuations for each vector component. All panels are plotted on the same scale in units of Gauss displayed by the color bar at the top. The black solid line represents the observed sunspot boundary.
In the text

Fig. 6. Comparison of the vector components (Bx, By, Bz) of the magnetic field. First row: HMI observation of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI), which serves as the reference for our sunspot model. Second and third rows: model with and without temporal fluctuations, respectively. Every column is plotted on the same scale in units of Gauss displayed by the color bar at the top of the column.

In the text

	Fig. 7. Comparison of αz maps measured from the original HMI observation (top), and from the model without (middle) and with (bottom) fluctuations of the magnetic field. The color bar at the top shows the αz values in Mm−1.
In the text

	Fig. 8. Temporal evolution of αz in the leading sunspot of active region NOAA 11072. The time relative to the first observation (2010.05.25 03:00:00 TAI) is shown in hours in the top left of each panel. The black solid line represents the umbra–penumbra boundary.
In the text

	Fig. 9. Temporal evolution of the sign of αz in the leading sunspot of active region NOAA 11072. The time relative to the first observation (2010.05.25 03:00:00 TAI) is shown in hours in the top left of each panel. The black solid line represents the umbra–penumbra boundary.
In the text

Fig. 10. HMI continuum image of the leading sunspot of active region NOAA 11072 (2010.05.25 03:00:00 TAI), which was used as a reference for the model presented in this work. The black solid line outlines the umbral area that was considered for most twist calculation methods. The white solid line and the black dotted line correspond to the areas that were used to get spatial averages of αz close to the center of the spot (${\alpha }_{\mathrm{av}}^{\mathrm{center}}$) and over the whole spot (${\alpha }_{\mathrm{av}}^{\mathrm{spot}}$), respectively.

In the text

Fig. 11. Monte-Carlo simulation results for methods described in Sect. 5. The q values have been multiplied by two in order to make them directly comparable to αz (see Appendix F). The black line indicates the reference value of the model, which is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the Monte-Carlo mean.

In the text

	Fig. 12. Monte-Carlo simulation results when αz is spatially averaged over the whole sunspot. The solid black line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the mean.
In the text

	Fig. B.1. Normalized PSF for HMI estimated by Yeo et al. (2014).
In the text

	Fig. B.2. Cuts through the center along the x-axis (top panel) and the y-axis (bottom panel) of the average spatial correlation of neighboring pixels (Fig. 4) and the normalized estimated PSF by Yeo et al. (2014) (Fig. B.1). The errorbars represent the standard error of the spatially averaged correlation.
In the text

Fig. D.1. Monte-Carlo simulation results for methods based on αz when Stokes’ theorem is used for calculating Jz. The black solid line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The dark/light gray shaded areas represent the range of 1, 2, or 3σ around the Monte-Carlo mean.

In the text

Fig. D.2. Monte-Carlo simulation results when αz is calculated with Stokes’ theorem and spatially averaged over the whole sunspot. The black solid line indicates the reference value of the model that is expected from the model without any fluctuations. The red dashed line shows the mean value from the Monte-Carlo simulations. The gray shaded areas represent the range of 1, 2, and 3σ around the mean.

In the text