Investigation of Surface Effects of Simple Flux Tubes Using Numerical Simulations

Sunspots play an important role in understanding the dynamical nature of the solar magnetic field. Although their surface appearance has been observed for over four centuries, little is known about the subsurface structure. They are known to strongly influence solar acoustic modes and there is a variety of possible interactions of the magnetic field with waves in its proximity. Sunspot seismology is the study of these waves and their interactions with the goal of understanding more about subsurface properties, such as magnetic field configurations, mass flows, and thermal structures.

Earlier studies revealed that most techniques of helioseismology break down in the presence of strong magnetic fields (Gizon et al. 2009; Moradi et al. 2010) and can therefore not be used in sunspot seismology. Since there are additional issues with velocity measurements within sunspots, such as the change of height due to the Wilson Depression, atomic lines being affected by magnetic fields and the suppression of oscillations (Braun et al. 1990), the analysis of data can be difficult. Oftentimes the study of isolated and simplified effects therefore rely on MHD simulations. Efforts have been made to find discernible surface signatures caused by flux tubes with different subsurface structures. Examples include investigations of acoustic Halos around sunspots (Rajaguru et al. 2013; Rijs et al. 2016) and scattered wave fields after interaction with magnetic fields (Zhao et al. 2011). Schunker et al. (2013) studied changes to observable travel times caused by abnormalities in the subsurface structure of a flux tube.

An important physical quantity to consider when wave fields are studied is the c_A = c_s layer where the sound speed c_s is equal to the Alfvén speed c_A (Rosenthal et al. 2002; Cally 2007). At this layer mode conversion is most prominent. Mode conversion is fundamental in understanding wave behavior in and around sunspots. Thus the (well established) p-mode absorption, for example, close to active regions (Braun et al. 1987, 1988) could eventually be contributed to conversion of waves in slow mode waves traveling along the flux tube (Cally & Bogdan 1997). Also, most interpretations nowadays of the formation of acoustic halos around active regions include mode conversion in some layer within the solar atmosphere (Hanasoge 2008; Khomenko & Cally 2012; Nutto et al. 2012).

When considering MHD simulations a number of simplifications and limitations need to be applied, so the computational expense remains reasonable. Especially when only atmospheric effects on wave fields are studied, it is justifiable to neglect radiative transfer and therefore convection and granulation. An appropriate code for solving the MHD equations for seismic propagation is the SPARC code, developed by Hanasoge (2007) and Hanasoge & Duvall (2007). It has been used extensively in the past. Rijs et al. (2016) used it to analyze the effect of the Alfvén limiter on the formation of the before mentioned acoustic halos. Shelyag et al. (2009), Przybylski et al. (2015), and Rijs et al. (2015) carried out single source excitation simulations to study how acoustic power is (distributed) along and around flux tubes.

In this work we first visualize the aforementioned slow mode waves traveling downwards into the interior of the simulation's domain. This downward propagation will eventually turn acoustic waves back up to the surface, where they can be measured. Assuming that these waves carry information, among other things, about the subsurface extent of the c_A = c_s layer, a rough image of magnetic field configuration with height can be obtained.

The simulation setup is described in Section 2, including the artificial atmosphere (Section 2.1) plus the flux tube model (Section 2.2). In Section 3 we analyze surface effects of the flux tube in the very simple scenario of one isolated source. Concerning more realistic simulations, we show in Section 4 how to potentially detect such surface effects in real data.

SPARC is a code that can be used to compute the interactions of waves with magnetic flux tubes, sound speed, and damping perturbations, and study the wave field in the presence of multiple/single sources or anomalies thereof. The linearized MHD and Euler equations in 3D Cartesian geometry are solved. The derivatives are computed using sixth-order compact finite differences (in all three directions) or fast Fourier transforms in the horizontal directions and an optimized second-order RK time stepping scheme is implemented (Hanasoge 2007).

As the goal of this work is to find information about subsurface structures at the surface carried by propagating waves, a solar-like background stratification including a magnetic field described by flux tube model needs to be set up.

