Radiatively accurate magnetohydrostatic sunspot model

Spectropolarimetrically accurate magnetohydrostatic sunspot model for forward modelling in helioseismology

D. Przybylski, S. Shelyag, P.S. Cally Monash Centre for Astrophysics, School of Mathematical Sciences, Monash University, Clayton, Victoria, Australia, 3800.

We present a technique to construct a spectropolarimetrically accurate magneto-hydrostatic model of a large-scale solar magnetic field concentration, mimicking a sunspot. Using the constructed model we perform a simulation of acoustic wave propagation, conversion and absorption in the solar interior and photosphere with the sunspot embedded into it. With the magnetically sensitive photospheric absorption line of neutral iron, we calculate observable quantities such as continuum intensities, Doppler velocities, as well as full Stokes vector for the simulation at various positions at the solar disk, and analyse the influence of non-locality of radiative transport in the solar photosphere on helioseismic measurements. Bisector shapes were used to perform multi-height observations. The differences in acoustic power at different heights within the line formation region at different positions at the solar disk were simulated and characterised. An increase in acoustic power in the simulated observations of the sunspot umbra away from the solar disk centre was confirmed as the slow magneto-acoustic wave.

Subject headings:
Sun: magnetic fields - Sun: oscillations - Sun: helioseismology - sunspots

1. Introduction

Techniques of local helioseismology are currently unable to unambiguously determine sub-surface structure of the flows and sound speed perturbations in and around large-scale solar magnetic field concentrations, such as sunspots and pores (Shelyag et al., 2007a; Gizon et al., 2009; Moradi et al., 2010). Due to the complexity of magnetohydrodynamic processes involved, our understanding of the behaviour of magnetoacoustic waves as they are absorbed, reflected and refracted by sunspots is far from complete.

Recently, it was demonstrated that solar magnetic fields and the process of magneto-acoustic wave mode conversion associated with them lead to significant changes in the wave travel times used in helioseismic inversions (Moradi & Cally, 2014; Hansen & Cally, 2014). Four processes associated with strong magnetism dominate wave behaviour in sunspots (Cally et al., 2015): fast/slow mode conversion at the Alfvén/acoustic equipartition level , allowing acoustic (slow) waves to transmit into the upper atmosphere if the ‘attack angle’ between wave vector and magnetic field is small but converting them to magnetic (fast) waves otherwise (Cally, 2006; Schunker & Cally, 2006); the “ramp effect” that reduces the effective acoustic cutoff frequency to , where is the magnetic field inclination from the vertical (Bel & Leroy, 1977); fast wave reflection around the height where the Alfvén speed matches the wave’s horizontal phase speed; and fast/Alfvén mode conversion that typically occurs over several scale heights near the fast wave reflection level, generating both upward and downward propagating Alfvén waves (Cally & Hansen, 2011). Fast/slow conversion is found to produce large negative travel time shifts, while fast/Alfvén conversion generates countervailing positive shifts provided the vertical plane containing the wave vector is nearly perpendicular to the vertical plane containing the magnetic field lines (Cally & Moradi, 2013). The ramp effect allows field-guided acoustic waves to enter the atmosphere in inclined field where frequency while normally the acoustic cutoff would prevent their propagation (). The complexity and sensitivity to magnetic field direction of wave motions above the equipartition level makes interpretation of observations difficult but potentially rewarding. Moradi et. al. (2015, submitted) studied the effects of directional time-distance helioseismology on the travel time measurements in the sunspot model.

There are also possible significant discrepancies in travel time measurements originating from the effects of non-locality of radiative transport in the solar atmosphere. Changes in spectral line formation heights due to magnetic field presence (see e.g. Shelyag et al., 2007b), systematic centre-to-limb variations in absorption line formation (Shelyag & Przybylski, 2014), as well as instrumental effects, such as stray light, and other processes involved in formation and measurement of radiation intensities and Doppler shifts result in our inability to unambiguously measure the travel time perturbations and, therefore, infer solar sub-surface structure (Rajaguru, 2011).

Rapid improvements in computational power already make it possible to perform forward modelling of magnetohydrodynamic wave propagation and mode conversion in “realistic” solar magnetic field structures (Shelyag et al., 2009; Moradi et al., 2009; Felipe et al., 2010; Cameron et al., 2011; Khomenko & Cally, 2012; Felipe, 2012; Zharkov et al., 2013; Felipe et al., 2014). Spectral line synthesis codes and radiative diagnostics tools also allow computations of mock observables from the simulated plasma parameters, allowing for direct comparison between simulations and observations in computational helioseismology.

Creating a sunspot that is both spectropolarimetrically accurate and magnetohydrostatically stable is inherently difficult, as the sound speed and temperature can change significantly with small changes in the density and pressure stratification. The sunspot model of Khomenko & Collados (2008) was created to allow empirical quiet and umbral solar models to be used in the near-surface layers in combination with a Schlüter-Temesvary flux tube model (Schlüter & Temesváry, 1958) in the interior. However, the model created this way is still not convectively stable. Convective instability is fatal to linear MHD simulations, but these codes are less expensive than full non-linear simulations, and ideal for the long time series required in helioseismology; for the study of fast and slow magneto-acoustic waves; and for simulating fast-slow and fast-Alfvén mode conversion in the photosphere and lower chromosphere. The effects of convective stabilisation on the eigenmodes of solar models for helioseismic simulations were studied by Schunker et al. (2011). A technique for stabilising the atmosphere is discussed in Sec. 2.

