© ESO 2018

1. Introduction

The geometric height at which unity continuum optical depth is reached is depressed within sunspots relative to the quiet Sun. This so-called Wilson depression (Wilson & Maskelyne 1774) is caused by a lower opacity within the sunspot due to the lower temperature and a reduced gas pressure. The magnetic pressure and the curvature force of the strong magnetic field of sunspots balance this reduced gas pressure with the gas pressure of the surrounding quiet Sun. The connection between the Wilson depression and the strength and geometry of the magnetic field makes the Wilson depression an important quantity for understanding the structure of sunspots. In particular, it is not known by how much the curvature force of the magnetic field contributes to stabilizing the sunspot.

Unfortunately, the Wilson depression remains one of the more poorly known parameters of sunspots. Several studies have tried inferring the Wilson depression by making use of the horizontal force balance between the sunspot and the surrounding quiet Sun. However, since we do not know by how much the curvature force contributes to the force balance, an accurate estimate of the Wilson depression with this method is not possible. Depending on the assumed influence of the curvature force, the derived Wilson depression lies in the range between 400 and 1000 km (Solanki et al. 1993; Martínez Pillet & Vázquez 1993; Mathew et al. 2004).

The Wilson depression can also be estimated geometrically when the sunspot approaches the limb (i.e., via the Wilson effect). However, this method is influenced by radiative transfer or changes of the size of the umbra and penumbra with height (see, e.g., the discussion in Solanki 2003). Hence, the Wilson depression cannot be inferred very accurately when using this method. Gokhale & Zwaan (1972) derived an average Wilson depression z_W of 600 ± 200 km based on the Wilson effect. In contrast, Prokakis (1974) measured a significantly larger Wilson depression of 950 − 1250 km with large spots having a higher Wilson depression (z_W = 1500 − 2100 km) than small spots (z_W = 700 − 1000 km). However, later results (z_W = 500 − 1000 km) obtained by Balthasar & Wöhl (1983) are in better agreement with those of Gokhale & Zwaan (1972).

Here we present an alternative method for measuring the Wilson depression that does not need to make any assumptions on the geometry of the magnetic field in sunspots or on the structure of the sunspot. Our method is based on imposing the divergence-free condition on the magnetic field vector deduced from the inversion of observed Stokes profiles. This approach has already been used by Puschmann et al. (2010) to derive the small-scale corrugation of the τ = 1 layer of a small patch within the penumbra of a sunspot. Here, we have modified this method to provide the large-scale corrugation of the τ = 1 layer within the entire sunspot. We first performed a test of the method: we use it to derive the Wilson depression from synthetic Hinode observations generated from two MHD simulations of sunspots (Rempel 2012, 2015). After that we applied it to one spot observed with Hinode. We also compared our results with the Wilson depression derived using the pressure method.

2. Data

2.1. Synthetic Hinode data

We generated synthetic Hinode data starting from two MHD simulations of sunspots by Rempel (2012, 2015) computed with the MURaM code (Vögler et al. 2005). Simulation run 1 has a size of 49.152 × 49.152 × 2.048 Mm³ with a resolution of 12 × 12 × 8 km³ (4096 × 4096 × 256 pixels). Simulation run 2 has a size of 61.44 × 61.44 × 2.976 Mm³ with a resolution of 48 × 48 × 24 km³ (1280 × 1280 × 124 pixels). Figure 1 shows continuum intensity images and Table 1 lists some parameters of the spots (the area of the spots, the maximum and the average magnetic field in the umbra at logτ = −0.9, and the minimum and average temperature in the umbra at logτ = 0). We defined the inner and outer boundary of the penumbra as 30% and 90% of the continuum intensity level of the quiet Sun, respectively, after smoothing the continuum images with a 2D Gaussian with σ = 812 km.

Fig. 1. Maps of the continuum intensity of the sunspots analyzed in this study. Panels A and B: sunspot simulations by Rempel (2012, 2015), panel C: NOAA AR 10923 observed by Hinode SOT/SP on 14 November 2006. The blue and green contours indicate the outer and inner penumbral boundaries (defined as 30% and 90% of the continuum intensity level of the quiet Sun, after smoothing the continuum images with a 2D Gaussian with σ = 812 km). The red boxes outline the regions that we used for deriving the Wilson depression in the two MHD simulations. We computed the Wilson depression across the entire FOV for AR 10923.

