3 Numerical experiment

As a second step in our approach to gain knowledge about the small effect of velocity fields on the solar oscillations we will use perturbation theory to carry out a numerical experiment.

We create artificial oscillation frequencies, that cover the whole $(l,~\nu)$ diagram in a range $0\le l\le 150$ for the harmonic degree and all respective azimuthal orders. As the availability of modes is restricted to the observation we also use only those multiplets which are observed, cf. Fig. 1. Finally these artificial frequencies are inverted for differential rotation.

3.1 Simulated oscillation data

The process of creating the data is the following: the frequencies are determined as eigenfrequencies of a standard solar model. We use the "Model S'' from J. Christensen-Dalsgaard to determine the degenerate oscillations.

Afterwards we perform a perturbation calculation to take the influence of differential rotation into account. As an input model for the differential rotation we use an inversion result, that is an average of six three-month periods of GONG data, beginning in March 1996 and ending in January 1998 (GONG months 10-27). The profile is similar to that shown in Howe et al. (2000). The result is the splitting of the degenerate oscillations into frequency multiplets, i.e. the frequency has now become a function of the azimuthal order m (Fig. 2).

In the next step we add one component $\vec{P}_s^t(\vec{r})$ of the poloidal flow field and use quasi-degenerate perturbation theory. The model for the flow component $\vec{P}_s^t(\vec{r})$ was given by Eq. (9). We assume that ${\rm d}\vec{P}_s^t/{\rm d}\tau=0$ and hence $\nabla\cdot(\rho_0(r)\vec{P}_s^t)=0$ so that v_s^t(r) can be expressed in terms of u_s^t(r)

$$v_s^t(r)=\frac{\partial_r(r^2\rho_0u_s^t(r))}{\rho_0 r s(s+1)}\cdot$$ (18)

The quantity u_s^t(r) is taken to be

$$u_s^t(r)= u_{\rm char}{\frac{4(R-r)(r- r_{\rm conv})}{(R-r_{\rm conv})^2}}\ \ \ {\rm for}\ \ r_{\rm conv}\le r \le R,$$ (19)

where $u_{\rm char}$ is a characteristic parameter that allows us to control the strength of the flow field component, R and $r_{\rm conv}$are the solar radius and the bottom of the convection zone. The resulting depth dependence of u_s^t(r) and v_s^t(r) is shown in Fig. 6. The radial component u_s^t(r) of the flow vanishes at the top and the bottom of the convection zone and has its maximum value $u_{\rm char}$ in the middle of the convection zone. In between we have chosen a simple parabolic dependence of u_s^t on the radius r.

The horizontal expansion coefficient v_s^t(r) changes its sign, so that a laminar flow on closed streamlines is established. The difference in amplitude of v_s^t between top and bottom of the convection zone reflects the density stratification. At the top the amplitude is large because of the low density, at the base of the convection zone the density is higher, therefore only a small amplitude in v_s^t is necessary to conserve the mass flow.

Figure 6: The depth-dependent expansion coefficients ust(r) and vst(r) of the radial and horizontal component of the poloidal velocity fields used (s=8).

In general such poloidal velocity fields have the properties of Bénard convection, i.e. cellular flows with upwards streaming matter in the center and downward streams at the borders of the cells. The geometry and the amplitude of the flow field is fixed by the parameters s, t, and $u_{\rm char}$. We discuss in Sect. 4 the results for the two extreme cases s=t (sectoral cells) and t=0 (zonal cells). A visualization of such field components is given in Fig. 7. The plots on the left and in the middle show streamlines of an equatorial and a meridional section through the poloidal flow field s=8, t=8, and the right plot shows the meridional section through the poloidal flow field s=8, t=0. Since a poloidal field with t=0 does not have components in direction of the rotation, an equatorial section through the flow is not useful to plot.

Figure 7: Streamlines of the equatorial (left) and meridional ($\phi =0$, middle) sections through the poloidal flow field s=8, t=8, and the meridional section through the poloidal flow field s=8, t=0(right). The counterclockwise flows are marked with solid lines, the clockwise flows with dashed lines.