In this paper, we present a model of a magneto-hydrostatic and spectropolarimetrically accurate sunspot. Our model is based on the sunspot-like model of Khomenko & Collados (2008). The model was adjusted to provide a more accurate replication of photospheric sunspot properties taken from semi-empirical models, while still maintaining a smooth transition of physical properties between the magnetic and non-magnetic regions required for stable numerical simulation. This technique makes it possible to obtain accurate photospheric absorption line formation heights as well as allowing the study of observational signatures of acoustic wave propagation in the simulated model at different positions on the solar disk. We perform a magnetohydrodynamic simulation of the propagation of a wave through this sunspot-like model and investigate the behaviour of acoustic waves in the simulated model using the synthesised radiation, as if it were observed. We also investigate effects of the centre-to-limb variation effects on Doppler velocity measurements and study the line bisector shapes to allow for a multi-height view in the line formation region, which can be used to observationally disentangle wave mode conversion process in the solar atmosphere.

The structure of the paper is as follows. In Section 2 we describe the background model. In Section 3 we explain the magnetohydrodynamic simulation and the spectral synthesis methods used to provide artificial observables. Section 4 provides results and description of the radiative effects on acoustic wave measurements in observations. In Section 6, we discuss our findings.

2. Model

To provide a convectively stable quiet Sun background model, the method described by Parchevsky & Kosovichev (2007) is used. The Brunt-Väisälä frequency must be positive for convective stability. Rearranging this equation in terms of , combining with the equation for hydrostatic stability (Equation 1) and introducing a free parameter gives,


where the gravity acceleration , Brunt-Väisälä frequency , and the adiabatic index are functions of depth. The equations are solved using a fourth-order Runge-Kutta method on a one-dimensional grid. The equispaced grid covers the height range from to and is resolved with the vertical step .

To enforce convective stability, the Brunt-Väisälä frequency dependence on is modified by setting the negative values in the convectively unstable solar interior to small positive ones. The free parameter must be greater than zero to enforce convective stability and is increased so to match the pressure in the model with the Standard Solar Model S (Christensen-Dalsgaard et al., 1996) pressure at a point in the interior . The quiet Sun model is created by smoothly joining the VALIIIC (Vernazza et al., 1981) to Model S and integrating from the top of the solar atmosphere downwards. Using the above, a convectively stable quasi-solar model can be created with negligible change to the photospheric line formation regions and the correct sound speed profile below the surface at the expense of a slightly reduced density and pressure in the deeper regions.

Following the method described by Khomenko & Collados (2008), an axisymmetric sunspot model is created in cylindrical geometry using three parameters , and which change the sunspot radius, magnetic field inclination and strength, respectively. A full description of the effects of these parameters on the magnetic field configuration is given by Khomenko & Collados. The model is defined on a two-dimensional - plane discretised into a domain from to in height, with a radius of and resolution of and .

Below Mm depth a Low-type magnetic flux tube (Low, 1980) is constructed using an extension of the exact Schlüter-Temesvary formulation (Schlüter & Temesváry, 1958).

For the near-surface layers and in the atmosphere a Pizzo-type magnetic flux tube is used (Pizzo, 1986). The Pizzo method creates a pressure-distributed magnetic field structure through an extension of the Low formulation used above. The magnetohydrostatic equation for a non-twisted, cylindrical structure can be simplified by introducing a magnetic vector potential (Low, 1975). This allows the problem to be reduced to a single equation for a scalar (Pizzo, 1986)


where is the gas pressure along the field lines.

The Pizzo method boundary conditions require both a quiet Sun (denoted with index ) and umbral (denoted with index ) pressure, density, temperature, and pressure scale height () and temperature distributions as functions of depth. The quiet Sun model (,,) generated above was used for the outer boundary condition. For the inner boundary the Avrett semi-empirical model (Avrett, 1981) is used, which is then joined to the pressure and density profiles at the axis of the self-similar flux tube using log-linear interpolation. This is then convectively stabilised using Equations (2) and (1) as described for the quiet Sun above. The Wilson depression can be prescribed by shifting the of the umbral model (, , ). The pressure and scale height are then distributed throughout the domain using the following:


The potential solution given by Equation (3) is used as an initial guess. The pressure distribution given by Equations (4) and (5) is iterated together with Equation (3) using a Gauss-Seidel method. Thus, the complete force balance is calculated with a specified precision, giving a final distribution of the potential and pressure.

The Pizzo and Low type flux tubes are then joined at and recalculated using Equations (4-5). The density and the radial and vertical components of the magnetic field vector are calculated according to:


To extend this model below a vertical flux tube with a constant and zero is used, and the pressure and density profiles are continued smoothly downwards.