We computed line profiles of the pair of Fe I lines at 6301.5 Å and 6302.5 Å for all Stokes parameters from these simulations using the SPINOR code (Frutiger et al. 2000); which uses the STOPRO routines for the forward calculation (Solanki 1987) for μ = 1. In order to save computation time, we only computed the line profiles for a narrow strip across the sunspot for simulation run 1 (indicated by the red lines in the left panel in Fig. 1). For simulation run 2, we used the full field-of-view (FOV). We then degraded the data to the spectral resolution of Hinode SOT/SP and rebined the data to the pixel size of Hinode SOT/SP (0.16″). We did not convolve the data with the optical point-spread function (PSF) since an inversion with the spatially coupled version of SPINOR (van Noort 2012; van Noort et al. 2013) would remove the influence of the PSF from the data anyway. We note that the sunspot AR 10923 observed by Hinode and analyzed here, was inverted with the spatially coupled version of SPINOR and we intend to apply this technique mainly to data inverted in the same manner. Finally, we added normally distributed noise to all four Stokes parameters with a signal-to-noise ratio that is proportional to the square root of the intensity and that reaches an average value of 1000 in the continuum.

2.2. Hinode observations of AR 10923

We used spectropolarimetric observations of AR 10923 made on 14 November 2006 in the Fe I line pair at 6301.5 Å and 6302.5 Å from the spectropolarimeter on the Solar Optical Telescope (SOT/SP, Kosugi et al. 2007; Lites 2007; Ichimoto et al. 2008; Lites et al. 2013) onboard the Hinode spacecraft. At that time, this active region was located close to disk center (x = 66.8″, y = −114.4″, μ = 0.99, see panel C in Fig. 1). The SOT/SP provides the full set of Stokes parameters along both spectral lines with a pixel sampling of ∼0.16″. This spot exhibits two filaments stretching far into the umbra, nearly splitting the umbra in two parts. This spot is significantly larger than the two simulated spots, however, the magnetic field in the umbra is weaker (see Table 1).

3. Deriving the atmospheric conditions

We derived the atmospheric parameters of both, the two simulated spots and AR 10923, by inverting the maps of the Stokes parameters with the SPINOR code under the assumption of local thermodynamic equilibrium (LTE). In case of the Hinode observations, we used the spatially coupled version of SPINOR (van Noort 2012; van Noort et al. 2013). The synthetic data were inverted using the non-coupled version of SPINOR since we did not take into account the PSF when generating the synthetic data. For all three spots, we set three nodes in optical depth, placed at logτ = −2.5, −0.9, 0. Figures 2–4 show the results of the inversion for the magnetic field strength B, the inclination γ, and the azimuth φ for all the spots at these optical depths. Using a spline interpolation, we then remapped the results of the inversion on an equidistant grid in logτ ranging from −6 to +1.5 with a sampling of 0.1.

Fig. 2. Inverted magnetic field vector along the analyzed stripe through simulation run 1. Top three rows: magnetic field intensity B, central three rows: inclination γ, and bottom three rows: azimuth of the magnetic field ϕ. We show each parameter at three different optical depths, from top to bottom: logτ = −2.5, −0.9, 0.

Fig. 3. Inverted magnetic field vector for simulation run 2. Top row: magnetic field intensity B, central row: inclination γ, and bottom row: azimuth of the magnetic field ϕ. From left to right: logτ = −2.5, −0.9, 0. The red box highlights the region, where we derived the Wilson depression.

Fig. 4. Same as Fig. 3 for the Hinode observation of AR 10923.

Afterward, we resolved the 180° azimuthal ambiguity by using the non-potential magnetic field computation method (NPFC, Georgoulis 2005). We could test the performance of the NPFC code with the help of the synthetic data. In most cases, the code resolved the ambiguity correctly within the sunspot. However, the code sometimes failed in regions where the magnetic field vector is almost exactly parallel to the line-of-sight. This occured predominantly at optical depths where the inversion is less accurate. We addressed this issue by applying the NPFC code only to the data at logτ = −0.9. We then successively resolved the ambiguity at greater and lower heights by demanding the magnetic field vector to vary smoothly with height. We selected the solution of the ambiguity, where the magnetic field vector is the closest to the one in the adjacent layer, where we have already resolved the ambiguity. This led to an accurate disambiguation within the sunspot. In the surroundings of the sunspots, however, the NPFC code failed in many cases to find the correct solution, even at logτ = −0.9.

4. Measuring the Wilson depression

4.1. Method

We determined the geometric height of the τ = 1 layer of the sunspot in a similar way as introduced by Puschmann et al. (2010), who could infer the geometrical height of the τ = 1 surface in an inverted atmosphere of a small patch of the penumbra of a sunspot observed with Hinode. Their method is based on the divergence-free condition of the magnetic field vector and on ensuring force balance. These two conditions are not fulfilled in an inverted atmosphere since the unknown height of the τ = 1 layer causes an offset in the geometric height scale of the individual pixels. However, both the divergence of the magnetic field and the deviations from force balance are minimized when shifting the atmosphere at each pixel by the corresponding height of the τ = 1 surface. Hence, this approach allows to infer the Wilson depression. Puschmann et al. (2010) use the following merit function:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{m,n}{w}_1{|\mathit{F}|}^2+w_2{(\mathbf{\nabla }·\mathit{B})}^2.}\end{array}$$(1)

