SCATTERING MATRIX FOR THE INTERACTION BETWEEN SOLAR ACOUSTIC WAVES AND SUNSPOTS. I. MEASUREMENTS

One of the major subjects in helioseismology is the use of solar acoustic waves to probe sunspots. Several data analysis methods have been developed for this purpose, such as the ring-diagram analysis (Hill 1988; Rajaguru et al. 2001; Nicholas et al. 2004; Bogart et al. 2008; Baldner et al. 2013; Jain et al. 2015), Hankel analysis (Braun et al. 1987; Bogdan et al. 1993; Braun 1995; Chen et al. 1996; Sun et al. 1997; Couvidat 2014), time–distance analysis (Duvall et al. 1993, 1996; Kosovichev 1996; Chou et al. 2000, 2009a, 2009b, 2009c; Kosovichev et al. 2000; Zhao et al. 2010, 2011b), and acoustic imaging (or acoustic holography; Chang et al. 1997; Chen et al. 1998; Braun & Lindsey 1999, 2000; Chou et al. 1999; Chou 2000; Lindsey & Braun 2000; Sun et al. 2002).

A different method has been developed recently. Cameron et al. (2008) calculated the cross correlation between a line-averaged Doppler signal outside a sunspot and the signal at the points in and around the sunspot for f-mode waves. The cross-correlation function, as a function of space and time, represents the propagation of an acoustic wave packet starting from the line through the sunspot. However, the cross-correlation function is not the wavefunction of the radial velocity of solar acoustic waves. The wavefunction of radial velocity is the radial component of velocity vector, which is the solution of the wave equation for solar acoustic waves. Zhao et al. (2011a) applied a deconvolution scheme on the cross-correlation function to remove the Fourier amplitude of the signal at the reference line, where the incident wave starts, and keep only its phase in the the cross-correlation function in order to obtain the wavefunctions of radial velocity of solar acoustic waves scattered by a sunspot for each radial order. Zhao & Chou (2013) further improved the method and were able to measure the scattered wavefunctions between adjacent radial orders. The scattered wavefunction is defined as follows: For an incident wave propagating toward a sunspot, after interacting with the sunspot, the wave is modified around the sunspot. The modified wave can be expressed as the sum of the incident wave and the scattered wave. The scattered wave satisfies an inhomogeneous wave equation. The inhomogeneous term, called the source term, is responsible for the generation of the scattered wave (Zhao et al. 2011a). Yang et al. (2012) set up a phenomenological model to describe this source term, and used the measured scattered wavefunctions to determine the parameters of the model by iterative fits up to a high-order Born approximation. Chou et al. (2012) used the measured total wavefunction (the sum of the incident and scattered wavefunctions) around the sunspot to compute a one-dimensional hologram, and demonstrated that this hologram could be used to reconstruct the two-dimensional wavefronts on the surface. Zhao & Chou (2016) used the measured scattered wavefunctions to compute the absorption cross section and scattering cross section for each radial order for the interaction between waves and sunspots.

Assessing the interaction between solar acoustic waves and sunspots can be considered a scattering problem. For a scattering problem, the most complete measured information is the scattering matrix. Each element of the scattering matrix is a complex number, representing the transition between two modes, or the conversion from one incident mode to another mode due to the interaction. The discussion of the general theory of the scattering matrix can be found in many physics textbooks (Bjorken & Drell 1964; Schiff 1971). For the interaction of acoustic waves with sunspots, the theory of the scattering matrix has been studied by various authors (Bogdan & Cattaneo 1989; Chou 1994; Fan et al. 1995; Chou et al. 1996; Hindman & Jain 2012). So far no correct measurement of the scattering matrix for the interaction between acoustic waves and sunspots has been reported. Hankel analysis has measured the amplitude change and phase change between an ingoing mode and an outgoing mode. However, these changes are not the amplitude and phase of the individual scattering matrix element. Instead, they are certain combinations of scattering matrix elements, as expressed by Equations (26) and (28) in Braun (1995). This is because the incident wave in Hankel analysis consists of many modes instead of a single mode. The multimode of the incident wave is such that each outgoing mode is a mixture of the scattered signals coming from various ingoing modes, owing to the mode mixing (see also Hindman & Jain 2012). If the mode mixing is negligibly small, the measurement in Hankel analysis can serve as an estimate for the scattering matrix elements. However, the mode mixing is non-negligible and even significant for some modes.

In this study, we devise a method of using a single-mode incident wave to measure the scattering matrix elements. The scattered wavefunctions measured in Zhao & Chou (2013) are the signals scattered from the incident wave of one radial order into the wave of another radial order, because these wavefunctions were measured with the cross correlation between the waves of two radial orders. Since the incident wave belongs to a particular radial order and propagates in a particular direction, each of its temporal Fourier components is a single mode (plane wave). For a time-independent sunspot, there is no mode mixing between different frequencies. Thus each temporal Fourier component of the scattered wave comes from the single-mode incident wave of the same frequency. Expanding the scattered wave of this frequency in terms of different spatial modes gives the scattering matrix elements between this incident mode and various scattered modes. In this study, we use the scattered wavefunctions measured in Zhao & Chou (2013) to compute the amplitude and phase of the scattering matrix elements for the transitions above the noise level.

In Section 2, the general definitions of transition matrix and scattering matrix for a single-mode incident wave are discussed. In Section 3, the method of computing transition matrix elements using the scattered wavefunctions measured in Zhao & Chou (2013) is discussed. Section 4 discusses the methods of determining the amplitude and phase of transition matrix elements. In Section 5, each step of data analysis in computing the transition matrix element from the measured scattered wavefunctions is described and the determination of the amplitude and phase of the measured transition matrix elements is discussed. Section 6 discusses the result of the amplitude and phase of the detectable transition matrix elements. The Appendix discusses the determination for the phase and amplitude of a single-mode signal if its wavenumber is not exactly equal to one of the Fourier wavenumbers.

Consider the interaction of the acoustic waves with a sunspot as a scattering problem. An incident wave, represented by the wavefunction of velocity vector ${{\boldsymbol{V}}}^{({\rm{i}})}({\boldsymbol{r}},t)$, propagates toward the sunspot, where the superscript i denotes the incident wave. It satisfies the wave equation of the quiet Sun. After interacting with the sunspot, the incident wavefunction ${{\boldsymbol{V}}}^{({\rm{i}})}({\boldsymbol{r}},t)$ is modified to ${\boldsymbol{V}}({\boldsymbol{r}},t)$ in and around the sunspot. The difference ${{\boldsymbol{V}}}^{({\rm{s}})}({\boldsymbol{r}},{\boldsymbol{t}})\equiv {\boldsymbol{V}}({\boldsymbol{r}},t)-{{\boldsymbol{V}}}^{({\rm{i}})}({\boldsymbol{r}},t)$ is defined as the scattered wavefunction (Morse & Feshbach 1953; Jackson 1999), where the superscript ${\rm{s}}$ denotes the scattered wave. The scattered wavefunction ${{\boldsymbol{V}}}^{({\rm{s}})}({\boldsymbol{r}},t)$ represents the signal generated by the sunspot that is perturbed by the incident wave.