Our simulation domain is enclosed in a box with 256 × 256 × 300 grid points. The dimensions are 373.76 × 373.76 × 40 Mm³. x and y dimensions are required to be chosen such that typical (solar-like) acoustic wavelengths are resolved. Thus, they were set to be comparable to available data sets of solar surface velocities, such as HMI (link to HMI). For the maximum depth, one has to usually find a trade-off between resolution (i.e., computational expense) and the existence of modes with deep turning point within the simulation domain. 40 Mm is hereby a reasonable choice.

When dealing with isolated sources as in Section 3, the boundaries in each direction included a perfectly matched layer to absorb outgoing waves as efficiently (with as little reflection) as possible (Hanasoge et al. 2010). This is reasonable for this particular scenario, since we consider single wave packets crossing each grid point only once. For the stochastic excitations considered in Section 4 the horizontal boundaries were chosen to be periodic, in order to keep the simulation as realistic as possible.

For every simulation run, the vertical velocities v_z at the surface (z = 0 Mm) are written out and stored at a cadence of Δt = 45 s (again, in order to resolve solar-like acoustic frequencies and to be similar to HMI data sets). In the following, vertical surface velocities stemming from simulations with the fully magnetized atmosphere are labeled as v_magn. In addition wave fields for quiet runs yielding v_quiet were obtained. Quiet in this case means an atmosphere without any magnetic field (effectively 1D). The idea is that the difference of the velocities v_diff = v_magn − v_quiet is basically a noise subtraction (where v_quiet is identified as noise), highlighting remnants of the full wave field influenced by the presence of the magnetic field. These remnants are mostly waves caused by gradients in pressure, density, and sound speed due to modifications of the background (by the magnetic flux tube), but also waves originating from mode conversion. This difference signal v_diff does not represent the scattered wave and might not have a trivial physical meaning, but it still lets us investigate the mode conversion within the flux tube (and any wave field that results from it) much easier than taking only the full wave field v_magn into account. Observation of v_diff in real data is generally not possible; however, in Section 4 we present a method of measuring it, to some extent, indirectly.

Additional simulations are set up, containing only the thermal perturbations of the flux tube within the atmosphere, but not the magnetic field itself (thus yielding v_therm). We use these thermal only runs as a qualitative measure for the direct influence of the magnetic field, similar to Rijs et al. (2015). We refrain from using the difference ${v}_{\mathrm{diff}}^{{\prime} }={v}_{\mathrm{magn}}-{v}_{\mathrm{therm}}$ for the quantitative analysis, since there is no trivial relation between ${v}_{\mathrm{diff}}^{{\prime} }$ and v_diff.

2.1. Stabilized Background Atmosphere

2.2. Monolithic Self-similar Flux Tube Model

For our flux tube model we use a monolithic, self-similar description as presented in Schlüter & Temesváry (1958) and Hanasoge (2008). The basic structure can be seen in Figure 1. The demonstrated magnetic field distribution is further denoted as model 1. To deal with the quickly growing Alfvén speed c_A in the z > 0 layers (and especially within the flux tube), we introduced an Alfvén limiter c_A^max = 90 km s⁻¹. Alfvén limiters should be set as high as possible, since they degrade the realism of the simulation further. Large values, however, make the computational expense large. ${c}_{{\rm{A}}}^{\max }=90$ km s⁻¹ is shown to be a good compromise in Rijs et al. (2016), and allows us to set the simulation time step to Δ_ts = 0.2 s. This work focuses more on the subsurface wave interactions, thus having a large ${c}_{{\rm{A}}}^{\max }$ is not crucial. Δ_ts = 0.2 s allows reasonable simulation wall times.

Zoom In Zoom Out Reset image size