Finally, the FreeEOS equation of state (Irwin, 2012) is applied to find the adiabatic index, temperature and sound speed at each grid cell in the model. The model is then converted to Cartesian geometry, giving the full set of physical parameters required for the MHD simulations and radiative transfer calculations.

Using the procedure explained above the magnetic field structure pictured in Figure 1 was constructed. The background image in the figure shows the modulus of magnetic field . The field lines are nearly vertical in the “umbral” region (), and show inclination of about in the “penumbral” region, , of the sunspot model.

In the figure, the dashed line shows the layer, while the dotted contours represent levels. As is evident from the figure, in the umbral region at the axis of the sunspot, the layer is positioned higher than layer, suggesting formation of the continuum radiation in the magnetically-dominated sunspot atmosphere.

Figure 1.— Magnetic field structure of the sunspot model. The magnetic field strength is shown with the magnetic field lines overplotted (solid). Also shown are the (middle dotted), 0.1 (upper dotted) and 10 (lower dotted) contours. The dashed line is the contour, representing the visible photosphere. Note that the aspect ratio is severely stretched.

The 6173Å photospheric absorption line of neutral iron is used for observations of the full solar disk with the Helioseismic Magnetic Imager (HMI) onboard the Solar Dynamic Observatory (SDO). Therefore, this line was chosen to carry out radiative diagnostics of the sunspot model using the SPINOR code (Solanki, 1987; Shelyag et al., 2007b). For each one-dimensional column of the model, continuum intensity and spectral line profile calculations are performed by solving the Unno-Rachovsky (Unno, 1956) radiative transfer equation for the Stokes vector . Off-disk centre observations are simulated by inclining the numerical domain and interpolating the density, temperature, magnetic field and velocities onto the new line of sight (los). The slanting is performed around the height and in the direction of positive (Figure 2). The velocity and magnetic field vectors are then projected into the new reference frame. The calculation uses 500 wavelength points with a to ensure the spectral line is highly resolved. The los velocity is given by . The magnetic field is recalculated using a similar relation.

Figure 2 shows the continuum images of the sunspot model calculated for , and angles between the surface and the los, which correspond to viewing cosine , respectively. We find that the model produces a realistic limb darkening dependence with a continuum value of 79% of the disk centre intensity at . This is only slightly higher than the 75% of the limb darkening curve determined by Foukal et al. (2004).

Figure 2.— continuum intensity of the sunspot model at the observational angles (top to bottom) , and to the vertical. The figures have been normalised to the quiet Sun value at inclination. Inclination is performed towards an observer displaced in the negative direction. The crosses show the two penumbral, and one umbral point used in Figure 3
Figure 3.— los velocity response functions of Stokes profile of Fei line computed for the models of quiet Sun (first column), near-side and far-side penumbrae (second and third columns, respectively), and umbra (fourth column) for positions at the solar disk. The y-axis represents the height along the los at which the perturbation is placed, where represents the layer for the quiet Sun photosphere. The optical depth axis is also shown, with a dashed line showing the observed depression of the photosphere. The corresponding Stokes- (top right) and Stokes- (top-left) profile shapes were over-plotted in white in each panel over the wavelength range shown in the figure.

The velocity response functions of the spectral line are shown in Figure 3 for the quiet Sun, two penumbral regions at , and in the centre of the sunspot umbra for the chosen positions at the solar disk. These locations have been marked with crosses in Figure 2. Since, for observations away from solar disk centre, two points at the same distance from the sunspot axis are not equivalent, the penumbral models have been chosen so that los of P1 crosses the umbral region, while the los of P2 inclines further into the penumbra. The response functions were calculated by computing a perturbed profile with a small positive (directed towards the observer) los velocity perturbation and subtracting from it an unperturbed profile for the same location.

The top row of Figure 3 shows the response functions of the four points in the model for the disk centre calculation. The perturbation is directed towards the observer, causing the line to be blue-shifted. The lobes of the response function tilt inwards towards the line core, marked as in the plots, clearly demonstrating dependence of the sensitivity of the profile on height; the regions closer to the line-core (on the -axis) are formed higher in the atmosphere.

The top-right figure shows a fully Zeeman-split profile in a strong umbral magnetic field. Notably, while the line formation height range is narrower compared to the quiet Sun, the response function lobes are wider in wavelength, suggesting higher sensitivity of the line to velocity perturbations.

The two penumbral points are identical in the solar disk centre simulation due to the symmetry of the sunspot. As the penumbral magnetic field is weaker, the line is not completely split. In the line core, the response function shows two smaller regions of sensitivity to the velocity perturbation.

Figure 3 demonstrates that the line formation height range increases with the inclination angle. In the case of the quiet Sun (left column of the figure), it increases from at the disk centre to at . As the line width does not change significantly, the wavelength range of the response function does not change with the inclination.

For the cases of magnetised penumbral and umbral atmospheres, the observed visible sunspot surface increases with the inclination angle. Between and 0.5, the level for the penumbral points P1 and P2 is shifted downwards by , and by for the umbra. The line sensitivity height range also increases further away from the disk centre, similarly to the quiet Sun.