In this section, a general definition for the transition matrix and the scattering matrix is given. First, the formal definition for the transition matrix in the three-dimensional space is discussed. Then its relation to a measurable quantity on the two-dimensional surface, called the "two-dimensional transition matrix," is discussed. The eigenfunctions of velocity (or displacement) vector of solar oscillations in the quiet Sun are orthonormal and form a complete set, denoted as $\{{{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })\}$, where β labels the spatial mode parameters and ${\omega }_{\beta }$ is the corresponding angular frequency. The eigenfunction ${{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$ is normalized as $\int {{\boldsymbol{V}}}_{\beta }^{(0)* }({\boldsymbol{r}},{\omega }_{\beta })$ $\cdot {{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$ $\rho (r){d}^{3}r={\delta }_{\beta \alpha }$, where the weighting function $\rho (r)$ is the mass density (Unno et al. 1989). Since we are working in a small area, the coordinates ${\boldsymbol{r}}$ can be approximated by the Cartesian coordinates $(x,y,z)$, where (x, y) are the horizontal directions and z is the vertical (radial) direction, and β labels $(n,{k}_{x},{k}_{y})$, where n is the radial order, and k_x and k_y are the horizontal wavenumbers. The frequency ${\omega }_{\beta }$ is connected with $(n,{k}_{x},{k}_{y})$ through the dispersion relation.

Suppose the incident wave is a single mode and its wavefunction is an eigenfunction of the quiet Sun ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha }){e}^{-i{\omega }_{\alpha }t}$. If the sunspot is time independent, the interaction with the sunspot would not modify the frequency ${\omega }_{\alpha }$. Since the sunspot has the spatial dependence, the interaction with the sunspot would modify the wavenumber $(n,{k}_{x},{k}_{y})$. The scattered wavefunction corresponding to this incident wave is denoted as ${{\boldsymbol{V}}}_{\alpha }^{\,({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha }){e}^{-i{\omega }_{\alpha }t}$, where the subscript α in ${{\boldsymbol{V}}}_{\alpha }^{\,({\rm{s}})}$ specifies this scattered wave coming from the incident mode α. The spatial part ${{\boldsymbol{V}}}_{\alpha }^{\,({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ is different from the incident mode ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$ and can be expanded in terms of the complete and orthonormal set of $\{{{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })\}$ as

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })=\displaystyle \sum _{\beta }\,{T}_{\beta \alpha }\,{{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta }),\end{eqnarray} \tag{ 1 }$$

where the expansion coefficient ${T}_{\beta \alpha }$, a complex number, is called the element of the transition matrix (T-matrix) with the eigenfunctions $\{{{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })\}$ as the basis (Bjorken & Drell 1964; Schiff 1971). It is noted that the T-matrix element is defined in the three-dimensional space. Since the scattered wavefunction ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ has the frequency ${\omega }_{\alpha }$, the coefficient ${T}_{\beta \alpha }$ is nonzero only when the mode β has the same frequency ${\omega }_{\alpha }$. If the T-operator $\hat{T}$ is defined as ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })\equiv \hat{T}{{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$, the T-matrix element ${T}_{\beta \alpha }$ is calculated as ${T}_{\beta \alpha }=\int {{\boldsymbol{V}}}_{\beta }^{(0)* }({\boldsymbol{r}},{\omega }_{\beta })\cdot \hat{T}{{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })\rho (r){d}^{3}r$. The scattering matrix (S-matrix) is the sum of the unity matrix and the T-matrix.

In helioseismology, the three-dimensional distribution of the vector wavefunction ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ in Equation (1) is unavailable. The observed helioseismic data are the Doppler signals which are the line of sight components of oscillatory velocity on the solar surface. For the modes of interest, the horizontal component is smaller than the radial component on the surface (Unno et al. 1989). Since the data used in this study are near the disk center, the contribution from the radial component is predominant in the observed Doppler signal. The radial component of the velocity is restored by a simple geometric correction, as most of data analysis (e.g., Zhao & Chou 2013). The measured wavefunction used in this study is the radial component of velocity on the solar surface.

If the radial components of ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ and ${{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$ in Equation (1) are denoted as ${v}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ and ${v}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$, respectively, the radial component of Equation (1) is

$$\begin{eqnarray}&&{v}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })=\displaystyle \sum _{\beta }\,{T}_{\beta \alpha }\,{v}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta }).\end{eqnarray} \tag{ 2 }$$

If Equation (2) is evaluated at the surface $z={z}_{0}$, ${v}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ on the lhs becomes ${v}_{\alpha }^{({\rm{s}})}(x,y,{z}_{0},{\omega }_{\alpha })\equiv {{\rm{\Psi }}}_{\alpha }^{({\rm{s}})}(x,y,{\omega }_{\alpha })$, which is the scattered wavefunction measured on the surface; ${v}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$ on the rhs at the surface becomes ${v}_{\beta }^{(0)}(x,y,{z}_{0},{\omega }_{\beta })\,\equiv {{\rm{\Phi }}}_{\beta }({z}_{0}){\phi }_{\beta }(x,y)$, where ${{\rm{\Phi }}}_{\beta }({z}_{0})$ is the radial eigenfunction evaluated at the surface and ${\phi }_{\beta }(x,y)$ is the two-dimensional orthonormal eigenfunction in the quiet Sun. In our study, ${\phi }_{\beta }(x,y)$ is the normalized plane wave on the surface, but it can in general be any two-dimensional orthonormal function. The proportional constant of ${{\rm{\Phi }}}_{\beta }(z)$ is determined by that of ${v}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$, which is determined by the normalization of the velocity vector ${{\boldsymbol{V}}}_{\beta }^{(0)}({\boldsymbol{r}},{\omega }_{\beta })$. With the notation defined previously, Equation (2) evaluated at the surface can be expressed as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\alpha }^{({\rm{s}})}(x,y,{\omega }_{\alpha })=\displaystyle \sum _{\beta }\,{T}_{\beta \alpha }\,{{\rm{\Phi }}}_{\beta }({z}_{0})\,{\phi }_{\beta }(x,y).\end{eqnarray} \tag{ 3 }$$

This equation is an expression on the two-dimensional (x, y) plane. It can be viewed as the scattered wavefunction measured on the surface, ${{\rm{\Psi }}}_{\alpha }^{({\rm{s}})}(x,y,{\omega }_{\alpha })$, expanded in terms of a set of two-dimensional orthonormal eigenfunctions $\{{\phi }_{\beta }(x,y)\}$. The expansion coefficient ${T}_{\beta \alpha }{{\rm{\Phi }}}_{\beta }({z}_{0})$, called the "two-dimensional T-matrix element," is calculated by ${T}_{\beta \alpha }{{\rm{\Phi }}}_{\beta }({z}_{0})\,=$ $\int {\phi }_{\beta }^{* }(x,y){{\rm{\Psi }}}_{\alpha }^{({\rm{s}})}(x,y,{\omega }_{\alpha }){dxdy}$ on the surface. Thus the T-matrix element ${T}_{\beta \alpha }$ can be obtained from dividing the measured two-dimensional T-matrix element by ${{\rm{\Phi }}}_{\beta }({z}_{0})$, which can be calculated with a standard solar model.

The measurement of the scattered wavefunction in Zhao & Chou (2013) does not use a single mode as the incident wave. Instead, the wave packet of a particular radial order n propagating in a particular horizontal direction is used as the incident wave. The scattered wavefunction is measured with the cross correlation between the signal of n and the signal of n', where n' could be equal to or different from n. That is, the measured scattered wavefunction represents the signal converted from the incident wave of n to the scattered wave of n'. Both the incident wave and the scattered wave consist of many modes of different frequencies. Therefore the formulas in Section 2 for the single-mode incident wave cannot be applied directly. Some more processes are needed.

Since the interaction does not change the frequency, part of the energy of an incident mode would be converted into other modes with different wavenumbers but the same frequency. Therefore the scattering problem can be worked in the ω domain for each ω separately. The incident wave of n, denoted as ${{\boldsymbol{V}}}_{n}^{{\text{}}({\rm{i}})}({\boldsymbol{r}},t)$, is a wave packet propagating in a particular direction on the surface. It consists of many modes of different frequencies propagating in the same horizontal direction, and can be expressed in terms of temporal Fourier series as

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{n}^{({\rm{i}})}({\boldsymbol{r}},t)=\displaystyle \sum _{\alpha }{A}_{\alpha }^{n}{{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha }){e}^{-i{\omega }_{\alpha }t},\end{eqnarray} \tag{ 4 }$$

where ${A}_{\alpha }^{n}{{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ is the Fourier coefficient and ${{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ is the eigenfunction of the quiet Sun with a radial order n, propagating in the particular horizontal direction. In our study, the incident wave propagates in the y direction. Thus only the modes with k_x = 0 can be present in each Fourier component ${{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ in Equation (4). Moreover, from the dispersion relation $\omega (n,{k}_{x},{k}_{y})$, each frequency corresponds to only one k_y for each n. Therefore there is only one mode (${\omega }_{\alpha },n,{k}_{x}=0,{k}_{y})$ present in each Fourier component ${{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$.

It is noted that each ${{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ in Equation (4) does not correspond to a real single mode; instead, it corresponds to one bin in the domain of $(\omega ,{k}_{x},{k}_{y})$. The size of the bin depends on the spatial size and temporal duration of the data used in the Fourier expansion. In the following discussion, the term "mode" is used for the bin in the domain of $(\omega ,{k}_{x},{k}_{y})$.

Now we come to the scattered waves. The scattered wavefunction from n to n' is denoted as ${{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},t)$. It can be expanded in terms of a temporal Fourier series as

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},t)=\displaystyle \sum _{\alpha }\,{{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha }){e}^{-i{\omega }_{\alpha }t},\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ is the Fourier component. Since the interaction does not cause any cross talk between the modes of different frequencies, each scattered wave ${{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ in Equation (5) comes from the corresponding incident wave ${A}_{\alpha }^{n}{{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ in Equation (4). Each incident mode ${{\boldsymbol{V}}}_{\alpha }^{(0)n}({\boldsymbol{r}},{\omega }_{\alpha })$ is converted into the scattered wave, denoted as ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})n^{\prime} }({\boldsymbol{r}},{\omega }_{\alpha })$. Thus we have ${{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })={A}_{\alpha }^{n}{{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})n^{\prime} }({\boldsymbol{r}},{\omega }_{\alpha })$. Since the scattered wave ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})n^{\prime} }({\boldsymbol{r}},{\omega }_{\alpha })$ propagates in multiple directions, it contains many modes with frequency ${\omega }_{\alpha }$ and radial order n'. Like Equation (1), the scattered wave ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})n^{\prime} }({\boldsymbol{r}},{\omega }_{\alpha })$ can be expanded in terms of the eigenfunctions of the quiet Sun with radial order n', $\{{{\boldsymbol{V}}}_{\beta }^{(0)n^{\prime} }({\boldsymbol{r}},{\omega }_{\beta })\}$, and the expansion coefficients are the T-matrix elements. Therefore we arrive at

$$\begin{eqnarray}&&{{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })={A}_{\alpha }^{n}\displaystyle \sum _{\beta }\,{T}_{\beta \alpha }^{n^{\prime} n}{{\boldsymbol{V}}}_{\beta }^{(0)n^{\prime} }({\boldsymbol{r}},{\omega }_{\beta }),\end{eqnarray} \tag{ 6 }$$

where ${\omega }_{\beta }={\omega }_{\alpha }$ and the expansion coefficient ${T}_{\beta \alpha }^{n^{\prime} n}$ is the T-matrix element between the incident mode α of radial order n and the scattered mode β of radial order n'. It is worth emphasizing that this expression is valid only for time-independent sunspots.

Like the discussion leading to Equation (3), we take the radial component of Equation (6) and evaluate it at the surface. The radial component of ${{\boldsymbol{V}}}_{n^{\prime} n}^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ evaluated at the surface is denoted as ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,{\omega }_{\alpha })$, which is the temporal Fourier component of the measured scattered wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$ from n to n' in Zhao & Chou (2013). The index α in ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,{\omega }_{\alpha })$ represents the mode parameters of the incident mode. The radial component of ${{\boldsymbol{V}}}_{\beta }^{(0)n^{\prime} }({\boldsymbol{r}},{\omega }_{\beta })$ evaluated at the surface is denoted as ${{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0}){\phi }_{\beta }(x,y)$, where ${{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0})$ is the radial eigenfunction of n' evaluated at the surface, and ${\phi }_{\beta }(x,y)$ is the orthonormal horizontal eigenfunction of the quiet Sun. With the notation defined previously, the radial component of Equation (6) evaluated at the surface becomes

$$\begin{eqnarray}{{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,{\omega }_{\alpha }) & = & \displaystyle \sum _{\beta }{A}_{\alpha }^{n}{T}_{\beta \alpha }^{n^{\prime} n}\,{{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0})\,{\phi }_{\beta }(x,y)\\ & \equiv & \displaystyle \sum _{\beta }{T}_{\beta \alpha }^{\prime n^{\prime} n}\,{\phi }_{\beta }(x,y).\end{eqnarray} \tag{ 7 }$$

The quantity ${T}_{\beta \alpha }^{\prime n^{\prime} n}\equiv {A}_{\alpha }^{n}{{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0}){T}_{\beta \alpha }^{n^{\prime} n}$ is the expansion coefficient as ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,{\omega }_{\alpha })$ being expanded in terms of the orthonormal horizontal eigenfunctions of radial order n', $\{{\phi }_{\beta }(x,y)\}$.

Similarly, for the incident wave, the radial component of Equation (4) evaluated at the surface yields

$$\begin{eqnarray}{{\rm{\Psi }}}_{n}^{(0)}(x,y,t) & = & \displaystyle \sum _{\alpha }{A}_{\alpha }^{n}{{\rm{\Phi }}}_{\alpha }^{n}({z}_{0}){\phi }_{\alpha }(x,y){e}^{-i{\omega }_{\alpha }t}\\ & \equiv & \displaystyle \sum _{\alpha }{A}_{\alpha }^{\prime n}{\phi }_{\alpha }(x,y){e}^{-i{\omega }_{\alpha }t},\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ is the radial component of the incident wavefunction ${{\boldsymbol{V}}}_{n}^{({\rm{i}})}({\boldsymbol{r}},t)$ evaluated at the surface and ${{\rm{\Phi }}}_{\alpha }^{n}({z}_{0})$ is the radial eigenfunction of n evaluated at the surface. The quantity ${A}_{\alpha }^{\prime n}\equiv {A}_{\alpha }^{n}{{\rm{\Phi }}}_{\alpha }^{n}({z}_{0})$ is the expansion coefficient as ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ being expanded in terms of the horizontal eigenfunctions of radial order n, $\{{\phi }_{\alpha }(x,y){e}^{-i{\omega }_{\alpha }t}\}$. The function ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ is the incident wavefunction measured in Zhao & Chou (2013).

From Equations (7) and (8), the T-matrix element ${T}_{\beta \alpha }^{n^{\prime} n}$ can be expressed as

$$\begin{eqnarray}&&{T}_{\beta \alpha }^{n^{\prime} n}=\left(\displaystyle \frac{{T}_{\beta \alpha }^{\prime n^{\prime} n}}{{A}_{\alpha }^{\prime n}}\right)\left(\displaystyle \frac{{{\rm{\Phi }}}_{\alpha }^{n}({z}_{0})}{{{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0})}\right)\equiv {\bar{T}}_{\beta \alpha }^{n^{\prime} n}\left(\displaystyle \frac{{{\rm{\Phi }}}_{\alpha }^{n}({z}_{0})}{{{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0})}\right),\end{eqnarray} \tag{ 9 }$$

where ${A}_{\alpha }^{\prime n}$ and ${T}_{\beta \alpha }^{\prime n^{\prime} n}$ can be obtained by expanding the measured wavefunctions ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ and ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$ in terms of the temporal Fourier series and orthonormal horizontal eigenfunctions, respectively. The quantity ${\bar{T}}_{\beta \alpha }^{n^{\prime} n}\equiv {T}_{\beta \alpha }^{\prime n^{\prime} n}/{A}_{\alpha }^{\prime n}$ is called the "two-dimensional T-matrix element." The goal of this study is to determine the complex ${\bar{T}}_{\beta \alpha }^{n^{\prime} n}$ from the measured incident and scattered wavefunctions. From ${\bar{T}}_{\beta \alpha }^{n^{\prime} n}$, one can obtain the T-matrix matrix element ${T}_{\beta \alpha }^{n^{\prime} n}$ with Equation (9). The conversion factor ${{\rm{\Phi }}}_{\alpha }^{n}({z}_{0})/{{\rm{\Phi }}}_{\beta }^{n^{\prime} }({z}_{0})$ can be obtained from the eigenfunctions calculated with a standard solar model.

In this study, the horizontal eigenfunction ${\phi }_{\beta }(x,y)=(1/2\pi )\exp ({{ik}}_{x}x+{{ik}}_{y}y)$ is the plane wave in the Cartesian coordinates. Following the setup in Zhao & Chou (2013), the incident wave ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ propagates in the y direction. Thus only the modes with the wavenumber $(n,{k}_{x}=0,{k}_{y})$ exist in the expansion of Equation (8), where k_y relates to ($n,\omega$) through the dispersion relation. Since the frequencies of incident mode α and scattered mode β are the same, the subscripts α and β in frequency ω are dropped hereinafter. If the mode parameters of the incident mode α are $(\omega ,n,0,{k}_{y})$ and the parameters of the scattered mode β are $(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime })$, the two-dimensional T-matrix element ${\bar{T}}_{\beta \alpha }^{n^{\prime} n}$ is explicitly expressed as $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ for clarity in the following discussion. The condition of equal frequency requires that the incident mode (n, k) and the scattered modes $(n^{\prime} ,k^{\prime} )$ correspond to the same frequency ω, where $k=| {k}_{y}|$ and $k^{\prime} ={({k}_{x}^{\prime 2}+{k}_{y}^{\prime 2})}^{1/2}$. The two-dimensional T-matrix element is calculated as

$$\begin{eqnarray}&&\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})\\ &&\quad =\displaystyle \frac{T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})}{A^{\prime} (\omega ,n,0,{k}_{y})},\end{eqnarray} \tag{ 10 }$$