The sum is over all the pixels in the images (indicated by the indices m and n). The first term consists of the residual force F = J × B + ρg − ∇P_g. Here, J is the current density, B is the magnetic field vector, ρ is the density, g is the surface gravity of the Sun, and P_g is the gas pressure. The influence of flows on the force balance is neglected. The second term in the merit function is the divergence of the magnetic field. The relative contributions of the two terms are balanced by the coefficients w₁ and w₂.

The magnetic field vector and all other atmospheric parameters in the merit function depend on the geometric height z. This height scale has two contributions. The first one is the z-scale z_rel of the stratification relative to the τ = 1 layer of each pixel, which is provided by the inversion. The second one is the height of the τ = 1 layer at each pixel (i.e., the Wilson depression z_W), which we want to determine:

$$\begin{array}{c}\hfill z(x,y,\tau )=z_{\mathrm{rel}}(x,y,\tau )+z_{\mathrm{W}}(x,y).\end{array}$$(2)

Changing z_W at the individual pixels corresponds to shifting the inverted atmosphere vertically. The merit function has its minimum when the inverted atmospheres for all pixels are aligned with respect to each other, that is, when z_W(x, y) is the corrugation of the τ = 1 surface. The alignment of the individual atmospheres can be evaluated at different heights z_rel relative to the τ = 1 surface. Puschmann et al. (2010) used z_rel = 200 km and set z_W = 0 as the height used to compute the merit function for an arbitrary pixel. They used a genetic algorithm for minimizing the merit function.

Here, we wanted to extend the method of Puschmann et al. (2010) to derive the corrugations of the τ = 1 layer of the entire spot. This is not straightforward since the Wilson depression at each pixel is a free parameter in the merit function. So, the number of free parameters becomes extremely large when applying this method to the entire spot. We reduced the number of free parameters by restricting the analysis to the large-scale corrugation of the τ = 1 layer. An efficient way to achieve this is to write the Wilson depression z_W(x, y) as a Fourier series of the horizontal wavenumbers k _x and k_y and then drop the terms above a maximum wavenumber in the Fourier series of z_W, therefore removing the small-scale corrugations. Now, the number of free parameters does not correspond to the number of pixels anymore, but it depends on the maximum wavenumber that is considered. This approach allowed us to reduce the number of free parameters in the merit function significantly without affecting the determination of the large-scale Wilson depression.

However, neglecting small-scale corrugations of the τ = 1 layer affects the divergence of the magnetic field, especially on small spatial scales. We addressed this problem by writing the merit function in Fourier space using Parseval’s theorem:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{m,n}{(\mathbf{\nabla }·\mathit{B})}^2=\frac{1}{N_x{N}_y}\sum _{k,l}{|\mathcal{F}(\mathbf{\nabla }·\mathit{B})|}^2.}\end{array}$$(3)

Here, the symbol ℱ indicates a 2D discrete Fourier transform (in the x- and in the y-direction), defined as

$$\begin{array}{c}\hfill {\widehat{a}}_{k,l}={\displaystyle \sum _{m=0}^{N_x-1}\sum _{n=0}^{N_y-1}{\mathrm{e}}^{-2\pi i(\frac{\mathit{mk}}{N_x}+\frac{\mathit{nl}}{N_y})}·a_{m,n}.}\end{array}$$(4)

The parameters N_x and N_y are the number of pixels along the x-axis and the y-axis, respectively. The indices k and l indicate the dimensionless wavenumber indices and are connected to the physical wavenumber by the relation $k_x=\frac{2\pi }{N_x\mathrm{\Delta }x}·k$ in case of the x-direction with Δx being the spatial resolution. Contrary to Puschmann et al. (2010), we did not include the force terms in our merit function. In our case, including the force terms leads to a less accurate determination of the Wilson depression. This is caused by neglecting advection (ρ (v · ∇)v) in the force balance (spectropolarimetry can only measure the line-of-sight velocity, not the full velocity vector). In the penumbra of the two simulated spots, the advection term reaches a similar magnitude as the Lorentz force due to the Evershed and the reverse Evershed flow. Hence, ignoring advection causes the force term in the merit function not to be sensitive to the large-scale Wilson depression anymore.

All terms in the merit function are positive, meaning that the divergence has to be zero at all spatial scales in order to minimize the merit function. So, one can restrict the analysis to the large spatial scales of the divergence only, which are less affected by neglecting the small-scale corrugations of the τ = 1 layer. Both for the Wilson depression and for the divergence of the magnetic field vector, we only considered spatial scales up to maximum dimensionless wavenumbers k_max and l_max. This decreased the effective spatial resolution significantly (given by the wavelength corresponding to the indices k_max and l_max). When using a maximum dimensionless wavenumber k_max = l_max = 3, for example, we retrieved an effective spatial resolution of 14 Mm for the part of simulation run 2 that is highlighted by the red box in Fig. 1.