Notably, the line profiles and the response function shapes for P1 and P2 are very different. The far-side umbra (second column of Figure 3, P1 in Figure 2) will have a formation range that extends into the highly magnetic umbra. This can be observed as an increasingly split profile as inclination increases in the 2nd and 3rd rows. The near-side penumbral pixel (third column, P2) will similarly form in a region of lower magnetic field. Due to the inclination of the magnetic field, P1 will measure a higher magnetic field strength along the los while P2 will incline against the direction of field line inclination.

The angle of the two ridges is seen to be larger in the umbral distributions than the quiet sun. For a small wavelength range in the quiet sun up to of the atmosphere will be measured. Comparitively, in the umbral distribution a similar filter would only sample around .

The large range of responses seen in the different penumbral and umbral positions will lead to larger uncertainty in the observation height of velocity measurements. The impact of Zeeman-split profiles on velocity measurements is only amplified at higher inclinations.

3. Numerical Simulations

We perform a simulation of acoustic wave propagation in the simulated model with the SPARC code. The code was designed to solve the linearised ideal MHD equations for wave propagation in a stratified solar environment (Hanasoge, 2011). The version of the code we employ for the simulations uses Message Passing Interface (MPI) to parallelise the computation and reduce computation time. It uses an implicit compact 6th order finite difference scheme applied to the horizontal and vertical derivatives. An explicit filter is used to prevent numerical instabilities in the solution. The boundary conditions used in the simulation include a Perfectly Matched Layer (PML) (Hanasoge et al., 2010) at the top and bottom boundaries, allowing for efficient absorption of the outgoing waves. A ‘sponge’ type absorbing layer is used on the side boundaries, which adds a linear friction term to the governing equations (Colonius, 2004). The code includes a Lorentz force limiter (Rempel et al., 2009), which is required due to the high Alfvén speed in the solar atmosphere. The limiter, although unphysical, prevents reduction of the time step and excessive computational times. The Alfvén speed limiter is set at , which is sufficiently high to allow fast MHD waves of interest, which have horizontal phase speed less than this, to propagate and refract correctly while still allowing us to use a manageable time step. Our cap is large enough for this to not be an onerous constraint. The implications of the limiter on helioseismic travel time shifts have been studied in detail by Moradi & Cally (2014).

The numerical grid has horizontal extent of grid points, with a physical size of , giving a resolution of in the horizontal directions. To deal with the large variation in physical parameters over the domain from the convective zone to chromosphere, the code uses a vertical grid spacing based on the sound speed. The grid has points between and and is distributed such that the acoustic travel times between each cell are the same for the quiet Sun, . This gives a resolution of around near the photosphere, and around in the lower solar interior. This means we do not resolve slow waves in the large regime, but these are effectively decoupled from the system anyway, so their neglect is not important. The following acoustic source, similar to that described by Shelyag et al. (2009), was used:


where s, s, s, Mm, . The position of the pulse is , .

The SPARC code solves the MHD equations for the perturbations around the MHS background model. A master-slave Open-MPI code has been written to take these perturbations, combine them with the background model and incline them as required. The SPINOR routines are then applied to each pixel to generate the full Stokes vector for each pixel. Using the generated Stokes- profiles, the los centre-of-gravity Doppler velocity is calculated by computing the position of the centre of gravity of the line profile and determining its shift from the unperturbed counterpart, computed for the background model, according to:


To calculate bisector Doppler velocities from the spectral line the relative intensity was determined by normalising the measured Stokes I between 0 and 1. Bisectors of the spectral line were calculated at 100 evenly spaced values between of . The bisectors were calculated for the background model and for each output snapshot. A Doppler velocity was then determined for each snapshot using the shift from the unperturbed background value, according to:


4. Results

Using the model described in Section 2 and methodology given in Section 3, a hour simulated observation of wave propagation through the sunspot model was performed. The top panel in Figure 4 shows a time-distance plot of the centre-of-gravity los Doppler velocity measured using Equation (10). The first three wave bounces can be easily seen. A shift in the wave arrival time can be observed as a flattening of the wavefront as it passes through the sunspot umbra at . The middle panel of Figure 4 shows a time-distance plot of the vertical component of velocity at the level of the simulation domain. A comparison between the top two panels shows that the vertical velocity in the domain and the los Doppler velocity are visually identical. Some reflection can be seen from the top and bottom PMLs, and from the side boundaries.

The bottom panel of Figure 4 shows the horizontal component of velocity, scaled by to provide a view of the slow magnetoacoustic wave in the strong magnetic field. A slice is taken through the centre of the simulated sunspot . The fast wave can be seen to propagate through the sunspot in the lower interior where plasma- is high. At around in the umbra, the incoming fast wave hits layer (See Figure 1) and undergoes partial transmission as a slow mode (effectively acoustic in ). The slow magnetoacoustic wave (now magnetic in ) can be seen to propagate back down into the sunspot as a flattening banding in the time distance plot. The wave amplitude in the atmosphere is low due to scaling by the very low densities, however, it still can be seen to continue to travel upwards above the photosphere and escapes through the absorbing upper boundary.