where the coefficient of the incident mode $A^{\prime} (\omega ,n,{k}_{x},{k}_{y})$ is calculated from the measured incident wavefunction of n by the Fourier transform,

$$\begin{eqnarray}&&A^{\prime} (\omega ,n,{k}_{x},{k}_{y})=\int {{\rm{\Psi }}}_{n}^{(0)}(x,y,t){e}^{i(\omega t-{k}_{x}x-{k}_{y}y)}\,{dt}\,{dx}\,{dy},\end{eqnarray} \tag{ 11 }$$

and the coefficient $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ is calculated from the measured scattered wavefunction of n' by the Fourier transform,

$$\begin{eqnarray}&&T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})\\ &&\quad =\displaystyle \int {{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t){e}^{i(\omega t-{k}_{x}^{\prime }x-{k}_{y}^{\prime }y)}\,{dt}\,{dx}\,{dy}.\end{eqnarray} \tag{ 12 }$$

Both $A^{\prime} (\omega ,n,{k}_{x},{k}_{y})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ are complex. The determination of their phase and amplitude will be discussed in the next section.

4.1. Method of Determining the Amplitude and Phase of $A^{\prime} (\omega ,n,{k}_{x},{k}_{y})$

4.2. Method of Determining the Amplitude and Phase of $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$

5.1. Incident and Scattered Wavefunctions

The incident wavefunction and scattered wavefunction used here are from Zhao & Chou (2013). These wavefunctions are computed from SDO/HMI data, which are 4096 × 4096 full-disk Dopplergrams taken at a rate of one image every 45 s. Two active regions, NOAAs 11084 and 11092, are studied. A time series of 12,402 frames (June 29–July 5, 2010) is used for NOAA 11084 and 12,118 frames (July 31–August 6, 2010) for NOAA 11092. Each of these two active regions has only one sunspot, which is approximately round and changes very little during the period of study (see Figure 1 in Zhao & Chou 2013). Thus the scattering problem here is considered time independent.

The total wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}(x,y,t)$ on the surface around the sunspot is measured using a deconvolution method, and the cross-correlation function computed between the signal averaged over a reference line outside the sunspot and the signal at the points in and around the sunspot on the surface. The line-averaged signal and the point signal are formed by the modes of radial orders n and n', respectively (n' can be equal to or different from n). The signal of radial order n is obtained by applying a ridge filter in the $({k}_{x},{k}_{y},\omega )$ domain (Cameron et al. 2008; Zhao et al. 2011a). The same method is applied to the signal of n in the quiet Sun to obtain the incident wavefunction ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$. For the case of n' = n, the scattered wavefunction from n to n is ${{\rm{\Psi }}}_{{nn}}^{({\rm{s}})}(x,y,t)\equiv {{\rm{\Psi }}}_{{nn}}(x,y,t)\,-{{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$. For the case of $n^{\prime} \ne n$, the scattered wavefunction from n to n' is ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)={{\rm{\Psi }}}_{n^{\prime} n}(x,y,t)$. A detailed description of the measurements of scattered wavefunctions is discussed in Zhao et al. (2011a) and Zhao & Chou (2013). As an example, the scattered wavefunctions from n = 4 to n' = 3, 4, 5 are shown in Figure 1. The dimension of each map is $23\buildrel{\circ}\over{.} 4\times 23\buildrel{\circ}\over{.} 4$ (400 × 400 pixels). The size of each pixel is $0\buildrel{\circ}\over{.} 0586$, or 0.712 Mm. The amplitude of the incident wavefunction on the reference line at t = 0 is set to unity. The animation of measured scattered wavefunctions for different n and n' ($n,n^{\prime} =0\mbox{--}5$) can be found in Zhao & Chou (2013).