We filtered both, the Wilson depression and the magnetic field vector in Fourier space to include only the information about large spatial scales. Hence, a map of the large-scale magnetic field at a fixed geometric height is affected by the filtering in two ways. Ignoring the small-scale corrugations of the τ = 1 layer causes an error in the derived magnetic field vector if the magnetic field varies strongly with height at a given pixel. Afterward, variations of the magnetic field itself on small spatial scales were removed. This filtering should remove most of the errors introduced by neglecting the small-scale corrugations of the τ = 1 layer.

This assumption can be tested using the simulated sunspots. We computed maps of all three components of the magnetic field vector at a fixed geometrical height (200 km above the mean height of the τ = 1 layer) from simulation run 2 for two different cases. In the first case, we derived the magnetic field on a grid in geometric height that includes small-scale corrugations of the τ = 1 surface, in the second case, we neglected corrugations above a maximum dimensionless wavenumber k_max = l_max = 3. Afterward, we also removed the signal at spatial scales above k_max or l_max from the maps of the magnetic field vector. The correlation coefficient between these two sets of maps of the magnetic field vector is 0.98 for all three components of the magnetic field vector, indicating that the small-scale corrugations of the τ = 1 surface do not have a strong influence on the large-scale magnetic field. Using this approach leads to the following expression for the merit function:

$$\begin{array}{cc}\hfill {\chi }^2& =\frac{1}{N_x{N}_y}{\displaystyle \sum _k^{k_{\max }}\sum _l^{l_{\max }}{|\mathcal{F}(\mathbf{\nabla }·\mathit{B})|}^2}\hfill \end{array}$$(5)

$$\begin{array}{cc}& =\frac{1}{N_x{N}_y}{\displaystyle \sum _k^{k_{\max }}\sum _l^{l_{\max }}{|i{k}_x{\widehat{B}}_x+i{k}_y{\widehat{B}}_y+\mathcal{F}\left(\frac{\partial B_z}{\partial z}\right)|}^2·}\hfill \end{array}$$(6)

Here, ${\widehat{B}}_x(k_x,k_y,z)$ and ${\widehat{B}}_y(k_x,k_y,z)$ are the Fourier transforms of B_x(x, y, z) and B_y(x, y, z), respectively. Again, the geometric height z consists of the relative height above the τ = 1 surface z_rel and the Wilson depression z_W. We evaluated the merit function around a fixed reference height z_rel and required the derived Wilson depression to have a mean value of zero. Like Puschmann et al. (2010), we minimized our merit function using a genetic algorithm. We ran the genetic algorithm ten times and used the solution which minimizes the merit function the most. After running the genetic algorithm, we defined z_W = 0 as the height of the τ = 1 surface averaged over the part of the FOV that has a distance of at least 7 Mm to the sunspot. This minimum distance is necessary because the height of the τ = 1 surface can be depressed in the close surroundings of the sunspot.

Our method relates the magnetic field vector in the quiet Sun and in the umbra at the same geometric height. This requires an accurate measurement of the magnetic field vector over a broad range in height since sunspots have a Wilson depression of a few 100 km. So, in order to derive reliable estimates of the Wilson depression, the magnetic field vector retrieved by the inversion needs to be reliable over a broad range in height around the reference height z_rel. Errors in the inversion occur predominantly at high optical depths (we did not consider data at heights z_rel < 85 km) and at very low optical depths. We chose z_rel to be as small as possible while ensuring that all inferred heights are greater than our threshold of 85 km. The inversion also suffers from systematic errors, such as the assumption of hydrostatic equilibrium. In the simulation data, this is a good approximation in granules and in the umbra but not so good in the intergranular lanes and in the penumbra. Errors occurring on small spatial scales should have limited influence though, due to our focus on large spatial scales.

We performed our analysis in Fourier space. This means that the magnetic field vector and the derived Wilson depression have periodic boundary conditions when expressing them in Fourier space. This is a good approximation when the field of view also contains the surroundings of the sunspot, where the magnetic field and the Wilson depression are significantly lower than within the spot.

4.2. Synthetic Hinode data

Figures 5 and 6 show the Wilson depressions derived from synthetic Hinode data for the two simulated sunspots. For comparison, we also show the true values of the Wilson depression for both spots which were inferred by resampling the MHD cubes on a grid in optical depth using SPINOR. We set z_W = 0 as the height of the τ = 1 surface averaged over the region of the FOV that has a distance of at least 7 Mm to the outer boundary of the penumbra (see dashed black contours in Figs. 5 and 6). In Table 2, we quantitatively compare the Wilson depression inferred by our method with the true solution and results from the pressure method (see Sect. 4.3). We list the maximum Wilson depression (z_{W, max}) and its average over the umbra (z_{W, U}) and over the penumbra (z_{W, PU}), both for the full spatial resolution (referred to as “full res.” in the table) and after degrading it to the resolution of our method (“filtered” in the table). We also derived error estimates of the Wilson depression for the umbra (Δz_{W, U}), the penumbra (z_{W, PU}), and the area outside the sunspot (Δz_{W, QS}). We derived these errors from maps of the absolute value of the difference between our solutions and the true Wilson depression (after degrading it to the spatial resolution of our derived Wilson depressions). We then defined the errors to be the average of this difference map over the respective region of the spot.