With the method described in Schlüter & Temesváry (1958), the c_A = c_s layer surfaces at r(z = 0) = 19.79 Mm away from the center of the flux tube. This is similar to a rather large, but still realistic sunspot. For the purpose of this work, two more (not necessarily physical) flux tube models were constructed, labeled Model 2 and Model 3. Model 2, shown in Figure 2, left panels, is inspired by Schunker et al. (2013), including a broadening of the flux tube at a depth of z = −2 Mm. A broadening like this might be a consequence of convective motions, fanning out the field lines. Model 3, see Figure 2, right panels, was infused with a sudden increase of the vertical magnetic field B_z at z = −1.5 Mm. This results in a stronger depression of the c_A = c_s layer, affecting the subsurface mode conversion, while still preserving the wine-glass structure. These toy models exhibit approximately the same surface parameters (i.e., B_z(z = 0) ≈ 2.8 kG, r(z = 0) = 19.79 Mm, etc.), which is a requirement in order to see if surface effects are affected only by subsurface properties.

Zoom In Zoom Out Reset image size

2.3. Mode Conversion and Downward Propagation

In order to keep the interactions at the flux tube boundary (that is c_A = c_s) as simple as possible, we employ a method similar to that of Shelyag et al. (2009). As mentioned before, only a single, quickly decaying oscillatory background displacement is used. It is quantified as:

$$\begin{eqnarray*}&&{v}_{z}=\sin \left(\displaystyle \frac{2\pi t}{{p}_{t}}\right)\exp \left(-\displaystyle \frac{{\left(t-{t}_{0}\right)}^{2}}{{\sigma }_{t}^{2}}\right)\exp \left(-\displaystyle \frac{| ({\boldsymbol{x}}-{{\boldsymbol{x}}}_{0}){| }^{2}}{| {{\boldsymbol{\sigma }}}_{x}{| }^{2}}\right)\end{eqnarray*}$$

where p_t is the oscillation period, t₀ is the starting time, σ_t is the temporal width (i.e., length), ${{\boldsymbol{x}}}_{0}$ is the location, and ${{\boldsymbol{\sigma }}}_{x}$ is the spatial width. We set p_t = 302 s, equating to ≈3.31 mHz. Also t₀ = 200 s, σ_t = 75 s (therefore quickly decaying), ${{\boldsymbol{x}}}_{0}\,=\left(-0.4,150.0,186.9\right)$ Mm (origin of the x–y-axis at the corner of the box), and ${{\boldsymbol{\sigma }}}_{x}=\left(0.6,1.5,1.5\right)$ Mm.

This displacement will cause the propagation of a wave packet in every direction, simulated by SPARC. In Figure 3 this propagation is visualized. The left panels show snapshots of the full wave field v_magn where the upper panels show the y-component, and the lower panels the z-component. The simulation of the wave field seen in the left-hand lower panel corresponds to the (theoretical) propagation of a single wave packet on the solar surface. For the chosen time t = 30 minutes, the separation of multiple skip branches (see Section 3) already becomes visible. Note that there are some artificial features, such as absorption at the boundaries and amplitude distortion due to the Cartesian geometry, which is weak enough to be negligible. In the right panels, the instantaneous difference v_diff is shown.

Zoom In Zoom Out Reset image size

Waves emerging from the source will behave as fast (fully) acoustic mode waves, once they cross the c_A = c_s layer, they can convert to the slow magnetoacoustic mode branch (Cally 2007). Slow mode waves will show two properties here: propagation preferably along magnetic field lines, and their transverse nature. Hence, showing the horizontal component will make these waves visible. As denoted in the top right panel of Figure 3 by black arrows, slow mode waves will start to trickle downwards along the flux tube, as long as they are within the c_A = c_s layer (and upwards out of the simulation domain). Once they cross that layer again, they can convert back to fast acoustic mode waves. These will now start to return to the surface after crossing their respective inner turning point, which is visualized in Figure 4 for t = 60 minutes and can best be seen for z < −5 Mm in the top right panel. Finally, the returning waves form a circular wave pattern centered on the flux tube, as seen in the bottom right panel of Figure 4.

Zoom In Zoom Out Reset image size