Zoom In Zoom Out Reset image size

5.2. Dissipation Correction

Both the incident wave and the scattered wave suffer from dissipation. their wave amplitude decreases as they propagate through the Sun. In determining the T-matrix elements, the dissipation effect needs to be corrected because the T-matrix should not contain information of dissipation. In this study, we use the measured dissipation rate of the incident wave in the quiet Sun to correct the dissipation of the scattered wave.

Since dissipation is frequency dependent, the dissipation correction has to be carried out for each frequency. The incident wave propagates only in the y direction. To reduce noise, the incident wavefunction ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ is averaged over x to obtain ${{\rm{\Psi }}}_{n}^{(0)}(y,t)$, which is shown in Figure 2. The wave packet, which starts from the reference line, exists in only a certain region in the (y, t) domain. The upper-left region in Figure 2 has no signal because it corresponds to the spatial region behind the wave packet. There is also no signal in the lower-right corner in Figure 2 because it corresponds to the spatial region ahead of the wave packet; that is, the wave packet has not yet reached this spatial region. The weak residual value in these two regions is noise. Thus we keep only the signals in the region between two red lines in Figure 2 when computing the temporal Fourier transform. Hereinafter, the region with signals is called the signal domain; the region outside the signal domain is called the noise domain. The range of the signal domain is determined with a criterion that the signal is greater than twice the noise level, whose value is computed with an area in the quiet Sun. The value of ${{\rm{\Psi }}}_{n}^{(0)}(y,t)$ is set to zero in the noise domain in computing the temporal Fourier transform. The boundary between the signal and noise domains is smoothed with a Hann window.

Zoom In Zoom Out Reset image size

The temporal Fourier transform transforms ${{\rm{\Psi }}}_{n}^{(0)}(y,t)$ into ${{\rm{\Psi }}}_{n}^{(0)}(y,\omega )$. The amplitude $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ as a function of $(y,\omega )$ is shown in Figure 3(a). It should be mentioned that the mask of the noise domain of ${{\rm{\Psi }}}_{n}^{(0)}(y,t)$ indeed reduces the noise in $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$. The distribution of $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ versus $(y,\omega )$ should be smooth, and its fluctuation is believed to be noise. Thus a two-dimensional smoothing in $(y,\omega )$ is carried out for $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$, with a box running mean of 7 × 7 points. The smoothed $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ is shown in Figure 3(b). For each ω, the smoothed $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ has an approximate exponential decrease in y. This smoothed $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ is used to correct the dissipation in the incident wavefunction and the scattered wavefunction by dividing them with the smoothed $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ if n' = n. For the case of $n^{\prime} \ne n$, the incident wavefunction ${{\rm{\Psi }}}_{n}^{(0)}(x,y,\omega )$ propagates from the reference line to the sunspot, and then part of the incident wave is converted into the scattered wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$. The dissipation correction curve for the scattered wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$ includes two different sections: the smoothed $| {{\rm{\Psi }}}_{n}^{(0)}(y,\omega )|$ is used for y between the reference line and the sunspot center, and the smoothed $| {{\rm{\Psi }}}_{n^{\prime} }^{(0)}(y,\omega )|$ is used for y beyond the sunspot center.

Zoom In Zoom Out Reset image size

5.3. Determination of the Dispersion Relation $k(n,\omega )$ from Incident Waves

For each $(n,\omega )$, there is a corresponding wavenumber $k\equiv {({k}_{x}^{2}+{k}_{y}^{2})}^{1/2}$ through the dispersion relation. As discussed in Section 4, we need a precise value of k in determining the phase of T-matrix elements. Usually the dispersion relation gives ω for given n and k; that is, the dispersion is expressed as a function $\omega (n,k)$. In this study, we need the dispersion relation expressed by a function $k(n,\omega )$ because n and ω are given in our case, as discussed in Section 4. In this subsection, we determine the function $k(n,\omega )$ from the incident wave in the quiet Sun. For the incident wave, $k={k}_{y}$ because k_x = 0. The dissipation-corrected x-averaged incident wavefunction discussed in Section 5.2, denoted as ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(y,\omega )$, is used to determine $k(n,\omega )$. Hereafter, the tilde sign denotes dissipation-corrected wavefunctions. As an example, the real part of ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(y,\omega )$ for n = 3 and $\nu =3.31$ mHz versus y is shown by the black dots in Figure 4, which is very close to a sinusoidal function with unity amplitude. For each ω, ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(y,\omega )$ is fitted to a single-mode wavefunction ${{ae}}^{{iky}+i\phi }$ to determine the wavenumber k. The curve of the fit is the red line in Figure 4. The degree of closeness between the data and the fit indicates that the wavefunction ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(y,\omega )$ of each frequency is very close to a single-mode wavefunction (plane wave) with a wavenumber k.

Zoom In Zoom Out Reset image size

Repeating this analysis for different n and ω yields the dispersion relation $k(n,\omega )$, shown as the black dots in Figure 5. The red line represents the smoothed result with a seven-point box running mean. This smoothed $k(n,\omega )$ is used in the following analysis to determine the y component of the wavenumber for the given $(\omega ,n,{k}_{x})$. For comparison, we include the dispersion relation determined by the conventional way, shown as the blue crosses, in Figure 5(b) (Bogart et al. 2011). The small discrepancy is caused by the conversion between the mode degree ℓ and our wavenumber, which is determined in a planar surface.

Zoom In Zoom Out Reset image size

5.4. Computation of T-matrix Elements

The T-matrix element $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ is computed with Equation (10), where $A^{\prime} (\omega ,n,0,{k}_{y})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ are computed with the Fourier transforms in Equations (11) and (12), respectively. However, the incident wavefunction ${{\rm{\Psi }}}_{n}^{(0)}(x,y,t)$ in Equation (11) and the scattered wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$ in Equation (12) need to be replaced by the dissipation-corrected wavefunction ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(x,y,t)$ and ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$, respectively.