Fig. 5. Wilson depression of simulation run 1. Panel A: true Wilson depression, panel B: Wilson depression derived from the inversion. The black and white contours indicate the outer and the inner penumbral boundary. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Panel C: Wilson depression along a cut across the sunspot (marked by the dashed lines in panels A and B). Red: Wilson depression derived from synthetic observations, blue: true Wilson depression, black: true Wilson depression filtered in Fourier space to the same range in wavenumber as the derived solution from the inversion, for direct comparison with the red curve.

Fig. 6. Wilson depression of simulation run 2. Panel A: true Wilson depression, panel B: Wilson depression derived from the inversion. The black and white contours indicate the outer and the inner penumbral boundary. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom row: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top row). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Blue: true Wilson depression, black: true Wilson depression filtered in Fourier space to the same range in wavenumber as the derived solution from the inversion (shown in red).

For spot 1, we computed the Wilson depression only along a slice across the sunspot to limit the computation time due to the large number of pixels. Within this slice, the magnetic field vector is not periodic in the y-direction, in particular the B_y component. Since this non-periodic boundary condition would affect the filtering in Fourier space, we computed the divergence using the full spatial resolution in the y-direction. However, we did apply Fourier filtering in the x-direction. We expressed the Wilson depression in Fourier space, though, in order to reduce the number of free parameters (k_max = 3, l_max = 1). The main differences between our solution and the true Wilson depression occur on small spatial scales, which we do not consider in our method. We degraded the true solution to the same spatial resolution as our derived Wilson depression (black curve in the bottom panel of Fig. 5) in order to remove the influence of the differences in spatial resolution. Our derived Wilson depression is generally in good agreement with the degraded true solution, but the maximum Wilson depression is ∼100 km larger than the degraded true solution (see Table 2). The averages of the Wilson depression over the umbra and the penumbra are in good agreement, though. In both cases, the mean absolute difference between our solution and the true Wilson depression is on the order of 60 km.

This error estimate is affected both by uncertainties introduced by the inversion of the Stokes parameters and by the stability of the genetic algorithm. We evaluated the stability of the minimization by computing the standard deviation of the Wilson depression inferred by running the genetic algorithm ten times. The resulting errors (34 km for the maximum Wilson depression and 17 km for the average over the umbra) are significantly smaller than the difference between the derived and the true Wilson depression. This suggests that the error is dominated by uncertainties in the inversion.

Spot 2 exhibits a very complex Wilson depression in the umbra with changes of several 100 km on small spatial scales (see Fig. 6). These small-scale corrugations cannot be resolved by our method (here we use k_max = l_max = 3). Again, for better comparison, we smoothed the true solution in order to account for this effect. Afterward, both the maximum Wilson depression and its averages over the umbra and the penumbra were in good agreement (see Table 2). A remarkable difference between our solution and the true solution occurs at x ≈ 24 Mm, y ≈ 20 Mm. At this position, there is a filament stretching to the umbra and the τ = 1 layer is located a few 100 km higher than in the surrounding umbra. Our method finds an extremely large Wilson depression at this location, in contradiction to the true solution. This error in the derived Wilson depression is probably caused by uncertainties in the inversion. Along the filament, the magnetic field changes strongly with height (by up to 2000 gauss over a range of 100 km at around τ = 1). This height dependency cannot be described correctly by the inversion. This feature not only affects the Wilson depression locally but also in the surrounding umbra. This is because the spatial resolution of our method is very low and because the Fourier transform is non-local (a perturbation of the signal at a given point also affects distant points after filtering in Fourier space). Again, we computed the mean absolute value of the difference between our derived Wilson depression and the true solution. The differences are roughly a factor of 1.5 larger than the ones for spot 1 and much larger than the uncertainty introduced by the genetic algorithm (13 km for the maximum Wilson depression and 9 km for the average over the umbra), in agreement with the findings from spot 1.

The inferred Wilson depression depends on the parameters of the merit function. Apart from the Wilson depression, the merit function (Eq. (6)) depends on the range of spatial scales that is considered (given by k_max and l_max) and the height above τ = 1 (z_rel) where the merit function is evaluated. As stated in the previous section, the height z_rel is a crucial parameter since it determines, what range in height (and optical depth) is used for inferring the Wilson depression. We chose the reference height z_rel to be as small as possible while ensuring that all inferred heights are greater than our threshold of 85 km (350 km for spot 1, 300 km for spot 2). For larger reference heights, the derived Wilson depression becomes less reliable. We also restrained the solution to low wavenumbers. Larger wavenumbers should be avoided since in that case the solution for the Wilson depression is more likely to exhibit an oscillatory behavior, especially in the quiet Sun.