Figure 4.— Response of the model to the acoustic source. Top panel - simulated Doppler velocities at inclination, measured at . Middle panel - simulated vertical component of velocity at a geometric height , measured at . The differences between these two can be seen to be small. Bottom pannel - through spot centre (0,0) Mm showing the propagation of slow modes down through the box. Two dashed horizontal lines in the top plots mark the position of the sunspot umbra.

Figure 5 shows a power spectrum plotted with azimuthally averaged wavenumber and frequency. The spectrum has been calculated from the hour time-series of centre-of-gravity velocities calculated from the synthesised line profile. The power ridges are well resolved, although there are gaps in power at and shifts seen for high . Similar gaps are found in the simulations of Parchevsky & Kosovichev (2007) with high top boundary (1.75 Mm; their Fig. 6c), and attributed to trapping of acoustic modes. We do not understand how acoustic trapping explains this phenomenon. The gaps are also present in quiet sun simulations (no magnetic field), but are largely removed when the top boundary is lowered to 500 km above the solar surface (their Fig. 6b). This suggests that there is some numerical dissipation mechanism operating in our model chromosphere () that we are yet to identify.

Figure 5.— A power spectrum for the sunspot box calculated from synthetic Doppler velocities.

Using the Fei spectral line data, velocity shifts were calculated for each bisector as the differences between the perturbed and background values (Equation (11)). From the 2.5 hour time series of these velocities acoustic power maps of wave propagation in the sunspot were created. The power maps are calculated for 100 bisector depth positions by taking the fast Fourier transform of the velocity time series for each of the pixels in the model. The acoustic power is binned into different frequency bins by applying a Gaussian filter with a FWHM of around the chosen central frequency. For each frequency band and inclination the acoustic power is normalised by its average in the simulation domain.

By taking the Doppler velocity measurements at multiple bisector depths, it was possible to disentangle the dependence of the wave behaviour on the height within the line formation region. As Figure 3 shows, bisectors taken higher in the line profile are formed deeper and closer to the continuum formation height, at a physical height of around , while bisectors taken deeper are closer to the absorption line core and formed higher at a physical height of around . As was already noted, these formation heights vary significantly in the penumbra and umbra of the sunspot, with an offset of due to the prescribed Wilson depression.

First we aim to investigate the horizontal distribution of acoustic power throughout the simulation domain. The region of the acoustic source at , has been masked in the power maps. Figures 6 and 7 show the acoustic power calculated from the Doppler shifts measured at the bisector positions of and , respectively. The acoustic power was binned into frequency 1 mHz bands centred at 3, 4, 5, 6 and 7 mHz (left to right in the figures), and the disk positions of , and were used (top to bottom). los velocity measurements off the disk centre are affected by a larger line formation region and the presence of horizontal velocity components in the los velocities. Immediately obvious in the figure is a series of concentric rings of power travelling out from the source. These rings can be thought of as the spatial analogue of the ridges seen in Figure 5 and occur at discrete values because a single point-like source was used in the simulation.

The differences between the acoustic power maps in Figure 6 and 7 represent changes in the wave behaviour over the formation height of the spectral line. This formation height can span almost a megameter close to the solar limb (see Figure 3), and the differences are substantially more pronounced in the magnetic regions.

In the acoustic power maps (Figure 6 and 7) a solid line representing the sunspot umbra is over-plotted. Outside the sunspot umbra the sunspot shadow is observed at all frequencies. In the sunspot shadow, the power from the pulse has been absorbed or reflected by the magnetic field structure. This is most obvious at lower frequencies (left two columns). As frequency increases, the ridges of increased acoustic power can still be seen behind the sunspot. Interestingly, the magnetic field perturbs the concentric rings seen at 7 mHz, as the fast wave propagation speed and turning height change. Behind the sunspot in the band (right column), a small region of increased power can be seen between the two outermost rings of the power ridges (around ). As this is seen in high frequencies and as a ring around the sunspot, rather than around the source, it appears to be a far-side acoustic halo. Acoustic halos are seen around active regions as an excess of acoustic power compared to the quiet sun. 111A more comprehensive look at physics of acoustic halos, based on three-dimensional simulations of this sunspot model can be found in Rijs et al. (2015), which expands on the previous works of Hanasoge (2008); Khomenko & Collados (2009). They are attributed to fast MHD waves reflected in active region atmospheres by the steep Alfvén speed gradients there.

Comparison of the weakly-magnetic regions of each column in Figure 6 shows little variation with inclination, regardless of frequency. In these regions the propagating fast wave dominates the simulated observations as there is little magnetic field to alter the fast-wave turning point or to allow for mode-conversion process to take place. Inside the umbral region, the acoustic power maps for the disk centre case show very little power in 3 and 4 mHz, and the power in the vertically-aligned oscillations is seen to be almost completely absorbed. As the inclination increases, the 3 and 4 mHz power in the umbra remains low, while the 5, 6 and 7 mHz power bands show a significant power enhancement.