Since the phase of Fourier components depends on the range of the data used in the Fourier transform, the spatial and temporal ranges of the data used in Equations (11) and (12) have to be identical. For the temporal Fourier transform, 825 frames (one frame every 45 s) of wavefunction are used. The corresponding frequency resolution is 0.0269 mHz. For the spatial Fourier transform, the data in an area of 250 × 250 pixels (178 × 178 Mm), indicated by a blue square box in Figure 1(a), are used. The lower boundary of this area is 130 pixels from the reference or 60 pixels (43 Mm) from the sunspot center. There are two reasons for selecting this area to compute the T-matrix elements: First, this area has the strongest scattered wave signals. Second, it is outside the source of the scattered waves, the sunspot. Thus no more scattered wave is generated in this area, and there is no magnetic contamination in measured Doppler signals. The corresponding resolution of k_x and k_y is Δk = 0.035 rad Mm⁻¹. In the following discussion, k_x and k_y are in units of ${\rm{\Delta }}k$.

For the computation of $A^{\prime} (\omega ,n,0,{k}_{y})$, the signal domain of the incident wave is the same as that of ${{\rm{\Psi }}}^{(0)}(y,t)$, discussed in Section 5.2. The dissipation-corrected wavefunction ${\tilde{{\rm{\Psi }}}}_{n}^{(0)}(x,y,t)$ is masked before computing $A^{\prime} (\omega ,n,0,{k}_{y})$ in Equation (11). The computation of $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ is more complicated because the signal domain of the scattered wave is different from that of the incident wave.

The first step in computing $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ with Equation (12) is the temporal Fourier transform. The signal domain of ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$ is determined by the x-averaged scattered wavefunction ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(y,t)$, without the dissipation correction. The method is similar to the determination of the signal domain of ${{\rm{\Psi }}}_{n}^{(0)}(y,t)$ in Figure 2. The noise domain of ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,t)$ is masked before computing the temporal Fourier transform; the Fourier component is denoted as ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$. As an example, the scattered wavefunctions from n = 3 to n' = 3 at two frequencies are shown in Figure 6.

Zoom In Zoom Out Reset image size

The second step is computing the spatial Fourier transform of ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$. To reduce noise, only the signal in the spatial signal domain is kept in the spatial Fourier transform. The spatial signal domain is also determined with the criterion that the signal is greater than twice the noise level. The noise level is the fluctuation (standard deviation) of the wavefunction determined from an area of 250 × 250 pixels in the quiet Sun for each $(\omega ,n,n^{\prime} )$. The signal is the amplitude of the wavefunction, and the signal domain is determined as follows. The x and y ranges of the signal domain are determined with the y-average of $| {{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega ){| }^{2}$ and the x-average of $| {{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega ){| }^{2}$, respectively. It is noted that the signal domain is different for different $(\omega ,n,n^{\prime} )$. For the example in Figure 6, the signal domain of ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$ is marked with a red rectangular box. The noise domain is again masked before performing the spatial Fourier transform. The result is $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$. It is noted that the amplitude of $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ is reduced because of masking in the y domain. This is corrected by multiplying a geometric factor to compensate the mask in the y domain. This simple correction is justified because the amplitude of dissipation-corrected ${\tilde{{\rm{\Psi }}}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$ is approximately constant in y.

In the measurement of the scattered wavefunctions, Zhao & Chou (2013) averaged the wavefunctions associated with the incident waves of different propagation directions to enhance the signal-to-noise ratio (S/N). Thus the direction information on the sunspot is lost, and the scattered wavefunction should be symmetric with respect to the incident direction, the y axis; that is, $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ should be symmetric with respect to ${k}_{x}^{\prime }=0$. The asymmetry in the measured $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ should be caused by noise. Therefore $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})$ is symmetried with respect to ${k}_{x}^{\prime }=0$ in the following discussion. For example, Figure 7 shows the amplitudes $| A^{\prime} (\omega ,n,{k}_{x},{k}_{y})|$ and $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})|$ at $\nu =3.07$ mHz for n = 3 and $n^{\prime} =$ 2, 3, 4. The incident wave amplitude $| A^{\prime} (\omega ,n,{k}_{x},{k}_{y})|$ is shown in Figure 7(a). The signal is significant only at k_x = 0. The very small and approximately constant value at ${k}_{x}\ne 0$ on the half ring is noise, since the incident wave propagates only in the y direction. The ring corresponds to $\nu =3.07$ mHz. For the scattered wave, $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})|$ from n = 3 to $n^{\prime} =2,3,4$ are shown in Figures 7(b), (c), (d), respectively. It is noted that the color scales are different in different panels. The amplitude from n = 3 to n' = 3 is much stronger than that from n = 3 to $n^{\prime} =2$ or 4. The amplitude peaks at ${k}_{x}^{\prime }=0$, and decreases rapidly with ${k}_{x}^{\prime }$ because the scattered wave propagates predominantly in the y direction. For most of the modes on the half ring in Figures 7(b), (c), (d), its value is noise. Only a few modes around ${k}_{x}^{\prime }=0$ are above the noise level; they represent the scattered signals. We will determine the amplitude and phase only for these ${k}_{x}^{\prime }$'s.

Zoom In Zoom Out Reset image size

5.5. Determination of the Amplitude of T-matrix Elements

From Equation (10), the amplitude of the T-matrix element $| \bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})|$ equals the ratio of $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{y})|$ to $| A^{\prime} (\omega ,n,0,{k}_{y})|$. As discussed in Section 4.1, the y component of mode wavenumber determined by the dispersion relation for the given $(n,\omega )$, denoted by ${k}_{{y}_{0}}$, may not be at a grid point. For example, the amplitude $| A^{\prime} (\omega ,n,0,{k}_{y})|$ corresponding to Figure 7(a) versus k_y is shown in Figure 8(a). From the dispersion relation shown in Figure 5, the y component of the mode wavenumber is ${k}_{{y}_{0}}=10.71{\rm{\Delta }}k$, a non-integer in units of the wavenumber resolution ${\rm{\Delta }}k$, for n = 3 and $\nu =3.07$ mHz. The value of ${k}_{{y}_{0}}$ is indicated by a vertical red line in Figure 8(a). In the Appendix, we discuss how to determine the amplitude and phase of a signal-mode signal if its wavenumber is not exactly equal to one of the Fourier wavenumbers. The method is demonstrated by a simulation using a one-dimensional signal of 250 points, the same number of data points in our spatial Fourier transform. The distribution in Figure 8(a) is similar to the simulation in Figure 16(a). The value of $| A^{\prime} (\omega ,n,0,{k}_{y})|$ peaks at the point nearest ${k}_{{y}_{0}}$, and decreases rapidly away from ${k}_{{y}_{0}}$.

Zoom In Zoom Out Reset image size

From the discussion in Section 4.1 and the Appendix, the power of the incident mode $(\omega ,n,0,{k}_{{y}_{0}})$ is $| A^{\prime} (\omega ,n,0,{k}_{{y}_{0}}){| }^{2}={\sum }_{{k}_{y}}| A^{\prime} (\omega ,n,0,{k}_{y}){| }^{2}$. Since $| A^{\prime} (\omega ,n,0,{k}_{y})|$ is significant only in the neighboring points of ${k}_{{y}_{0}}$, in practice, only those amplitudes above one-tenth of the peak amplitude are included in the summation, to avoid including noise in the summation. The result of $| A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})|$ for this case is 0.973, very close to unity. This is consistent with the fact that the amplitude of a dissipation-corrected incident wave is unity, as discussed in Section 5.2. In the following analysis, $| A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})|$ is set to unity.

For the scattered mode, the value of $k^{\prime}$ is determined from the dispersion relation $k^{\prime} (n^{\prime} ,\omega )$. For a given ${k}_{x}^{\prime }$ at the grid point, the y component of the mode wavenumber ${k}_{{y}_{0}}^{\prime }$ is computed as ${k}_{{y}_{0}}^{\prime }={(k{^{\prime} }^{2}-{k}_{x}^{\prime 2})}^{1/2}$, which may not be at a grid point. As an example, for the scattering from n = 3 to n' = 3, $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{{y}_{0}})|$ corresponding to Figure 7(c) vs. ${k}_{y}^{\prime }$ for ${k}_{x}^{\prime }=0$, ${\rm{\Delta }}k$, $2{\rm{\Delta }}k$ are shown in Figures 9(a), (c), (e), respectively. For ${k}_{x}^{\prime }=0$, ${\rm{\Delta }}k$, $2{\rm{\Delta }}k$, ${k}_{{y}_{0}}^{\prime }=10.71{\rm{\Delta }}k$, $10.67{\rm{\Delta }}k$, $10.52{\rm{\Delta }}k$, respectively. The value of ${k}_{{y}_{0}}^{\prime }$ is marked by the vertical red line in each panel. The value of $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{{y}_{0}})|$ peaks at the point nearest ${k}_{{y}_{0}}^{\prime }$, and decreases rapidly away from ${k}_{{y}_{0}}^{\prime }$. The peak value for each ${k}_{x}^{\prime }$ decreases with ${k}_{x}^{\prime }$ and is rapidly lost in noise. Thus we can determine the T-matrix elements only for small ${k}_{x}^{\prime }$ values, such as ${k}_{x}^{\prime }=0-3{\rm{\Delta }}k$, for most cases.

