Helioseismology of sunspots: how sensitive are travel times to the Wilson depression and to the subsurface magnetic field?
Key Words.:Sun: helioseismology – Sun: sunspots – Sun: interior – Sun: oscillations
In order to assess the ability of helioseismology to probe the subsurface structure and magnetic field of sunspots, we need to determine how helioseismic travel times depend on perturbations to sunspot models. Here we numerically simulate the propagation of f, p, and p wave packets through magnetic sunspot models. Among the models we considered, a 50 km change in the height of the Wilson depression and a change in the subsurface magnetic field geometry can both be detected above the observational noise level. We also find that the travel-time shifts due to changes in a sunspot model must be modeled by computing the effects of changing the reference sunspot model, and not by computing the effects of changing the subsurface structure in the quiet-Sun model. For p modes the latter is wrong by a factor of four. In conclusion, numerical modeling of MHD wave propagation is an essential tool for the interpretation of the effects of sunspots on seismic waveforms.
1 A numerical approach
Elucidating the subsurface processes of sunspot formation and evolution is important for understanding solar magnetism (Moradi et al. 2010; Rempel & Schlichenmaier 2011). Helioseismology is the only available tool for imaging the subsurface structure, magnetic field and dynamics of sunspots (Gizon et al. 2010), all of which are largely unknown.
Early approaches interpreted the seismic signature of sunspots in terms of small perturbations to a quiet-Sun reference model (Fan et al. 1995; Kosovichev 1996). However, the effects of the surface magnetic field and associated structural changes cause large perturbations in the wave speed (e.g. Lindsey & Braun 2005; Gizon et al. 2006). Recently, direct numerical simulations of MHD wave propagation through sunspot models have been used to capture the large effects of the magnetic field on the waves (Gizon et al. 2010; Cameron et al. 2011).
A study by Liang et al. (2013) measured the travel-time shifts caused by the sunspot in AR 9787. In an effort to understand and interpret these measurements we used numerical simulations. We simulated the interaction of surface-gravity (f) and acoustic (p and p) wave packets with several sunspot models. With spectroscopy the height of the Wilson depression in sunspots is not known to better than 50 km (Solanki et al. 1993). Here we investigate if seismic diagnostics can provide a stronger constraint. Further, we study the seismic consequences of two subsurface magnetic field geometries below sunspots: a monolithic flux tube versus a weaker, spread-out magnetic field. Finally, we explore the sensitivity of waves to sound-speed perturbations at different depths along the axis of a sunspot model.
These parametric studies enable us to test the assumptions of existing approximation methods used in the interpretation of sunspot seismology. By comparing with measured noise in seismic travel times we are able to draw conclusions about the detectability of different aspects of sunspot structure.
2 Sunspot models
2.1 Background quiet-Sun model (QS)
Our sunspot models are embedded in the background solar model CSM_A described by Schunker et al. (2011). The CSM_A model is based on Model S from Christensen-Dalsgaard et al. (1996), but was modified to be convectively stable so that linear waves can be propagated using the SLiM wave simulation code (Cameron et al. 2008). An indication of the sensitivity of the modes of oscillation in CSM_A as a function of height, , is given by their kinetic energy density (see Fig. 1). Throughout this paper we define as height measured upward from the bottom of the photosphere.
2.2 Reference sunspot model (RS)
Cameron et al. (2011) constructed a parametric sunspot model that causes a seismic response consistent with that observed for the sunspot in active region NOAA 9787 (Liang et al. 2013, , Paper I). The vector magnetic field of the sunspot model is cylindrically symmetric, i.e. it depends only on the horizontal distance from the axis of the sunspot, , and on height, . The model of Cameron et al. (2011) is a self-similar and monolithic magnetic field model specified by the strength of the field on the axis of the sunspot () as a function of height. Here we adopt the same sunspot model with the parameters given in Table 1 of Cameron et al. (2011), except that the parameters controlling the inclination of the magnetic field are modified slightly ( Mm and Mm) so that the magnetic field inclination at the umbra-penumbra boundary () better matches observed field inclinations in sunspots (Schunker et al. 2008). The solid black curves in Fig. 1 show a sample of the magnetic field lines. Throughout this paper we refer to this model as the reference sunspot model abbreviated as ‘RS’.
2.3 Changes to the Wilson depression
The Wilson depression is a physical depression of the surface (optical depth unity) caused by strong magnetic fields in the umbra of a sunspot (Loughhead & Bray 1958). The Wilson depression of our reference sunspot model is at km in the center of the umbra and at km in the middle of the penumbra as in Cameron et al. (2011).
In order to assess the ability of helioseismology to measure the height of the Wilson depression, we consider two additional sunspot models in which the Wilson depression is shifted by km with respect to the reference sunspot model. In these modified models, we translate the density, pressure, and sound-speed profiles in height by km, while keeping the magnetic field unchanged (models are not in magnetohydrostatic equilibrium).
2.4 Alternative magnetic field model
The monolithic reference sunspot model described in Sect. 2.2 is only one possibility among many. We considered a sunspot model with an alternative magnetic field below Mm, which fans out as a function of depth. This alternative model is motivated by, e.g., the dynamical disconnection model of Schüssler & Rempel (2005) whereby the flux tube is disrupted by convective motions at depths greater than 2 Mm. In the alternative magnetic sunspot model, we keep the sound speed, pressure and density profiles unchanged from the reference sunspot model.
More precisely, the alternative magnetic field is modified so that the magnetic field becomes weaker below Mm (above, it is identical to the reference sunspot model), smoothly transitioning to weak field strengths below Mm. Along the axis of the sunspot, the -component of the magnetic field, , is obtained from the reference field, , as follows:
where is the fifth order polynomial in such that is continuous and has continuous first- and second-order derivatives, i.e. where . The same field lines as for the reference sunspot model are shown for the alternative magnetic field model in Fig. 1.
2.5 Subsurface sound-speed perturbations
We also consider a set of sunspot models in which we modify the subsurface sound speed. In these models, the sound speed is given by
where is the sound speed in the reference sunspot model and
is a 3D Gaussian perturbation to the sunspot model. In the above expression, the relative amplitude of the perturbation is , the height of the centre of the perturbation is , and the width of the perturbations are fixed to Mm, Mm. For one experiment (Sect. 5.1) we use and and Mm. For a second experiment (Sect. 5.2) we use to and Mm.
3 Numerical simulations of wave propagation
We used the SLiM code (Cameron et al. 2007, 2011) to simulate the propagation of wave packets through the sunspot models described above. We use a simulation domain that covers Mm in each of the horizontal directions and extends from Mm below to Mm above the bottom of the photosphere. The depth of the simulation domain is sufficient to study acoustic modes with radial orders . The spatial resolutions are Mm per pixel horizontally and Mm per pixel vertically. The time resolution is s. We did a resolution study which showed that for the waves studied here the numerical solutions are converged.
We carry out simulations for individual wave packets with modes of fixed radial order, i.e. with either f, p, or p modes. The initial conditions () for the simulations are specified by the displacement, , and the velocity,. This requires specifying the wave amplitude distribution, , as a function of wavenumber, (see Cameron et al. 2008, 2011). For each wave packet is chosen in such a way that the vertical velocity in the simulation () is directly comparable to an observed cross-covariance function. This comparison is justified because the cross-covariance is closely related to the Green’s function in the far field (e.g., Gizon et al. 2010).
The wave packets in our simulations are comparable to the Solar and Heliospheric Observatory’s Michelson Doppler Imager (SOHO/MDI) observed cross-covariance functions presented by Liang et al. (2013). The peak of wave power occurs at and mHz for the f modes, at and mHz for the p modes, and at and mHz for the p modes. At , each wave packet peaks at Mm (at the left edge of Fig. 2), from where they propagate in the direction (to the right) and interact with the sunspot centred at Mm.
4 Travel times
4.1 Measuring travel times
Travel-time measurements are very useful to characterize the sensitivity of seismic waveforms to subsurface perturbations. In this paper we explore how changes to the reference sunspot model affect travel times.
We define the travel-time shift, , as the time by which a sliding reference wavelet, , must be shifted to best match the seismic waveform . Following Gizon & Birch (2002) and Liang et al. (2013) is the time, , that maximises the function
In practice, the reference wavelet was chosen to be from either the quiet-Sun simulation or the reference sunspot simulation. To denote these two kinds of measurements (which give similar results, see Sect 5.2), we use the notations and respectively. Furthermore, we use the notation to indicate the background model through which the seismic waveform is propagated, where ‘MODEL’ is, e.g., ‘RS’, ‘RS+pert’, or ‘QS+pert’. The time integrals in Eq. 4 are approximated by summations with a time discretisation of s. With this time discretisation travel times are accurate to s, which is sufficient for our purpose.
4.2 Estimating travel-time noise
While the numerical simulations performed in this paper are used to estimate the expectation value of a travel-time measurement (the signal), we also need to estimate the noise level in order to determine if a signal is detectable.
We use the travel-time maps measured by Liang et al. (2013) covering seven independent days of SOHO/MDI observations of AR9787 (Gizon et al. 2009; Moradi et al. 2010). In order to decrease noise, we further average these travel-time maps over a region behind the sunspot where scattering phase shifts are large. This elliptical portion of the travel-time maps is shown in Fig. 2. We use angle brackets to denote this spatial averaging of the travel times.
For one day of observations, the spatially averaged travel-time shift introduced by the sunspot is s for the p modes. In order to estimate the noise in these measurements, we compute for each of the seven days of observations. We define the observed noise level as the standard error in the mean of the spatially averaged travel-time shifts. The observed noise level decreases with radial order, it is s for the f modes, s for p modes and s for p modes (Table 1).
We carried out two sets of numerical experiments. In the first set, we propagated f, p and p waves independently through different sunspot models to predict the travel-time shifts caused by various perturbations to the reference sunspot model, and to determine if they are above the observed noise level. In the second set of experiments we looked at the sensitivity of p wave travel times to changes in the amplitude of the subsurface sound-speed. We considered sound-speed perturbations to both the quiet-Sun model and to the reference sunspot model. The results of these experiments are described in the sections below.
5.1 Is it possible to distinguish between sunspot models?
One way to measure the effect of a change in the sunspot model on the waves is to look at the change in the travel-time shifts measured with respect to the waves in the quiet-Sun reference model, i.e. by computing the difference .
Table 1 shows these changes in travel-time shifts for the perturbations to the reference sunspot model described in Sects. 2.3, 2.4 and 2.5. For the case where the Wilson depression is moved deeper (shallower) by km, we find that the f wave packet travels slower (faster) and the p and p wave packets travel faster (slower). Because the travel-time shifts for these two cases do not have the same amplitude, perturbation theory that is linear in the height of the Wilson depression will not fully capture this result. The physical cause of these effects is not clear. We note, however, that a change in the Wilson depression implies a change in the density, sound-speed and Alfvén speed which all affect the wave speed (Hanasoge et al. 2012). For the case of the p wave packet the change in the travel-time shift is roughly twice the noise level, implying that the p travel times are useful for constraining the height of the Wilson depression to an accuracy better than km (assuming the other sunspot parameters are well constrained).
Reducing the strength of the subsurface magnetic field causes a change in the travel-time shift of roughly s for all wave packets. The physical process governing the interaction of the wave with the subsurface magnetic field is not clear. Recall that this change in the magnetic field is only below Mm. As for the case of changes to the height of the Wilson depression, for the p wave packet the change in the travel-time shift is about twice that of the observed noise level. Assuming that the Wilson depression was well constrained (and all other sunspot parameters), the subsurface magnetic field may also be constrained.
Inserting subsurface sound-speed enhancements causes the p and p mode wave packets to travel faster than expected. The f wave packet is not sensitive to the subsurface sound-speed at all, which is also not surprising since it is a surface gravity wave. Unless the observed signal to noise can be increased for the p and p wave packets by a factor of at least , sound-speed perturbations of this type will not be able to be detected with the spatially averaged travel-time shifts of these modes.
|Sunspot Perturbation||Travel-time Shift (s)|
|50 km deeper|
|50 km shallower|
|Alternative magnetic field|
5.2 Can we avoid using numerical simulations?
In the second experiment, we further investigate the sensitivity of travel-time shifts to changes in the subsurface sound-speed. Fan et al. (1995), Kosovichev (1996) and Birch et al. (2004) explored this problem using perturbation theory around quiet-Sun reference models. Here we go beyond small amplitude perturbation theory and look at finite amplitude perturbations to the subsurface sound-speed. In addition, we look at the effects of identical sound-speed changes to a reference sunspot model and in the quiet-Sun model.
We insert Gaussian sound-speed perturbations of various amplitudes into the reference sunspot model (i.e. Eq. 3 with to and Mm). In order to assess the significance of the reference model, we also insert identical sound-speed perturbations into the background quiet-Sun model so that
where is the quiet-Sun sound speed, and is the perturbation to the reference sunspot model defined by Eq. 3. We simulated the propagation of p wave packets through these models.
Figure 3 shows the change in the spatially averaged travel-time shifts, , as a function of the amplitude of the sound-speed perturbation. The ensuing changes to the travel-time shift are not linear with the amplitude of the sound-speed perturbation. The deviation from linearity for a sound-speed perturbation with a 10% amplitude is 5% – a parabolic fit describes the change in the travel-time shifts with respect to sound-speed perturbation amplitude better. Notice that for all amplitudes, even up to 50%, the change in the travel-time shift is well below the observed noise level. Thus, these subsurface sound-speed perturbations are not detectable using only the p wave packet.
Another way to measure the travel-time shift due to the sound-speed perturbations is to directly compute the travel-time shift with respect to the reference sunspot model, i.e. (see Fig. 3). This gives a different, but qualitatively similar, result to the previous method. This shows that the general behaviour of the change in the travel-time shifts with amplitude is not just a consequence of the measurement procedure.
Figure 3 also shows the travel-time shifts, , caused by imposing the same sound-speed perturbations in the quiet-Sun model. In this case, the linear approximation is incorrect by more than 10% when the amplitude of the sound-speed perturbation is 10%. More importantly, these travel-time shifts are roughly four times larger than the travel-time shifts caused by imposing the same sound-speed perturbation in the reference sunspot model. As a result, the quiet-Sun model cannot be used to predict the travel-time shifts caused by sound-speed perturbations to a sunspot model because the sunspot substantially changes the p wave packet, and thus its sensitivity to subsurface perturbations. Therefore, it is necessary to use numerical simulations for a better understanding of the effects of a sunspot on the waves.
6 Summary and prospects
We have shown that it should be possible to constrain the height of the Wilson depression to within km and to distinguish between the two models for the subsurface magnetic field structure (a monolithic and a weaker spreading field) using seven days of SOHO/MDI observations of one sunspot. Using only the spatial average of the travel-times shifts it is not possible to detect 10% sound-speed changes with a length scale of 5 Mm at depths between and Mm with the same observations. We expect that including more information, e.g. the spatial distribution of the travel-time shifts, the amplitude of the cross-correlations, and higher order modes, will improve the ability of helioseismic measurements to constrain the subsurface structure of sunspots.
We have also demonstrated that the wave travel-time shifts caused by changes in sound speed to our sunspot model are different than the wave travel-time shifts caused by the same change in sound speed of the quiet-Sun model. This is due to the effect of the sunspot on the wavefield, and this effect is not captured by the traditional approximations of local helioseismology (Fan et al. 1995; Kosovichev 1996; Birch et al. 2004). Computational helioseismology, however, provides a clear path forward. These simulations may also be used to test the magnetoghydrodynamic ray approximation (e.g. Cally 2000; Schunker & Cally 2006).
Other avenues to explore are the use of vector magnetograms and intensity images to further constrain the surface properties of sunspots and the use of SDO-HMI Doppler observations (Scherrer et al. 2012) for helioseismic measurements.
Acknowledgements.H.S. and L.G. acknowledge research funding by Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1 “Astrophysical flow instabilities and turbulence” (Project A18, WP “Seismology of magnetic activity”). The German Data Center for SDO, funded by the German Aerospace Center (DLR), provided the IT infrastructure.
- Birch et al. (2004) Birch, A. C., Kosovichev, A. G., & Duvall, Jr., T. L. 2004, ApJ, 608, 580
- Cally (2000) Cally, P. S. 2000, Sol. Phys., 194, 189
- Cameron et al. (2007) Cameron, R., Gizon, L., & Daiffallah, K. 2007, Astronomische Nachrichten, 328, 313
- Cameron et al. (2008) Cameron, R., Gizon, L., & Duvall, Jr., T. L. 2008, Sol. Phys., 251, 291
- Cameron et al. (2011) Cameron, R. H., Gizon, L., Schunker, H., & Pietarila, A. 2011, Sol. Phys., 268, 293
- Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286
- Fan et al. (1995) Fan, Y., Braun, D. C., & Chou, D.-Y. 1995, ApJ, 451, 877
- Gizon & Birch (2002) Gizon, L. & Birch, A. C. 2002, ApJ, 571, 966
- Gizon et al. (2010) Gizon, L., Birch, A. C., & Spruit, H. C. 2010, ARA&A, 48, 289
- Gizon et al. (2006) Gizon, L., Hanasoge, S. M., & Birch, A. C. 2006, ApJ, 643, 549
- Gizon et al. (2009) Gizon, L., Schunker, H., Baldner, C. S., et al. 2009, Space Sci. Rev., 144, 249
- Hanasoge et al. (2012) Hanasoge, S., Birch, A., Gizon, L., & Tromp, J. 2012, Phys. Rev. Lett.
- Kosovichev (1996) Kosovichev, A. G. 1996, ApJ, 461, L55
- Liang et al. (2013) Liang, Z. ., Gizon, L., & Schunker, H. 2013, submitted to A&A, Paper I
- Lindsey & Braun (2005) Lindsey, C. & Braun, D. C. 2005, ApJ, 620, 1107
- Loughhead & Bray (1958) Loughhead, R. E. & Bray, R. J. 1958, Australian Journal of Physics, 11, 177
- Moradi et al. (2010) Moradi, H., Baldner, C., Birch, A. C., et al. 2010, Sol. Phys., 267, 1
- Rempel & Schlichenmaier (2011) Rempel, M. & Schlichenmaier, R. 2011, Living Reviews in Solar Physics, 8, 3
- Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
- Schunker et al. (2008) Schunker, H., Braun, D. C., Lindsey, C., & Cally, P. S. 2008, Sol. Phys., 251, 341
- Schunker & Cally (2006) Schunker, H. & Cally, P. S. 2006, MNRAS, 372, 551
- Schunker et al. (2011) Schunker, H., Cameron, R., Gizon, L., & Moradi, H. 2011, Sol. Phys.
- Schüssler & Rempel (2005) Schüssler, M. & Rempel, M. 2005, A&A, 441, 337
- Solanki et al. (1993) Solanki, S. K., Walther, U., & Livingston, W. 1993, A&A, 277, 639