From the response function (first column in Figure 3) we expect there to be larger differences in power between the line core (Fig. 6) and line wings (Fig. 7) as we observe further away from the solar disk centre. There is little difference found in the disk centre cases (top row) between the two figures. However, at inclination significant differences can be seen between the power maps at the line core and line wings. Particularly, (1) the umbral power increase is only seen at the line core (Fig. 6); (2) the ring-like structure (marked by the dashed circle in the left panel of Fig. 6) found at around in the bottom left two panels is somewhat wider at the line core than in the line wings (Fig. 7). This structure is most apparent in the frequency bands (first two columns), and the power in it decreases with increasing frequency. While the inclination of the magnetic field at the surface at the radius of is degrees from vertical the magnetic field strength is low, and the equipartition layer is located above the line formation region. Therefore, the ring is of acoustic nature and cannot be related to the slow magneto-acoustic mode, as it is found at the source side in both and inclinations corresponding to the los direction which is either parallel or highly inclined to the magnetic field.

Figure 6.— Acoustic power map calculated from the shifts in the bisector of the Fei 6173Å line at a bisector height of . The columns, from left to right, show power in the 3, 4, 5, 6 and 7 mHz bands. The rows, from top to bottom, show measurements made at , and inclination from the vertical, where the field of view has been inclined in the -direction. The power at each inclination angle has been normalised by its average in each frequency band. This represents a velocity measurement made near the line-core, showing the behaviour higher in the atmosphere. The sunspot umbra has been marked with a solid circle, while the dashed line in the bottom right panel represents the slice taken in Figure 8. The low frequency ring has been marked in the top left panel.
Figure 7.— Acoustic power map calculated from the shifts in the bisector of the Fei 6173Å line at a bisector height of 0.9  . The layout of the columns and rows is as seen in Figure 6. This figure represents measurements made near the line wings.

As demonstrated, the umbral and penumbral acoustic power structures are mostly seen near the line core (Figure 6). In the los velocities measured from bisectors in the line wings (Figure 7) only a faint structure can be seen at high inclinations, again more obvious in the 7 mHz power band (bottom right).

Figure 8.— Bisector power map at x = 0 Mm. The layout of the columns is as in Figure 6. The rows, from top to bottom, show measurements made at , and inclination from the vertical. The Y axis in this figure represents the line depth where the bisector wavelength and Doppler velocity are measured and covers a large part of the line formation region of in length. The formation region will depend on the magnetic field strength and inclination (Figure 3).

To better understand the three-dimensional structure of the umbral power increase, in accordance with the response functions shown in Figure 3, a multi-height observation is made by computing the Doppler velocities from bisector shifts measured at different line depths within the line formation region. In Figure 8, a Doppler velocity map for a slice marked by a dashed line in Figure 7 is plotted for each bisector depth in the range .

For the and inclinations (top and middle rows of Figure 8), the changes in the structure and magnitude of acoustic power over the line formation range are limited to an increase in power in the high frequency bands for observations made closer to the line core. This matches the observations of acoustic halos by Rajaguru et al. (2013), where the acoustic power was weaker using filtergrams close to the line-wings. At inclination (bottom row, left two columns) the faint radius low frequency ring can be made out, increasing in radius at larger heights.

The vertical extent of the umbral power structure, seen at high inclinations at 6–7 mHz (bottom right two panels in Figure 68) shows a significant increase in power higher in the formation region. The formation of this acoustic power structure seems to start mid-way up the line formation region, suggesting a highly localised phenomenon. The inclination of the field lines at the centre of this power increase () is approximately from the vertical. Taking into account the observation angle of , the field line is almost perpendicular to the line-of-sight. Comparison with Fig. 1 shows the continuum formation height of the spectral line and the layer cross at and radius from the sunspot centre, allowing for direct observation of conversion of fast (parallel to the magnetic field line, perpendicular to the los) waves to slow (perpendicular to the magnetic field line, and parallel to the los) waves. The increased power in this region corresponds to the slow magneto-acoustic wave in the region where the magnetic field is close to perpendicular to the los. This result confirms previous findings by Zharkov et al. (2013) that the observed acoustic power increase in the sunspots is a signature of slow magneto-acoustic waves.

5. Discussion and Conclusion

In this paper we have: (1) described a modification of the Khomenko sunspot model we developed to provide a more accurate line formation region, allowing for accurate spectral synthesis to be performed; (2) analysed response functions in our model to understand our synthesised centre-to-limb observations of the Fei spectral line in a model sunspot; (3) have investigated wave propagation through the model using a three-dimensional simulation of magnetoacoustic wave propagation; (4) computed the spectral line profiles and provided a time series of the simulated observations of the sunspot model using the simulation data; (5) used spectral line bisector measurements to perform multi-height simulated observations over the line formation region; (6) studied maps of the acoustic power in the sunspot to better understand the absorption line response to the oscillations in sunspots and the effects of non-locality of radiative transport on helioseismic measurements.