4.3. Comparison with the pressure method

For comparison, we also derived the Wilson depression by demanding horizontal force balance between the sunspot and the surrounding quiet Sun (Martínez Pillet & Vázquez 1993; Solanki et al. 1993; Mathew et al. 2004). The reduced gas pressure inside the umbra is compensated by the Lorentz force of the strong magnetic field. After a few assumptions, such as radial symmetry of the magnetic field vector, no magnetic field in the azimuthal direction and a negligible Evershed flow (Maltby 1977), this can be expressed as

$$\begin{array}{c}\hfill P_{\mathrm{gas}}(r=a,z)-P_{\mathrm{gas}}(r,z)=B_z^2(r,z)/8\pi +F_{\mathrm{c}}(r,z)/8\pi .\end{array}$$(7)

Here, r = 0 refers to the center of the umbra and r = a to a point in the quiet Sun. The second term on the right-hand side, F_c is the curvature integral:

$$\begin{array}{c}\hfill F_{\mathrm{c}}=2{\int }_r^a{B}_z(r^{\prime },z)\frac{\partial B_{\mathrm{r}}(r^{\prime },z)}{\partial z}\mathrm{d}{r}^{\prime }.\end{array}$$(8)

Equation (7) can be used to infer the Wilson depression of a sunspot by comparing the measured total pressure (gas pressure plus magnetic pressure) at some optical depth with the pressure stratification of a reference model atmosphere of the quiet Sun (which extends at least to the depth of the Wilson depression). The derived Wilson depression depends on the curvature integral, which cannot directly be inferred from the observations. Commonly, one assumes F_c = 0 (see e.g., Mathew et al. 2004).

We applied this method to both simulated sunspots. We used a horizontal average of a quiet-sun region from simulation run 2 as the reference atmosphere and extract the pressure and the magnetic field from the inversions at τ = 1. The resulting Wilson depressions are in the range 420 − 490 km, when averaging over the umbra (see Table 2), in good agreement with the results from previous studies for F_c = 0, which lay in the range 400–450 km (Martínez Pillet & Vázquez 1993; Solanki et al. 1993; Mathew et al. 2004). Figure 7 compares the Wilson depression inferred from this method for simulation run 2 with our results based on the divergence of the magnetic field and the true solution.

Fig. 7. Comparison of the Wilson depression derived from the divergence method and the pressure method for synthetic Hinode data generated from the simulation run 2. Panel A: Wilson depression derived using the pressure method and panel B: Wilson depression derived using the divergence of the magnetic field (the same image as shown in panel B in Fig. 6). The black and white contours indicate the outer and the inner penumbral boundary, respectively. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom panels: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top panels). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Blue: true Wilson depression (zW), red: Wilson depression derived using the divergence method (zW, d), black: Wilson depression inferred from the pressure method (zW, p), green: Wilson depression inferred from the pressure method degraded to the spatial resolution of the divergence method (zW, p, filt).

The pressure method results in a Wilson depression which is about 110 km smaller than the true value for the simulations. This suggests that the curvature integral is positive in both spots. We obtained the curvature integral by assuming that the difference between the true Wilson depression and the one obtained from the pressure method is due to the curvature integral neglected by the pressure method. For the spot in simulation 2, an error in the Wilson depression of 110 km corresponds to F_c/8π ≈ 2.8 × 10⁵ dyn cm⁻² in the umbra. This is about half of the value of the magnetic pressure in the umbra (4.9 × 10⁵ dyn cm⁻²).

The inferred value of F_c/8π is in good agreement with the curvature integral derived directly from the MHD cube. The curvature integral depends on height, between z = −700 km and z = −400 km, F_c/8π varies between 1.5 × 10⁵ dyn cm⁻² and 5.9 × 10⁵ dyn cm⁻². At greater heights, it becomes negative, reaching −1.2 × 10⁶ dyn cm⁻² at z = −240 km. This change of sign is caused by the vertical derivative of the radial component of the magnetic field in Eq. (8). The strength of B_r reaches its maximum at around τ = 1 and decreases toward larger or smaller heights, causing a change of sign of its vertical derivative.

In the penumbra, the pressure method significantly underestimates the Wilson depression (by more than 100 km for the spot in simulation run 2). This is caused by using only the vertical component of the magnetic field when estimating the magnetic pressure. In the penumbra, the magnetic field is predominantly horizontal, meaning that the magnetic pressure and, hence, also the Wilson depression is underestimated when using B_z only.

4.4. Real Hinode data: case study of AR 10923