Figure 8: The investigated inversions for poloidal flow fields s=8 with azimuthal order t=8 (sectoral cells, left columns) and t=0 (zonal cells, right columns). The respective characteristic parameters, $u_{\rm char}$, are given in the middle column. The inversion is sensitive to the sectoral cells (strong dependence on the velocity amplitude), but insensitive to zonal cells.

3.2 Inversion for angular velocity

After generating the frequencies we use these numerically determined frequencies as helioseismic data and invert them for differential rotation. For this purpose the frequencies within a multiplet, i.e., frequencies of oscillations with identical n and l, are parameterized in terms of the central frequency $\omega_{nl}$and less than 2l+1 coefficients a_s such that

$$\omega_{nl}(m)=\omega_{nl} + \sum_{s=1}^{s_{\rm max}} a_s(n,l) P_s^{(l)}(m) .$$ (20)

The basis functions are orthogonal polynomials, which are related to Clebsch-Gordon coefficients C_s0lm^lm (Ritzwoller & Lavely 1991) by

$$P_s^{(l)}(m)=(-1)^s \sqrt{\frac{2s+1}{2l+1}} C_{s0lm}^{lm}.$$ (21)

The Wigner 3j symbols that we used in Eqs. (11) and (17) are related to the Clebsch-Gordon coefficients by

$$\left(\begin{array}{ccc} {l_1}&{l_2}&{l_3}\\ {m_1}&{m_2}&{m_... ..._1-l_2+m_3}}{(2 l_3 + 1)^{1/2}}C_{l_1 m_1 l_2 m_2}^{l_3 m_3} .$$ (22)

We are interested in the rotation rate. This is independent of the central frequency $\omega_{nl}$ and hence it does not matter whether our numerical model is able to reproduce frequencies that agree with the observed frequencies. This would be important if we were interested in a model of the Sun or the excitation of acoustic modes. The coefficients are combinations of the splittings and independent of $\omega_{nl}$.

In mathematical terms the rotation rate $\bar{\Omega}(r_0,\theta_0)$at a given radius r₀ and a colatitude $\theta_0$ can be written as

$$\bar{\Omega}(r_0,\theta_0)=\sum_{i=1}^M c_i(r_0,\theta_0) d_i ,$$ (23)

where d_i are the M data values and where we assume that they are linear functionals of the true underlying rotation rate $\Omega (r,\theta )$ only

$$d_i=\int\limits_0^R\int\limits_0^\pi K_i(r,\theta) \Omega(r,\theta)~ {\rm d}r {\rm d}\theta .$$ (24)

The kernels K_i are presumed to be known functions. The c_iare appropriately chosen coefficients. Equations (23) and (24) yield

  

$$\bar{\Omega}(r_0,\theta_0) \int\limits_0^R\int\limits_0^\pi \sum_{i=1}^M c_i(r_0,\theta_0) K_i(r,\theta) \Omega(r,\theta)~ {\rm d}r {\rm d}\theta$$ = (25)

$$\equiv \int\limits_0^R\int\limits_0^\pi {\cal K}(r_0,\theta_0;r,\theta) \Omega(r,\theta)~ {\rm d}r {\rm d}\theta.$$ (26)

The "averaging kernel'' ${\cal K}$ for a given location $(r_0, \theta_0)$describes the weighting of the true rotation rate to give an estimate of the inversion solution for that point. A full set of such averaging kernels would therefore cover the solar interior completely and describe the inverted rotation rate by an approximation $\bar{\Omega}(r,\theta)$.

Several inversion techniques have been developed. There exist two major families, the one is known as "Regularized Least Squares'' or RLS (Tikhonov & Arsenin 1977; Craig & Brown 1986), the other as "Optimally Localized Average'' (OLA; Backus & Gilbert 1968). The inversion of Eq. (13) was carried out with a two-dimensional RLS fitting, which is described by Schou et al. (1998).