To summarize: waves that cross the flux tube boundary get partially converted and travel downwards, along field lines. They convert back again and travel to the surface. This is especially interesting, since the subsurface properties of the flux tube (at least within the c_A = c_s layer) will have an influence on the wave field that can be observed at the surface. Indeed we expect that the travel time, due to time being spent subsurface, and the shape, due to interference on the surface of the flux tube boundary, of the re-emerging wave packet will change. This change is thus related to the topology of the c_A = c_s layer. A quantitative analysis is done in the next section.

The surface signal v_diff as seen in the bottom right of Figure 5 is rather weak (about 10% of the full wave field amplitude), but is still contained in v_magn (as v_magn = v_diff + v_quiet). It is however superimposed by the wave packet of the original source, making it difficult to detect. In Sections 3 and 4 it is shown how the two can be separated and potentially measured in real data.

Zoom In Zoom Out Reset image size

Using the difference v_diff' we observe a similar wave pattern, but with an amplitude of about 1% of v_magn. From this behavior we learn that the re-emerging wave packet that is observed at the surface cannot only stem from thermal modifications to the background atmosphere (due to the presence of a magnetic flux tube), but must also carry contributions from the subsurface magnetic field itself. This is an important insight, meaning that v_diff is in fact sensitive to subsurface magnetic fields.

The method we use to describe the wave propagation quantitatively is the time–distance analysis (Duvall et al. 1993; Gizon & Birch 2005). By cross-correlating the point x₁ in which the single source is located, with an arc ${\bar{x}}_{2}({\rm{\Delta }})$ at different times t and distances Δ, we can measure quantities like the group travel time t_g, the amplitude A, and the central frequency ν₀ of the propagating wave packet. For this analysis the source is now put in the center of the simulation box ${{\boldsymbol{x}}}_{0}=(-0.4,186.88,186.88)$ Mm, and thus x₁ is put directly within the flux tube. This positioning eliminates the initial travel time delay, meaning that the wave packet will instantly interact with the magnetic field, making the analysis more simple. The cross-correlation $C(t,{\rm{\Delta }})$ calculated as a function of Δ is then shown in Figure 5. For the two simulation runs v_magn and v_quiet we define:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\mathrm{magn}}(t,{\rm{\Delta }}) & \equiv & \mathrm{Cross} \mbox{-} \mathrm{correlation}\ \mathrm{for}\ \mathrm{full}\ \mathrm{magnetic}\ \mathrm{run}\\ {C}_{\mathrm{quiet}}(t,{\rm{\Delta }}) & \equiv & \mathrm{Cross} \mbox{-} \mathrm{correlation}\ \mathrm{without}\ \mathrm{magnetic}\ \mathrm{field}\\ {C}_{\mathrm{diff}}(t,{\rm{\Delta }}) & \equiv & \mathrm{Cross} \mbox{-} \mathrm{correlation}\ \mathrm{of}\ \mathrm{the}\ \mathrm{difference}\end{array}\end{eqnarray} \tag{ 1 }$$

where by difference, again the instantaneous difference ${v}_{\mathrm{diff}}\,={v}_{\mathrm{magn}}-{v}_{\mathrm{quiet}}$ is meant. Note that the initial point x₁ for the cross-correlation C_diff is taken from the full wave field v_magn, since we want to analyze the correlation of the waves contained in v_quiet with those of the initial displacement. Since the amplitude of v_diff is only about 10% of that of v_magn, it is expected that this is also the case for the amplitude of C_diff, as can be seen in the right panel of Figure 5.

The broad contributions that can be seen for t > 90 minutes and Δ < 210 Mm can be assigned to a (unphysical but weak) reflection at the bottom boundary. Since these are restricted to an area that is generally not of interest, they can be ignored. Furthermore, the usual multiple skip branches are seen very clearly in both cases.