The sunspot is found to absorb or scatter the incoming acoustic wave energy in all but the 6 and 7 mHz frequency band. There are signatures of a slight power enhancement seen around the sunspot, similar to those seen in observations of acostic halos.

Small ridges of acoustic power can be seen in the 7 mHz band in the shadow of the sunspot. Zharkov et al. (2013) showed in a simple 2D simulation that there was a significant power enhancement as a far-side acoustic halo. We can see slight power enhancements in the range of 25-40 Mm from the sunspot. Comparing the different rows in Figures 6 and 7 shows that outside the sunspot umbra, the horizontal and vertical velocities behave similarly, with only minor differences in the power maps.

The appearance of a high-frequency anomalous acoustic power excess in the sunspot centre – the umbral “belly button” – can be seen predominantly in the case of inclination with faint traces at lower inclinations. This geometry suggests that it is driven by slow magneto-acoustic waves, as it is seen when the horizontal velocity component dominates the los velocity (bottom rows Figure 6 and 7). It can be seen as a crescent-like structure, which is most dominant towards the line core (higher in the atmosphere, Figure 6) and very faint in the spectral line wings (lower in the atmosphere, Figure 7). Umbral power enhancements are seen in the space-based (HMI, Zharkov et al. (2013)) and ground based (Balthasar et al., 1998) observations of sunspots. Notably, no power excess was observed in the G-band (Nagashima et al., 2007). The appearance of this power increase in a simulated sunspot suggests a magneto-acoustic phenomenon, rather than photon noise (Donea & Lindsey, 2015). There are many differences in both sunspot properties and the radiative effects on HMI measurements that could explain the lack of the umbral power increase in all acoustic observations of sunspots. These include changes in the Wilson depression, the wide range of the velocity-response function in a magnetic structure (Figure 3), low resolution and high-noise measurements off the disk-centre or issues with using discrete filters on highly split profiles.

Current measurements of acoustic travel times in computational helioseismology are largely performed using measurements at the geometric heights in the simulation domain (Moradi & Cally, 2013), or on a surface roughly representing the continuum formation height determined by the layer in the simulation domain (Khomenko & Cally, 2012; Moradi et al., 2015). Despite the fact that the physical velocity in the simulation matches reasonably well to the los velocities calculated from the simulated spectral lines, this method misses a lot of information that can otherwise be gained from the range of formation of the spectral lines. As we show, this range also changes substantially if the simulation is performed for positions at the solar disk away from the centre. Using the model we described, artificial observables mimicking the HMI and MDI pipelines can be made (Scherrer et al., 2012; Fleck et al., 2011), as well as comparisons to ground based observations.

The multi-height Doppler measurements made by Nagashima et al. (2014), using the HMI filter-grams provide a similar approach to multi-height measurements as the bisectors used in this study. Rather than velocity measurements made using shifts in Stokes-, HMI uses measurements of both Stokes and Stokes . As the next step, it will be important to fully simulate the HMI data pipeline response to a variable magnetic atmosphere before a direct comparison to the HMI measurements can be made.