Zoom In Zoom Out Reset image size

Like the incident mode, the amplitude $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})|$ equals the the square root of the summation of $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{y}^{\prime };\omega ,n,0,{k}_{{y}_{0}}){| }^{2}$ over ${k}_{y}^{\prime }$ in the neighborhood of ${k}_{{y}_{0}}^{\prime }$. Since $| A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})| =1$, the amplitude of the T-matrix element $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ equals $| T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})|$ from Equation (10).

For case of $n\ne n^{\prime}$, the signals are much weaker and the results are nosier. As an example, Figure 10 shows the plots similar to Figure 9, but from n = 3 to $n^{\prime} =4$ at $\nu =3.29$ mHz. The amplitudes are much smaller than those in Figure 9, and the phases are also nosier.

Zoom In Zoom Out Reset image size

5.6. Determination of the Phase of T-matrix Elements

The phase of a spatial Fourier component depends on the spatial range of the data used in computing the Fourier transform. The phase of the Fourier component is the phase of the corresponding mode signal at the first spatial point of the data. If a different range of data is used in computing the Fourier component, its phase will be shifted by $k{\rm{\Delta }}x$, where k is the wavenumber of the Fourier component and ${\rm{\Delta }}x$ is the shift of the first point of the data used in the computation. Thus, for different wavenumbers, the phase shifts are different. In our analysis, an area indicated by the blue square box in Figure 1(a) is used to compute the Fourier transform. The phase of $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ or $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ determined in Section 5 is the phase of the mode at $(x=0,y=0$), which is located at the lower-left corner of the blue square box in Figure 1(a). A question arises: Since the phase determined previously depends on the area used in the computation, how do we interpret the phase measured in Section 5?

For the diagonal elements of the T-matrix, the scattered mode and the incident mode are the same; that is, n' = n, ${k}_{x}^{\prime }={k}_{x}$, and ${k}_{{y}_{0}}^{\prime }={k}_{{y}_{0}}$. The phase shifts for $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ are equal, even if a different area is used. Thus the phase of the T-matrix element $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ would be independent of the area used. For the off-diagonal elements, the scattered mode and the incident mode are different; the phase shifts for $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ are different because their wavenumbers are different. Thus the phase of $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ is different for a different area. How can we resolve this issue? A reasonable approach is to shift the phases of $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ to a reasonable common reference point before computing the phase difference.

The scattered waves are generated by the sunspot. Although every part of the sunspot contributes to the generation of the scattered waves, we can considered the sunspot center a reasonable reference point for all scattered modes for the purpose of comparing the phases of different modes. Thus the phases of $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ and $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ are shifted to the sunspot center; that is, the phase measured in Section 5 is extrapolated to the sunspot center. After extrapolation, the phase difference between $T^{\prime} (\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ and $A^{\prime} (\omega ,n,0,{k}_{{y}_{0}})$ is computed to obtain the phase of $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$.

After being generated, the scattered wave propagates inside the sunspot with a speed different from that in the quiet Sun before it leaves the sunspot. Similarly, the incident mode also propagates with a different speed inside the sunspot before it generates the scattered mode. Therefore the phase of $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ obtained here indicates the following: as the incident mode propagates through the sunspot, various scattered modes are generated inside the sunspot. The phase of $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ measured in this study does not solely reflect the phase difference between the scattered and incident modes as the scattered mode is generated. It also includes the phase change due to the change in wave speed as the incident and scattered modes propagate inside the sunspot. Note that the phase can always be added or subtracted by 360°.

The amplitude and phase of various elements $\bar{T}(\omega ,n^{\prime} ,{k}_{x}^{\prime },{k}_{{y}_{0}}^{\prime };\omega ,n,0,{k}_{{y}_{0}})$ for NOAA 11092 are shown in Figures 11–14. For the case of $n=0\mbox{--}2$, the amplitude of $n^{\prime} \ne n$ is too small to be detected. For $n=3\mbox{--}5$, only $n^{\prime} =n\pm 1$ is detectable. For n' = n, ${k}_{x}^{\prime }=0-3{\rm{\Delta }}k$ are plotted, while for $n^{\prime} \ne n$, only ${k}_{x}^{\prime }=0-2{\rm{\Delta }}k$ are plotted, because the results of other greater ${k}_{x}^{\prime }$'s are noisy. For NOAA 11084, we show only one example of n = 4 in Figure 15 for comparison. The error bar of amplitude is estimated with the amplitude in the noise domain. The error bar of phase is estimated as the arcsine of the ratio of amplitude error to measured amplitude. Overall, the errors of amplitude and phase are consistent with the fluctuations of amplitude and phase. In general, the amplitude error increases with frequency, while the phase error is larger at the lower frequency. For the case of n' = n, both amplitude and phase decrease with ${k}_{x}^{\prime }$. However, for the case of $n^{\prime} \ne n$, this trend is less apparent because the differences among different ${k}_{x}^{\prime }$ are less than the error bars. Figures 12–15 show that although the amplitude of $n^{\prime} \ne n$ is smaller than that of n' = n, they are of the same order of magnitude. This indicates that the mode mixing between adjacent radial orders is significant for some radial orders.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The penumbra radii of NOAAs 11084 and 11092 are about 15 and 20 pixels, respectively. The comparison of Figures 13 and 15 shows that the amplitude of NOAA 11092 is larger than that of NOAA 11084, and the amplitude ratio is approximately equal to the square of sunspot radius ratio. For the phase, the values of two sunspots are rather close. This indicates that the amplitude correlates with the sunspot size, while the phase is insensitive to the sunspot size.

The T-matrix element can be obtained from these measured two-dimensional T-matrix element with Equation (9). Since the radial wavefunction ${{\rm{\Phi }}}_{\alpha }^{n}(z)$ depends only on n and k, the conversion factor is ${{\rm{\Phi }}}_{{nk}}({z}_{0})/{{\rm{\Phi }}}_{n^{\prime} k^{\prime} }({z}_{0})={{\rm{\Phi }}}_{n\omega }({z}_{0})/{{\rm{\Phi }}}_{n^{\prime} \omega }({z}_{0})$, where the equality is based on the dispersion relation and $\omega ^{\prime} =\omega$. The conversion factor is a real number. Thus the phase of the T-matrix element is equal to the phase of the two-dimensional T-matrix element. For the case of n' = n, the conversion factor is unity for any ${k}_{x}^{\prime }$. For the case of $n^{\prime} \ne n$, one needs to compute the radial eigenfunctions of incident mode and scattered mode at the surface from a standard solar model to obtain the conversion factor. The S-matrix is the sum of the unity matrix and the T-matrix.

One cannot make a direct comparison of the amplitude and phase of T-matrix elements measured here with the results of other methods, such as Hankel analysis (Braun 1995), time–distance analysis (Chou et al. 2009a, 2009b, 2009c), and phenomenological-model analysis (Yang et al. 2012), because the definitions of measured quantities in different methods are different. It is also difficult to find a clear and useful relation between the measured quantities of different methods for comparison.

In our analysis, the sunspot is assumed to be time independent, so there is no mode mixing between different frequencies. If the sunspot is time-dependent with a frequency Ω, a mode would mix with other modes with a frequency shift of Ω after interacting with the sunspot (Schiff 1971). Since our two sunspots change little in the observing period (about 7 days), the characteristic timescale for the sunspots is longer than 10 days; the corresponding frequency ${\rm{\Omega }}\lt {10}^{-3}$ mHz. The frequency bin in our analysis is 0.0269 mHz, much greater than Ω. Therefore the assumption of time independence for the sunspots is justified.

The S-matrix (or T-matrix) elements provide the most complete observational or experimental information on the interaction responsible for the scattering and on the properties of the scatterer. In this study, we use the scattered wavefunctions measured with cross correlation and deconvolution to compute the various T-matrix elements. These T-matrix elements provide rich information on the interaction of the acoustic waves with the sunspots and the properties of the sunspots. The properties of these measured T-matrix elements and their implication for the sunspots and the interaction will be the subject of a follow-up paper to avoid making this paper too long. Here we discuss only one implication from Figures 11–15.

All amplitudes of two-dimensional T-matrix elements in Figures 11–15 have a general trend: They increase with frequency for each n and ${k}_{x}^{\prime }$. Since the horizontal phase speed anti-correlates with frequency for each n, the measured amplitude decreases with the horizontal phase speed for each n. This property provides information on the interaction region, which is responsible for generating the amplitude of the T-matrix element. It suggests that the interaction strength decreases with depth from a depth of about 1 Mm, which corresponds to the the lower turning point of the shallowest mode cavity in our measurements, or a horizontal phase speed of $2\pi \nu /{\ell }=0.02$ mHz. This is explained as follows. The upper turning point of mode cavity changes little with horizontal phase speed, while the depth of its lower turning point increases significantly with horizontal phase speed; that is, a mode samples more and more depth range as its horizontal phase speed increases. If the interaction strength decreases with depth, the influence of the interaction on the mode would decrease as the mode cavity size increases. Therefore, the amplitude of the T-matrix element decreases with the horizontal phase speed. The detailed depth distribution of the interaction strength requires the inversion of measured amplitude. In Figures 11–15, the dependence of the phase on frequency or horizontal phase speed is very different from that of the amplitude and has more complicated features. It suggests that the depth distribution of the interaction region responsible for the phase is different from that responsible for the amplitude. The discussion on this issue is left for a follow-up paper.

In this study, we compute the T-matrix elements with the incident modes of k_x = 0 only. Since our measured scattered wavefunctions are for an axisymmetric sunspot, using the incident waves of ${k}_{x}\ne 0$ would not yield additional information. In fact, for the given $(\omega ,n,n^{\prime} )$, the T-matrix elements with an incident mode of ${k}_{x}\ne 0$ can be obtained from those with the incident mode of k_x = 0 by a rotation transformation, for the case of an axisymmetric scatterer.

Since the data used in our analysis are a time series of 825 frames (about 10.5 hr), and each frame is 250 × 250 pixels (178 × 178 Mm), the resolutions in the frequency and wavenumber domains are 0.0269 mHz and 0.035 rad Mm⁻¹, respectively. Each incident "mode" or scattered "mode" discussed in this study is not really a single mode; instead, it is some kind of average over several adjacent modes within a bin. Each T-matrix element represents the transition between two such averaged modes. Since the wavefunctions used in this study are measured with cross correlation, their phase is relative to the phase of the reference line (initial location of the incident wave), and no random phase is involved as the observed Doppler data (Zhao et al. 2011a). The phase and amplitude of adjacent modes are expected to vary smoothly. Thus the averaged mode should represent the properties of the modes within the bin.

In this study, the T-matrix elements are calculated using the measured scattered wavefunctions through the Fourier transform, as in Equations (11) and (12). It is worth mentioning that the Fourier components in Equations (11) and (12) can also be directly computed from the observed Doppler signals, instead of using the measured scattered wavefunctions. The scattered wavefunction is computed with cross correlation and deconvolution, as in Equation (2) in Zhao & Chou (2013). Thus the temporal Fourier component of the scattered wavefunction from n to n' can be computed with

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )={A}^{-1}\displaystyle \frac{{\bar{{\rm{\Psi }}}}_{n}^{* }({y}_{0},\omega )}{| {\bar{{\rm{\Psi }}}}_{n}({y}_{0},\omega )| }{{\rm{\Psi }}}_{n^{\prime} }(x,y,\omega ),\end{eqnarray} \tag{ 13 }$$