Again, we repeat this analysis with the subtraction of thermal runs, i.e., v_therm', in which we replaced the difference signal accordingly to calculate the correlation C_diff. As expected, the main distinction is the amplitude. Due to the similarity of both correlation functions, we qualitatively conclude that C_diff is in fact sensitive to the magnetic field.

As described in Section 2.3, the wave packet caused by the oscillatory background displacement will travel a certain distance (depending on the depth of the inner turning point) before re-emerging at the surface. It is therefore expected that the group travel times t_g,magn of C_magn and t_g,diff of C_diff differ. This can be measured by employing a fit to the data, to quantify the group travel times separately. For this, Gabor wavelets (Gizon & Birch 2005) of the form

$$\begin{eqnarray}&&C(t)=A\exp \left(-{\left(2\pi d\nu \right)}^{2}{\left(t-{t}_{{\rm{g}}}\right)}^{2}\right)\cos \left(2\pi {\nu }_{0}(t-{t}_{{\rm{p}}}\right)\end{eqnarray} \tag{ 2 }$$

are used. With A being the amplitude, dν the width and ν₀ the central frequency of the spectral distribution of the wavelet, t_g the group speed, and t_p the phase speed. The results of ${t}_{{\rm{g}},\mathrm{magn}}$ and ${t}_{{\rm{g}},\mathrm{diff}}$ for this fit are shown for the first three branches in the left panel of Figure 6. A shift ${\rm{\Delta }}{t}_{{\rm{g}}}={t}_{{\rm{g}},\mathrm{magn}}-{t}_{{\rm{g}},\mathrm{diff}}$ can clearly be seen for all branches and is highlighted with an annotation for the one-skip branch. Fitting the branches becomes more difficult starting from three-skip, since for such distances, multiple skips start to overlap. Since ${\rm{\Delta }}{t}_{{\rm{g}},1}$ only varies slightly with the distance Δ, we estimate it by taking the average over Δ. This yields:

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{{\rm{g}},1}=282.6\,{\rm{s}}.\end{eqnarray} \tag{ 3 }$$

The same analysis for ${\rm{\Delta }}{t}_{{\rm{g}},2}$ and ${\rm{\Delta }}{t}_{{\rm{g}},3}$ gives similar results, but is less precise, since the fit does not converge as consistently and is generally less robust. We conclude that ${\rm{\Delta }}{t}_{{\rm{g}}}\equiv {\rm{\Delta }}{t}_{{\rm{g}},1}$ is the most accurate estimation of the time delay.

Zoom In Zoom Out Reset image size

Shown in the center panels is the phase travel time t_p. Here a similar behavior is observed, although the time delay

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{{\rm{p}},1}=311.1\ {\rm{s}}\end{eqnarray} \tag{ 4 }$$

is more than twice as large as t_g,1. However, the fit is much more erratic and becomes less reliable for small and large distances. For the analysis presented in the next section, we will be using t_g,1 only, in principle it can nevertheless be performed with any quantity yielded by the fit.

Since one of the goals of this work is learning about the subsurface configuration of the simulated flux tube, this study is repeated for the other two flux tube models (see Figure 2). It turns out that changing the model affects Δt_g only weakly. However, as shown in the right panel of Figure 6 other parameters like the central frequency ν₀ (see Equation (2)) of the wavelet will show significant changes Δν₀. While small changes in the model (model 1 to model 2) almost do not affect ν₀ (Δν₀ ≈ 0.04 mHz), more drastic changes (model 1 to model 3) may increase the frequency by up to Δν₀ ≈ 0.5 mHz. This will allow us to distinguish different subsurface magnetic field configurations, as long as the signal C_diff is available, or can be reconstructed reasonably well (more on this in Section 4).

The fact that the group travel time delay Δt_g is similar for all three models can be explained by considering the radius of the c_A = c_s layer. Wave emergence visible in v_diff will start at the boundary of said layer, leading to similar subsurface travel times for similar radii r(z). When slow mode waves are converted back to fast mode acoustic waves at the boundary, interference between them happens, depending on the shape of the c_A = c_s layer. This interference will influence the wave pattern observed at the surface. That explains why ν₀(Δ) (being an indicator for the spatial frequency distribution of the wave packet) is very similar for model 2 and the original model 1, but shows significant changes of up to 0.5 mHz for model 3, when the c_A = c_s-layer shape is changed drastically.

Using real data, the instantaneous difference v_diff is generally not available. Any surface waves, that were caused by the mode conversion mechanism as described in Section 2.3 (seen in the right panels of Figure 4) are expected to be a feature of real sunspots as well. They will, however, have a small amplitude compared to other acoustic signals and will be swamped by them. In this section we will show one possible way to achieve a separation of the desired signal from the background nonetheless.

For this we ran two additional simulations with a solar-like stochastic forcing function, similar to the one described in Hanasoge (2007) again with and without magnetic field. Since the simulation does not include radiative transfer and the background is artificially stabilized, the resulting velocity fields will not exhibit granulation. Thus the simulation loses some degree of realism. For most acoustic analyses, however, granulation noise is unwanted anyway. The results are therefore treated as if granulation has been removed beforehand (by, for example, averaging or filtering processes). Looking again at the instantaneous difference v_diff, we see a circularly shaped wave pattern emerging from around the flux tube, equivalently to the behavior observed in Figure 4. Note that this is also briefly described in Hanasoge (2007).

We employ a time–distance analysis, similar to what we did in Section 3, considering the fact that the cross-correlation signal originating from within the flux tube C_diff is delayed by the time Δt_g. This time, however, the stochastic sources add a lot of unwanted contributions and noise. To still be able to construct a time–distance diagram, we employ a filtering process described in Gizon & Birch (2005) and rely on a point-to-arc average. A phase-speed filter that prefers waves with an average (one skip) travel distance of 14.5 Mm was applied. It is preferable to consider waves emerging close to the flux tube, since the amplitude of v_diff is higher, due to its circular nature. Furthermore, we use points in a circle with a radius of r_C = 29.31 Mm around the flux tube (instead of directly inside it) to be correlated. Altogether, 616 point-to-arc cross-correlation functions are averaged.

As mentioned, the flux tube signal will be weak and superimposed by the cross-correlation signal of the original wave packet C_magn, making it difficult to fit the desired signal C_diff directly. The superposition will, however, modify the wave packet in an asymmetric way, since it only contributes at specific times (t + Δt_g, no contribution at t − Δt_g). This can be quantified by using the Gizon & Birch (2004) fitting method. It requires essentially only a reference function C⁰(t,Δ) that predicts the expected travel time at a certain distance Δ. In our case this function is easily obtained by using single source simulations as in Section 3. For real data one can either use theoretical predictions (i.e., the simulations done for this work, see Figure 5), or a cross-correlation function from heavily averaged time–distance diagrams. The result of this fitting method is given in Figure 7 at an exemplary distance Δ = 13.2 Mm, where the travel time t_g,magn is annotated within the according image. A window with a width of ${\sigma }_{\mathrm{window}}\,=13.5\,\mathrm{minutes}$ (shown as dashed orange line) was applied the interval of estimated travel time.

Zoom In Zoom Out Reset image size

The challenge is now to detect a change in the travel time due to the emerging wave signal of the flux tube. We do this by "correcting" the data ${C}_{\mathrm{magn}}(t,{\rm{\Delta }})$, and the reference function ${C}^{0}(t,{\rm{\Delta }})$. In principle the time–distance diagram on the right in Figure 7 exhibits contributions as shown on the right in Figure 5. This can be modeled, as long as we can estimate ${t}_{{\rm{g}},\mathrm{diff}}$ and the amplitude A_diff of the cross-correlation C_diff with t_g,Model and A_Model. The actual cross-correlation functions for the full wave field C_magn and instantaneous difference C_diff differ in more than just the amplitude and travel time (see Figure 6, right panel), which is, however, neglected for the simplicity of this analysis. For the correction, we simply subtract the contribution of the flux tube:

$$\begin{eqnarray}&&{C}_{\mathrm{magn}}^{{\prime} }={C}_{\mathrm{magn}}-{C}_{\mathrm{Model}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{C}^{0^{\prime} }={C}^{0}-{C}_{\mathrm{Model}}\end{eqnarray} \tag{ 6 }$$

where primed quantities means corrected. Correcting the data like this assumes that the difference signal v_diff contributes to the correlation C_magn in a linear fashion. Generally the construction of C_magn from v_magn and v_diff is more complicated, therefore the correction in Equation (6) is based on a first-order approximation. C_Model is then constructed via

$$\begin{eqnarray}&&{C}_{\mathrm{Model}}(t,{\rm{\Delta }})={A}_{\mathrm{Model}}\cdot {C}^{0}\left(t+{\rm{\Delta }}{t}_{{\rm{g}}},{\rm{\Delta }}\right),\end{eqnarray} \tag{ 7 }$$

with, in this case, ${\rm{\Delta }}{t}_{{\rm{g}}}={t}_{{\rm{g}},\mathrm{Model}}-{t}_{{\rm{g}}}^{0}$. Fitting the corrected data ${C}_{\mathrm{magn}}^{{\prime} }$ will then yield a slightly modified travel time ${t}_{{\rm{g}}}^{{\prime} }$, shown in Figure 8. The difference for the two results ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }$ is not equal to ${\rm{\Delta }}{t}_{{\rm{g}}}$, as it is only a slight modification due to the contribution of the flux tube signal C_diff. It is, however, related to Δt_g, by how well the model attempt C_Model agrees with C_diff. An exemplary attempt at constructing the model C_Model is shown in Figure 9. Generally it is expected that ${t}_{{\rm{g}}}^{{\prime} }\lt {t}_{{\rm{g}}}$, since ${t}_{{\rm{g}},\mathrm{magn}}\,\gt {t}_{{\rm{g}},\mathrm{quiet}}$. Also, if C_Model does not agree well with C_diff, it is expected that ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }\lt {\rm{\Delta }}{t}_{{\rm{g}}}$. Summarizing, we can make the first constraint on the measurement of t_g', and therefore for the agreement between C_Model and C_diff:

$$\begin{eqnarray}&&{t}_{{\rm{g}}}-{\rm{\Delta }}{t}_{{\rm{g}}}\lt {t}_{{\rm{g}}}^{{\prime} }\lt {t}_{{\rm{g}}}.\end{eqnarray} \tag{ 8 }$$

In our exemplary analysis, we find ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }=159.7\,{\rm{s}}$, where we set (${t}_{{\rm{g}},\mathrm{Model}}-{t}_{{\rm{g}}}^{0}=$) ${\rm{\Delta }}{t}_{g}=270.0\,{\rm{s}}$ and ${A}_{\mathrm{Model}}=0.20$. Since for simulations, C_diff is available, we can calculate the theoretical value of ${t}_{{\rm{g}}}^{{\prime} }$ and deduce Δt_g from fits, as shown in Figure 6, for comparison. Here we find:

$$\begin{eqnarray}&&{t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }=191.7\ {\rm{s}}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{g}=282.6\ {\rm{s}}.\end{eqnarray} \tag{ 10 }$$

As can be seen, the initial estimate of ${\rm{\Delta }}{t}_{g}=270.0\,{\rm{s}}$ leads to a ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }$ that is already close to the theoretical value.

Zoom In Zoom Out Reset image size

The constraint (8) yields a broad estimate on how to choose the parameters A_Model and t_g,Model for the model attempt C_Model, but the exact value of Δt_g will remain unknown for real data. One would need to do this analysis for several sunspots, to further narrow down Equation (8), or rely on the simulations done in this work, to have a reference for the required estimate of Δt_g.