We also tested the performance of our method on the Hinode observations of AR 10923 (see Fig. 8 and Table 2). Again, we used a maximum dimensionless wavenumber k_max = l_max = 3, which corresponds to an effective spatial resolution of ∼19 Mm in the x-direction and ∼21 Mm in the y-direction. The resulting Wilson depression looks similar to the one of the simulated spots, although the average Wilson depression is by ∼40 km deeper. It displays two distinct local maxima, where the Wilson depression exceeds 700 km. These maxima occur in regions with strong magnetic field (see top row of Fig. 4). These two regions are separated by a filament, along which the Wilson depression is reduced by ∼60 km. The Wilson depression smoothly decreases from the center of the umbra toward the outer penumbral boundary. The height of the τ = 1 surface is more or less constant outside the sunspot, which is to be expected at the large spatial scales considered here.

Fig. 8. Comparison of the Wilson depression derived from the divergence method and the pressure method for AR 10923. Panel A: Wilson depression derived using the pressure method and panel B: Wilson depression derived using the divergence method. The black and white contours indicate the outer and the inner penumbral boundary, respectively. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom panels: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top panels). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Black: Wilson depression inferred from the pressure method, red: Wilson depression inferred from the pressure method degraded to the spatial resolution of the divergence method, green: Wilson depression derived using the divergence method.

Unfortunately, we could not derive a reliable estimate of the error of the Wilson depression for this spot. As found from the simulated spots (Sect. 4.2) the main errors in z_W arise from uncertainties in the inversions. Obtaining the errors in inverted values is never straightforward, also from the spatially coupled inversion scheme of van Noort (2012) and van Noort et al. (2013), which was used to invert this spot. We expect that the uncertainty is similar to or larger than for the synthetic data (i. e., roughly 95 km or larger).

The Wilson depression derived using the pressure method (see Fig. 8 and Table 2) is significantly smaller (by about 180 km) than the one based on our method, suggesting that the curvature integral is positive in this spot, as well. Assuming the z_W deduced via the divergence method to be correct, we retrieved F_c/8π ≈ 4.21 × 10⁵ dyn cm⁻². This is higher than the pressure from the vertical component of the magnetic field (3.5 × 10⁵ dyn cm⁻²) in the umbra at τ = 1. The high value of F_c of AR 10923 is consistent with typical values found in other spots, while the results for the simulated spots are lower than what has been reported before. Solanki et al. (1993) estimate the curvature integral to be within the range 3.5 × 10⁵ dyn cm⁻² ≲ F_c/8π ≲ 1.6 × 10⁶ dyn cm⁻² for the sunspot that they studied. The value we found lies within this range.

5. Discussion

We could successfully reproduce the Wilson depression of two simulated sunspots by minimizing the divergence of the magnetic field. We also applied this method to a sunspot observed with Hinode (AR 10923) and derived a Wilson depression that is consistent with the results statistically provided by some studies of the Wilson effect (∼600 km, Gokhale & Zwaan 1972). Our results suggest that the pressure due to the vertical component of the magnetic field and the curvature integral contribute by a similar amount to the horizontal force balance within sunspots, as has already been suggested by previous investigations (Martínez Pillet & Vázquez 1993; Solanki et al. 1993; Mathew et al. 2004). The Wilson depression of the sunspot simulations by Rempel (2012, 2015) is smaller by ∼50 km than the one of AR 10923. In the simulated spots, the curvature integral is significantly lower than the pressure due to the vertical component of the magnetic field, leading to a lower Wilson depression. This might be caused by the different geometry of the simulated spots compared to AR 10923 (see Table 1). The simulated spots exhibit a slightly stronger magnetic field but they are significantly smaller than the AR 10923 spot.

The main limitation of our method is that it requires a reliable inversion of the full magnetic field vector over a broad range in height. In the umbra, our method requires an estimate of the magnetic field up to ∼800 km above the τ = 1 surface. The lines used in this study are not very sensitive at these height and so, the magnetic field is predominantly extrapolated from lower layers. Observations of multiple spectral lines with a broad range of formation heights are needed in order to retrieve a better inversion (as is planned, e.g., with the upcoming Sunrise III mission). Systematic errors also arise from the assumption of LTE and hydrostatic equilibrium in the inversion, although these are likely small compared with the uncertainty resulting from the inversions. Also, modeling the height dependency of atmospheric parameters with splines is not always a good representation of the true atmosphere, especially when there are strong gradients with height. In addition, the 180°-ambiguity of the Zeemann-effect needs to be resolved accurately across the entire spot at all optical depths. Fortunately, some of these inaccuracies of the inversion occur on small spatial scales, so that their influence on large spatial scales should be limited. For example, as shown in Sect. 4.1, neglecting the small-scale corrugations of the τ = 1 layer has only a small influence on the large-scale divergence of the magnetic field. A detailed determination of the errors in the spatially coupled inversion is beyond the scope of this paper. We therefore use the tests made with the synthetic data as a rough guide to the total error.

The Wilson depression derived from the synthetic Hinode data exhibits an error on the order of ∼95 km. In case of real observations, the error is probably somewhat higher. The synthetic data were generated using the same radiative transfer code and the same line parameters as were used for the inversion. Hence, the inversion of the synthetic data is likely to be more accurate than for real observations. Bearing these arguments in mind, we assign a preliminary error of 100 km to the Wilson depression derived from the Hinode observations.