where ${\bar{{\rm{\Psi }}}}_{n}({y}_{0},\omega )$ is the temporal Fourier component of the line-averaged observed Doppler signal of n at the reference line $y={y}_{0}$, and ${{\rm{\Psi }}}_{n^{\prime} }(x,y,\omega )$ is the temporal Fourier component of the observed Doppler signal of n' at the point (x, y). Since the spatial dependence of ${{\rm{\Psi }}}_{n^{\prime} n}^{({\rm{s}})}(x,y,\omega )$ is from ${{\rm{\Psi }}}_{n^{\prime} }(x,y,\omega )$, the Fourier component in Equation (12) can be directly computed from the spatial Fourier transform of ${{\rm{\Psi }}}_{n^{\prime} }(x,y,\omega )$, which is the temporal Fourier component of the observed Doppler signal. A question arises: Why do we go through the extra step in this study to calculate the T-matrix elements using the measured scattered wavefunctions instead of directly using the observed Doppler signals? The reason is that the measured scattered wavefunction has information on the propagation of the scattered wave packet in spacetime. The scattered wave packet occupies only a limited domain in spacetime. It is easy to identify the signal domain and mask the noise domain before the Fourier transform, as discussed in Sections 5.2 and 5.4. The masking of the noise domain significantly reduces the noise in the computed Fourier components to allow for the determination of the amplitude and phase of the T-matrix elements.

The interpretation of measured T-matrix elements will be one of the main topics in a follow-up paper. Here we give only a brief discussion on the connection between the measured T-matrix elements and the interaction of the wave with a sunspot. First we formulate equations for the interaction between the acoustic waves and a sunspot, which is described by a source function. Then the T-matrix element is expressed in terms of an integral of the source function and the quiet-Sun Green's function. The source function, which relates to the subsurface properties of the sunspot, can be inferred from the measured T-matrix elements.

For the quiet Sun, the spatial part of the wavefunction for mode α, ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$, is the solution of the time-independent wave equation

$$\begin{eqnarray}&&[{{ \mathcal L }}^{(0)}({\boldsymbol{r}})-{\omega }_{\alpha }^{2}]{{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })=0,\end{eqnarray} \tag{ 14 }$$