Zoom In Zoom Out Reset image size

The mechanism of acoustic fast mode waves being converted into downward propagating slow mode waves is a known process (Cally et al. 2003; Rijs et al. 2015, 2016), but it has never been studied in particular. Moreover the ramping effect (Cally 2007), which is similar in its nature (upward propagation instead of the downward propagation considered here), was investigated more extensively, due to its possible contribution to acoustic halos. In this work, it was shown that flux tubes behave as sources of acoustic power, as long as they are being excited from the outside (Section 2.3). The consequence is that emerging waves from a sunspot contribute to the acoustic power in its vicinity. Of course, the absorption of p-mode power (Braun et al. 1987, 1988) makes detecting this power excess in real data nontrivial.

The time delay Δt_g shown in Figure 6 is a way to distinguish the signals C_magn and C_diff. Setting up simulations with solar-like sources, it was shown that although the amplitude of the C_diff is low, Δt_g can still be estimated. As shown in Section 3 (left panel of Figure 6), a decent reconstruction of C_diff via C_Model using an estimate for Δt_g may reveal subsurface properties of the investigated flux tube. It is necessary, however, to obtain the wave pattern of C_diff to get dν and ν₀ (see Equation (2)). Since t_g does not vary for different flux tube models, dν and ν₀ are needed to relate the surface signal to subsurface properties. It also appears that slight changes to the model that do not affect the c_A = c_s-layer do not affect the surface signal C_diff. For real data this will require additional effort in creating C_Model.