The estimated error of our method is mainly statistical in nature. Any systematic component is significantly lower than the one of the Wilson effect or of the pressure method. As explained in the introduction, the values of the Wilson depression derived using the Wilson effect vary by more than 1000 km between different studies in a systematic manner. In case of the pressure method, the inferred Wilson depression depends on the assumed value of the curvature integral.

Our method is limited to the large-scale corrugations of the τ = 1 surface. Measuring the Wilson depression with a higher spatial resolution requires including more Fourier coefficients in the Fourier series of the Wilson depression, which affects the minimization of the merit function. In addition, the divergence of the magnetic field is dominated by the largest spatial scales. When evaluating the divergence of the magnetic field for AR 10923 on a corrugated grid in geometrical height (300 km above τ = 1 at each pixel), 76% of the total variance of the divergence of the magnetic field occur within the range of spatial scales that we considered in Sect. 4.4 (k_max = l_max = 3). Hence, the merit function is not very sensitive to smaller spatial scales, at least when considering entire sunspots.

Acknowledgments

This work benefited from the Hinode sunspot database at MPS, created by Gautam Narayan. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 695075) and has been supported by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea.

References

All Tables

All Figures

Fig. 1. Maps of the continuum intensity of the sunspots analyzed in this study. Panels A and B: sunspot simulations by Rempel (2012, 2015), panel C: NOAA AR 10923 observed by Hinode SOT/SP on 14 November 2006. The blue and green contours indicate the outer and inner penumbral boundaries (defined as 30% and 90% of the continuum intensity level of the quiet Sun, after smoothing the continuum images with a 2D Gaussian with σ = 812 km). The red boxes outline the regions that we used for deriving the Wilson depression in the two MHD simulations. We computed the Wilson depression across the entire FOV for AR 10923.

In the text

	Fig. 2. Inverted magnetic field vector along the analyzed stripe through simulation run 1. Top three rows: magnetic field intensity B, central three rows: inclination γ, and bottom three rows: azimuth of the magnetic field ϕ. We show each parameter at three different optical depths, from top to bottom: logτ = −2.5, −0.9, 0.
In the text

	Fig. 3. Inverted magnetic field vector for simulation run 2. Top row: magnetic field intensity B, central row: inclination γ, and bottom row: azimuth of the magnetic field ϕ. From left to right: logτ = −2.5, −0.9, 0. The red box highlights the region, where we derived the Wilson depression.
In the text

	Fig. 4. Same as Fig. 3 for the Hinode observation of AR 10923.
In the text

Fig. 5. Wilson depression of simulation run 1. Panel A: true Wilson depression, panel B: Wilson depression derived from the inversion. The black and white contours indicate the outer and the inner penumbral boundary. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Panel C: Wilson depression along a cut across the sunspot (marked by the dashed lines in panels A and B). Red: Wilson depression derived from synthetic observations, blue: true Wilson depression, black: true Wilson depression filtered in Fourier space to the same range in wavenumber as the derived solution from the inversion, for direct comparison with the red curve.

In the text

Fig. 6. Wilson depression of simulation run 2. Panel A: true Wilson depression, panel B: Wilson depression derived from the inversion. The black and white contours indicate the outer and the inner penumbral boundary. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom row: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top row). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Blue: true Wilson depression, black: true Wilson depression filtered in Fourier space to the same range in wavenumber as the derived solution from the inversion (shown in red).

In the text

Fig. 7. Comparison of the Wilson depression derived from the divergence method and the pressure method for synthetic Hinode data generated from the simulation run 2. Panel A: Wilson depression derived using the pressure method and panel B: Wilson depression derived using the divergence of the magnetic field (the same image as shown in panel B in Fig. 6). The black and white contours indicate the outer and the inner penumbral boundary, respectively. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom panels: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top panels). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Blue: true Wilson depression (zW), red: Wilson depression derived using the divergence method (zW, d), black: Wilson depression inferred from the pressure method (zW, p), green: Wilson depression inferred from the pressure method degraded to the spatial resolution of the divergence method (zW, p, filt).

In the text

Fig. 8. Comparison of the Wilson depression derived from the divergence method and the pressure method for AR 10923. Panel A: Wilson depression derived using the pressure method and panel B: Wilson depression derived using the divergence method. The black and white contours indicate the outer and the inner penumbral boundary, respectively. The dashed black contour shows the region that was used for defining zW = 0 (distance of at least 7 Mm to the sunspot). Bottom panels: Wilson depression along cuts across the sunspot (marked by the dashed lines in the top panels). Panel C: a cut along the x-direction, panel D: a cut along the y-direction. Black: Wilson depression inferred from the pressure method, red: Wilson depression inferred from the pressure method degraded to the spatial resolution of the divergence method, green: Wilson depression derived using the divergence method.

In the text