where ${{ \mathcal L }}^{(0)}({\boldsymbol{r}})$ is the time-independent wave operator of the quiet Sun, which involves vectors (Unno et al. 1989). The incident wave is one of the solutions of Equation (14)—for example, ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$. In the presence of the sunspot, the operator ${{ \mathcal L }}^{(0)}({\boldsymbol{r}})$ is modified to ${{ \mathcal L }}^{(0)}({\boldsymbol{r}})+{{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$, where the perturbation of wave operator ${{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$ describes the interaction between the wave and the sunspot. After interacting with a time-independent sunspot, the frequency remains unchanged and the incident wave ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$ is converted into ${{\boldsymbol{V}}}_{\alpha }({\boldsymbol{r}},{\omega }_{\alpha })={{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })+{{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$, which satisfies

$$\begin{eqnarray}&&[{{ \mathcal L }}^{(0)}({\boldsymbol{r}})+{{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})-{\omega }_{\alpha }^{2}][{{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })+{{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })]=0.\end{eqnarray} \tag{ 15 }$$

Equation (15) can also be expressed as an inhomogeneous equation for the scattered wavefunction ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$

$$\begin{eqnarray}&&[{{ \mathcal L }}^{(0)}({\boldsymbol{r}})-{\omega }_{\alpha }^{2}]{{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })={{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}}),\end{eqnarray} \tag{ 16 }$$

where ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}})\equiv -{{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}}){{\boldsymbol{V}}}_{\alpha }({\boldsymbol{r}},{\omega }_{\alpha })$ is called the source function which generates the scattered wave ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$. For the first-order Born approximation, ${{\boldsymbol{V}}}_{\alpha }({\boldsymbol{r}},{\omega }_{\alpha })$ in the expression for the source function is replaced by the incident wave ${{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$.

Equation (16) can be solved with the Green's function of the quiet Sun. Since the operator ${{ \mathcal L }}^{(0)}({\boldsymbol{r}})$ involves vectors, the solution ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ and the source ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}})$ have, in general, different directions. Thus the Green's function is a second-rank tensor in three-dimensional space (Jackson 1999). In this study, we use the Greek alphabet to label the modes and the Latin alphabet to denote the spatial directions. The solution ${{\boldsymbol{V}}}_{\alpha }^{({\rm{s}})}({\boldsymbol{r}},{\omega }_{\alpha })$ in direction i can be expressed as

$$\begin{eqnarray}&&{V}_{\alpha }^{({\rm{s}})i}({\boldsymbol{r}})=\int \displaystyle \sum _{j}\,{G}^{{ij}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){J}_{\alpha }^{j}({\boldsymbol{r}}^{\prime} ){d}^{3}r^{\prime} ,\end{eqnarray} \tag{ 17 }$$

where ${J}_{\alpha }^{j}({\boldsymbol{r}}^{\prime} )$ is the source function in direction j, and the second-rank tensor ${G}^{{ij}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )$ is the Green's function satisfying

$$\begin{eqnarray}&&[{{ \mathcal L }}^{(0)}({\boldsymbol{r}})-{\omega }_{\alpha }^{2}]({\boldsymbol{r}}){G}^{{ij}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )={\delta }_{{ij}}\,\delta ({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} ),\end{eqnarray} \tag{ 18 }$$

where ${\delta }_{{ij}}$ is the Kronecker delta, and $\delta ({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} )$ the Dirac δ-function.

With Equation (17), the T-matrix element ${T}_{\beta \alpha }$ can be expressed in terms of the Green's function as

$$\begin{eqnarray}{T}_{\beta \alpha } & = & \displaystyle \int \rho ({\boldsymbol{r}})\displaystyle \sum _{i}\,{V}_{\beta }^{(0)i* }({\boldsymbol{r}}){V}_{\alpha }^{({\rm{s}})i}({\boldsymbol{r}})\,{d}^{3}r\\ & = & \displaystyle \int \rho ({\boldsymbol{r}})\displaystyle \sum _{{ij}}\,{V}_{\beta }^{(0)i* }({\boldsymbol{r}}){G}^{{ij}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){J}_{\alpha }^{j}({\boldsymbol{r}}^{\prime} )\,{d}^{3}{{rd}}^{3}r^{\prime} .\end{eqnarray} \tag{ 19 }$$

The eigenfunction ${V}_{\beta }^{(0)i}({\boldsymbol{r}})$ and the Green's function ${G}^{{ij}}({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )$ of the quiet Sun can be solved numerically with a solar model. With the measured ${T}_{\beta \alpha }$, Equation (19) can be used to infer the source function ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}}^{\prime} )$. For the forward study, given a modeled ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}}^{\prime} )$, the T-matrix elements can be computed to compare with the measured values. One can also apply inversion on Equation (19) to estimate the source function ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}}^{\prime} )$. From ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}})=-{{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}}){{\boldsymbol{V}}}_{\alpha }^{(0)}({\boldsymbol{r}},{\omega }_{\alpha })$, one can infer ${{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$ from the known ${{\boldsymbol{J}}}_{\alpha }({\boldsymbol{r}})$.

The ultimate goal is to determine the subsurface properties of the sunspot, which include the magnetic field, the thermal perturbation, and the flow in the sunspot region. These properties are the perturbations to the quiet-Sun structure. They would appear in the basic equations, including the mass-conservation equation, the momentum equation, and the energy equation, and in turn would appear in the perturbation of wave operator ${{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$. Thus one needs a theory to connect these properties with ${{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$ in order to learn these properties from the inferred ${{ \mathcal L }}^{({\rm{s}})}({\boldsymbol{r}})$.

The SDO/HMI data are provided by the HMI team at Stanford University. M.H.Y. and D.Y.C. are supported by the MOST of ROC under grant NOST-105-21120-M-007-001-MY3.

In this appendix we discuss the determination of the phase and amplitude of a single-mode signal whose wavenumber is not exactly equal to one of the Fourier wavenumbers. The method is demonstrated with an example of a one-dimensional single-mode discrete signal of N points,

$$\begin{eqnarray}&&{f}_{q}=a\,{e}^{i2\pi {p}_{0}q/N}{e}^{i{\phi }_{0}},\end{eqnarray} \tag{ 20 }$$

where a is the amplitude and ${\phi }_{0}$ is the phase; q and p₀ label the position and the wavenumber, respectively; and q is an integer, but p₀ may not be an integer. We discuss how the amplitude a and phase ${\phi }_{0}$ can be determined from the Fourier transform. The discrete Fourier transform of f_q is

$$\begin{eqnarray}&&{F}_{p}=\displaystyle \frac{1}{N}\displaystyle \sum _{q=0}^{N-1}\,{f}_{q}\,{e}^{-i2\pi {pq}/N}=\displaystyle \frac{1}{N}\,a\,{e}^{i{\phi }_{0}}\displaystyle \sum _{q=0}^{N-1}\,{e}^{i2\pi ({p}_{0}-p)q/N}.\end{eqnarray} \tag{ 21 }$$

Here the sign of the Fourier transform is consistent with that in Equations (11) and (12). If p₀ is an integer (at a grid point), ${F}_{p}=a\,{e}^{i{\phi }_{0}}{\delta }_{p,{p}_{0}}$. The amplitude and phase of f_q are simply the amplitude and phase of the Fourier component ${F}_{{p}_{0}}$.

If p₀ is a non-integer (not at a grid point), F_p in Equation (21) becomes

$$\begin{eqnarray}{F}_{p} & = & \displaystyle \frac{a}{N}\,{e}^{i{\phi }_{0}}\displaystyle \frac{1-{e}^{i2\pi ({p}_{0}-p)}}{1-{e}^{i2\pi ({p}_{0}-p)/N}}\\ & = & \displaystyle \frac{a}{N}\,\displaystyle \frac{\sin [\pi ({p}_{0}-p)]}{\sin [\pi ({p}_{0}-p)/N]}{e}^{i{\phi }_{0}}{e}^{i\pi ({p}_{0}-p)(1-1/N)}.\end{eqnarray} \tag{ 22 }$$

The Fourier component F_p is nonzero everywhere in the domain of p, but its amplitude peaks at the integer nearest p₀. For example, Figure 16(a) shows $| {F}_{p}|$ vs. p for the case of a = 1, ${\phi }_{0}=30^\circ$, p₀ = 11.4, and N = 250 based on Equation (22). The vertical red line marks $p={p}_{0}$. Note that the data analysis in Section 5.4 also has N = 250.

Zoom In Zoom Out Reset image size

The summation of power in the Fourier domain is

$$\begin{eqnarray}&&\displaystyle \sum _{p=0}^{N-1}| {F}_{p}{| }^{2}=\displaystyle \sum _{p=0}^{N-1}\,\displaystyle \frac{{a}^{2}}{{N}^{2}}\,\displaystyle \frac{{\sin }^{2}[\pi ({p}_{0}-p)]}{{\sin }^{2}[\pi ({p}_{0}-p)/N]}={a}^{2}.\end{eqnarray} \tag{ 23 }$$

This is consistent with Parseval's theorem. Therefore the amplitude of f_q is equal to the square root of the summation of power over the Fourier domain.

The phase of F_p is the argument of F_p (between $-\pi$ and π), denoted as $\arg ({F}_{p})$. Assume p₀ is between two integers $p^{\prime}$ and $p^{\prime} +1$ (i.e., $p^{\prime} \lt {p}_{0}\lt p^{\prime} +1$). From Equation (22), the phases of F_p at $p^{\prime}$ and $p^{\prime} +1$ are

$$\begin{eqnarray}&&\arg ({F}_{p^{\prime} })={\phi }_{0}+\pi ({p}_{0}-p^{\prime} )(1-1/N),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\arg ({F}_{p^{\prime} +1})={\phi }_{0}+\pi ({p}_{0}-p^{\prime} -1)(1-1/N).\end{eqnarray} \tag{ 25 }$$

The phase difference between $p^{\prime}$ and $p^{\prime} +1$ is $\arg ({F}_{p^{\prime} +1})\,-\arg ({F}_{p^{\prime} })=-\pi (1-1/N)\approx -\pi$, if $N\gg 1$. For our case, N = 250, the approximation creates an error of 4 × 10⁻³. For $p\lt p^{\prime}$, $\arg ({F}_{p})\approx \arg ({F}_{p^{\prime} })$, and $\arg ({F}_{p})\approx \arg ({F}_{p^{\prime} +1})$ for $p\gt p^{\prime} +1$. That is, $\arg ({F}_{p})$ is approximately constant on either side of p₀, and the value drops approximately π across p₀.

From Equations (24) and (25), ${\phi }_{0}$ is the linear interpolation between $\arg ({F}_{p^{\prime} +1})$ and $\arg ({F}_{p^{\prime} })$ at p₀. Figure 16(b) shows the phase of F_p vs. p for the case of a = 1, ${\phi }_{0}=30^\circ$, and p₀ = 11.4, based on Equation (22). The vertical red line marks $p={p}_{0}$. The blue cross sign marks the intersection of the black line and the red line; that is, it is exactly the interpolation of phase at p₀. Its value is equal to the phase of f_q, ${\phi }_{0}=30^\circ$.

In summary, for a non-integer wavenumber p₀, the signal in the Fourier domain spreads around p₀, and there is a 180^◦ phase drop from the grid point smaller than p₀ to the grid point greater than p₀. The phase of f_q (the phase of F_p at p₀) can be computed with the linear interpolation between the phases of F_p at two grid points nearest p₀. The amplitude of f_q is the square root of the summation of $| {F}_{p}{| }^{2}$ over p. In practice, we need to sum $| {F}_{p}{| }^{2}$ over only a small range around p₀, because $| {F}_{p}{| }^{2}$ decreases rapidly away from p₀.