With the proposed method in Section 4, a measurement of ${t}_{{\rm{g}}}^{{\prime} }$ needs to be done reliably, which will be difficult in the case of real data. In principle, the method can be done for more sets of filters, as described in Gizon & Birch (2005) and thus, more distances Δ. We tested this here, and got similar results, but not for all distances Δ. This is expected, since, as mentioned, the amplitude of C_diff becomes weaker for large Δ, decreasing the signal-to-noise ratio, due to its circular wave behavior. Again, another hurdle in estimating Δt_g reliably.

Using the SPARC code to simulate the interaction of different kinds of waves with simple flux tube models, the effects of mode conversion have been visualized (see Figure 3). Slow mode waves traveling downwards along magnetic field lines of the flux tube convert back to acoustic fast mode waves, that deflect back up and are measurable at the surface. It is demonstrated how these surface waves are altered from subsurface changes in the flux tube model (see Figure 5). Moreover, the fact that these waves spend time traveling within the flux tube, a time delay Δt_g (see Equation (3)) between re-emerging wave and original (as caused by the initial source) wave can be measured. We find that ${\rm{\Delta }}{t}_{{\rm{g}}}=282.6\,{\rm{s}}$. It was also found that subsurface changes in the flux tube models, especially changes to the shape of the c_A = c_s-layer, influence the frequency distribution of the surface wave pattern (see Figure 6).

A method to estimate Δt_g for real data is presented in Section 4. Although the method might become unreliable for real data due to many sources of noise, for our simulations a value of ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }=159.7\,{\rm{s}}$ with an assumption of Δt_g = 225 s proved to be accurate when compared to the theoretically predicted value of ${t}_{{\rm{g}}}-{t}_{{\rm{g}}}^{{\prime} }=191.7\ {\rm{s}}$ and the associated ${\rm{\Delta }}{t}_{g}\,=282.6\ {\rm{s}}$.

The reconstruction of C_Model with quantities available in real data, in order to estimate C_diff will need some additional effort, to be reliable. Also, measuring properties of the wave pattern, like the central frequency ν₀, might require direct detection of C_diff, which again will be difficult due to its comparatively low amplitude.

This work serves as a theoretical basis for a new method with the potential of adding to the knowledge of subsurface sunspot properties. The next step for further analysis regarding this topic is executing this study for real data and making these simulations more realistic by, for example, tuning the flux tube model, the background model etc. This will include fine tuning the proposed method as, for example, in Equation (8).

We thank Shravan M. Hanasoge for making the SPARC code, being the basis for this work, publicly available at http://www2.mps.mpg.de/projects/seismo/sparc/.