Sergiy Shelyag’s research is supported by Australian Research Council Future Fellowship. Damien Przybylski and Sergiy Shelyag’s research is performed with the computational resources provided by the Multi-modal Australian ScienceS Imaging and Visualisation Environment (MASSIVE) (www.massive.org.au) and through Astronomy Australia Limited by the NCI National Facility systems at the Australian National University and the Centre for Astrophysics and Supercomputing of Swinburne University of Technology (Australia).


  • Avrett (1981) Avrett, E. H. 1981, The Physics of Sunspots, 235
  • Balthasar et al. (1998) Balthasar, H., Martínez Pillet, V., Schleicher, H., Wöhl, H. 1998, Sol. Phys., 182, 65
  • Bel & Leroy (1977) Bel, N., & Leroy, B. 1977, å, 55, 239.
  • Cally (2006) Cally, P. S. 2006, Royal Society of London Philosophical Transactions Series A, 364, 333
  • Cally & Hansen (2011) Cally, P. S., & Hansen, S. C. 2014, ApJ, 738, 119
  • Cally & Moradi (2013) Cally, P. S., & Moradi, H. 2013, MNRAS, 435, 2589
  • Cally et al. (2015) Cally, P. S., Moradi, H., & Rajaguru, S. P. 2015, in Low Frequency Waves in Space Plasmas, ed. A. Keiling, Dong-Hun Lee, K.-H. Glassmeier and V. Nakariakov, American Geophysical Union Geophysical Monograph Series, submitted.
  • 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
  • Colonius (2004) Colonius, T. 2004, Annual Review of Fluid Mechanics, 36, 315
  • Felipe et al. (2010) Felipe, T., Khomenko, E., Collados, M., & Beck, C. 2010, ApJ, 722, 131
  • Felipe (2012) Felipe, T. 2012, ApJ, 758, 96
  • Felipe et al. (2014) Felipe, T., Socas-Navarro, H., & Khomenko, E. 2014, ApJ, 795, 9
  • Foukal et al. (2004) Foukal, P., Bernasconi, P., Eaton, H., & Rust, D. 2004, ApJ, 611, L57
  • Fleck et al. (2011) Fleck, B., Couvidat, S., & Straus, T. 2011, Sol. Phys., 271, 27
  • Gizon et al. (2009) Gizon, L., Schunker, H., Baldner, C. S., et al. 2009, Space Sci. Rev., 144, 249
  • Hanasoge et al. (2010) Hanasoge, S. M., Komatitsch, D., & Gizon, L. 2010, A&A, 522, A87
  • Hanasoge (2008) Hanasoge, S. M. 2008, ApJ, 680, 1457
  • Hanasoge (2011) Hanasoge, S. M. 2011, Astrophysics Source Code Library, 5006
  • Hansen & Cally (2014) Hansen, S. C., & Cally, P. S. 2014, Sol. Phys., 289, 4425
  • Irwin (2012) Irwin, A. W. 2012, Astrophysics Source Code Library, 11002
  • Khomenko & Collados (2008) Khomenko, E., & Collados, M. 2008, ApJ, 689, 1379
  • Khomenko & Collados (2009) Khomenko, E., & Collados, M. 2009, A&A, 506, L5
  • Khomenko & Cally (2012) Khomenko, E., & Cally, P. S. 2012, ApJ, 746, 68
  • Donea & Lindsey (2015) Donea, A.-C., Lindsey, C., 2015, ApJ, submitted
  • Low (1975) Low, B. C. 1975, ApJ, 197, 251
  • Low (1980) Low, B. C. 1980, Sol. Phys., 67, 57
  • Moradi et al. (2009) Moradi, H., Hanasoge, S. M., & Cally, P. S. 2009, ApJ, 690, L72
  • Moradi et al. (2010) Moradi, H., Baldner, C., Birch, A. C., et al. 2010, Sol. Phys., 267, 1
  • Moradi & Cally (2013) Moradi, H., & Cally, P. S. 2013, Journal of Physics Conference Series, 440, 012047
  • Moradi & Cally (2014) Moradi, H., & Cally, P. S. 2014, ApJ, 782, LL26
  • Moradi et al. (2015) Moradi, H., Cally, P.S., Przybylski, D., et al. 2015, MNRAS, submitted
  • Nagashima et al. (2007) Nagashima, K., Sekii, T., Kosovichev, A. G., et al. 2007, PASJ, 59, 631
  • Nagashima et al. (2014) Nagashima, K., Löptien, B., Gizon, L., et al. 2014, Sol. Phys., 289, 3457
  • Parchevsky & Kosovichev (2007) Parchevsky, K. V., & Kosovichev, A. G. 2007, ApJ, 666, 547
  • Pizzo (1986) Pizzo, V. J. 1986, ApJ, 302, 785
  • Rajaguru (2011) Rajaguru, S. P. 2011, Astronomical Society of India Conference Series, 2, 181
  • Rajaguru et al. (2013) Rajaguru, S. P., Couvidat, S., Sun, X., Hayashi, K., & Schunker, H. 2013, Sol. Phys., 287, 107
  • Rijs et al. (2015) Rijs, C., Moradi, H., Przybylski, D., & Cally, P. S. 2015, arXiv:1501.01074
  • Rempel et al. (2009) Rempel, M., Schüssler, M., & Knölker, M. 2009, ApJ, 691, 640
  • Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  • Schlüter & Temesváry (1958) Schlüter, A., & Temesváry, S. 1958, Electromagnetic Phenomena in Cosmical Physics, 6, 263
  • Schunker & Cally (2006) Schunker, H., & Cally, P. S. 2006, MNRAS, 372, 551
  • Schunker et al. (2011) Schunker, H., Cameron, R. H., Gizon, L., & Moradi, H. 2011, Sol. Phys., 271, 1
  • Shelyag et al. (2007a) Shelyag, S., Erdélyi, R., & Thompson, M. J. 2007, A&A, 469, 1101
  • Shelyag et al. (2007b) Shelyag, S., Schüssler, M., Solanki, S. K.,  Vögler, A. 2007, A&A, 469, 731
  • Shelyag et al. (2009) Shelyag, S., Zharkov, S., Fedun, V., Erdélyi, R., & Thompson, M. J. 2009, A&A, 501, 735
  • Shelyag et al. (2010) Shelyag, S., Mathioudakis, M., Keenan, F. P., & Jess, D. B. 2010, A&A, 515, A107
  • Shelyag & Przybylski (2014) Shelyag, S., & Przybylski, D. 2014, PASJ, 66, S9
  • Socas-Navarro et al. (2014) Socas-Navarro, H., de la Cruz Rodriguez, J., Asensio Ramos, A., Trujillo Bueno, J., & Ruiz Cobo, B. 2014, arXiv:1408.6101
  • Solanki (1987) Solanki, S. K. 1987, Ph.D. Thesis,
  • Unno (1956) Unno, W. 1956, PASJ, 8, 108
  • Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635
  • Zharkov et al. (2013) Zharkov, S., Shelyag, S., Fedun, V., Erdélyi, R., & Thompson, M. J. 2013, Annales Geophysicae, 31, 1357
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description