Solar Interior Rotation and its Variation

Abstract

This article surveys the development of observational understanding of the interior rotation of the Sun and its temporal variation over approximately forty years, starting with the 1960s attempts to determine the solar core rotation from oblateness and proceeding through the development of helioseismology to the detailed modern picture of the internal rotation deduced from continuous helioseismic observations during solar cycle 23. After introducing some basic helioseismic concepts, it covers, in turn, the rotation of the core and radiative interior, the “tachocline” shear layer at the base of the convection zone, the differential rotation in the convection zone, the near-surface shear, the pattern of migrating zonal flows known as the torsional oscillation, and the possible temporal variations at the bottom of the convection zone. For each area, the article also briefly explores the relationship between observations and models.

1 Introduction

The internal rotation of the Sun is intimately related to the processes that drive the activity cycle. Brown et al. (1989) stated that, “Knowledge of the internal rotation of the Sun with latitude, radius, and time is essential for a complete understanding of the evolution and the present properties of the Sun,” and this remains true today.

The Sun rotates on its axis approximately once every 27 days; however, the rotation is not uniform, being substantially slower near the poles than at the equator. This superficial aspect of the solar differential rotation was well known from sunspot observations as early as the 17th century. However it is only within the last 30 years that it has become possible to observe the rotation profile in the solar interior, and mostly within the most recent solar cycle that its subtle temporal variations have become evident. Helioseismology — the study of the waves that propagate within the Sun and the inference from their properties of the solar interior structure and dynamics — is the most important tool we have to measure this internal rotation.

In this review, we start by introducing some of the basic concepts of helioseismology (Section 2) and the inversion problem (Section 3) as it applies to the internal solar rotation. Next, after a brief historical overview (Section 4) of the observations, we consider what we have learned from helioseismology about the rotation profile and its variation with depth.

We consider first the time-invariant part of the solar rotation profile. The main features of interest are (Figure 1):

• 1. the radiative interior and core, which appear to rotate approximately as a solid body, though the innermost core may behave differently (Section 5);

• 2. the tachocline, a relatively thin zone of shear between the differentially rotating convection zone and the radiative interior, which is believed to play an important role in the solar dynamo (Section 6);

• 3. the differential rotation in the bulk of the convection zone (Section 7); and

• 4. the subsurface shear layer between the fastest-rotating layer at about 0.95 R_⊙ and the surface (Section 8).

Figure 1:

A section through the interior of the Sun, showing the contours of constant rotation and the major features of the rotation profile, for a temporal average over about 12 years of MDI data. The cross-hatched areas indicate the regions in which it is difficult or impossible to obtain reliable inversion results with the available data.

We will consider each of these in turn, working outwards from the core to the surface, and then discuss the time-varying part of the rotation — the torsional oscillation (Section 9) and the possible variations at the base of the convection zone (Section 10). We attempt to place the observations in the context of models; however, this is a review from an observer’s point of view, and an exhaustive examination of the models themselves is beyond its scope.

2 Acoustic Modes

2.1 Introduction

The raw data of helioseismology consist of measurements of the photospheric Doppler velocity — or in some cases intensity in a particular wavelength band — taken at a cadence of about one minute and generally collected with as little interruption as possible over periods of months or years; the measurements can be either imaged or integrated (“Sun as a Star”). An overview of the observation techniques can be found in Hill et al. (1991a). Figure 2 shows a typical single Doppler velocity image of the Sun, and Figure 3 a portion of an l = 0 time series, derived by averaging the velocity over the visible disk for each successive image in a set of observations. The five-minute period and the rich beat structure are clearly visible in the time series. For an example of an integrated-sunlight spectrum from a long series of observations, see Figure 15.

Figure 2:

A single Doppler velocity image of the Sun from one GONG [Global Oscillation Network Group] instrument (left), and the difference between that image and one taken a minute earlier (right), with red corresponding to motion away from, and blue to motion towards, the observer. The shading across the first image comes from the solar rotation.

Figure 3:

A segment of an l = 0 time series of Doppler velocity observations, showing the dominant five-minute period and the rich beat structure.

As was first discovered by Deubner (1975), the velocity or intensity variations at the solar surface have a spectrum in k − ω or l − ν space that reveals their origin in acoustic modes propagating in a cavity bounded above by the solar surface and below by the wavelength-dependent depth at which the waves are refracted back towards the surface. These “p modes” can be classified by their radial order n, spherical harmonic degree l, and azimuthal order m; as discussed, for example, in Section 2.2 of Birch and Gizon (2005), the radial displacement of a fluid element at time t, latitude θ and longitude ϕ can be written in the form

$$\delta r(r,\theta, \phi, t) = \sum\limits_{m = -l}^l {{a_{nlm}}} {\xi _n}l(r)Y_l^m(\theta, \phi){e^{i{\omega _{nlm}}t}},$$

(1)

where ξ_nlm is the radial eigenfunction of the mode with frequency ω_nlm and $Y_{l}^{m}(\theta,\phi)$ is a spherical harmonic. As seen in Figure 4, the power in the spectrum falls along distinct “ridges” in the l − ν plane, each ridge corresponding to one radial order. The modes making up the n = 0 ridge are the so-called f modes, which are surface gravity waves. The p modes, so called because their restoring force is pressure, are excited at the surface and have their largest amplitudes there. Another class of modes, the g modes with gravity as the restoring force, excited in the core and with amplitudes vanishing at the surface, are hypothesized to exist but have so far not been definitely observed (Section 5.9).

Figure 4:

Typical l − ν spectrum from one day of GONG observations (image courtesy NSO/GONG).

The longer the horizontal wavelength — and the lower the degree — the more deeply the modes penetrate, with the radial l = 0 mode going all the way to the core of the Sun (but providing no rotational information), while modes with l ≥ 200 or so penetrate only a few megameters below the surface and are generally too short-lived to form global standing waves; these are the modes used for local helioseismology. The lower turning point radius, r_t, is a monotonic function of the phase speed ν/L, where $L = \sqrt {l(l + 1)} \approx l + 1/2$, as shown in Figure 5. The varying penetration depth with degree, as illustrated in Figure 6, makes it possible to deduce the rotation and other properties of the solar interior profile as a function of depth.

Figure 5:

Lower turning point of acoustic modes as a function of phase speed ν/L, based on Model S of Christensen-Dalsgaard et al. (1996).

Figure 6:

Locations of modes in the l, ν plane for a typical MDI mode set. The colored shading represents the radial regions in which the modes have their lower turning points; the core, r ≤ 0.2 R_⊙, the radiative interior, 0.2 ≤ r/R_⊙ ≤ 0.65, the tachocline, 0.65 ≤ r/R_⊙ ≤ 0.75, the bulk of the convection zone, 0.75 ≤ r/R_⊙ ≤ 0.95, and the near-surface shear layer, r/R_⊙ ≥ 0.95; the dashed line on the lower right corresponds to r/R_⊙ = 0.99.

The lowest-degree modes are observed in integrated sunlight, but for the purposes of measuring the interior rotation profile we are mostly concerned with what are termed medium-degree (l ≤ 300) modes, which can be observed with imaging instruments of relatively modest (≈ 10 arcsec) resolution. The power in the modes peaks at about 3 mHz, or a period of 5 minutes; useful measurements can be made for modes between about 1.5 and 5 mHz, with the frequency determination becoming more challenging at the extremes due to signal-to-noise issues and, at the high-frequency end, to the increasing breadth of the peaks.

2.2 Differential rotation and rotational splitting

The Sun’s rotation lifts the degeneracy between modes of the same l and different m, resulting in “rotational splitting” of the frequencies as waves propagating with and against the direction of rotation (prograde and retrograde) have higher and lower frequencies. To first order, the splitting δν_m,l ≡ ν_−m,l − ν_+m,l is proportional to the rotation rate multiplied by m.

Figure 7 shows a typical m − ν acoustic spectrum of GONG data at l = 100. The effect of the rotation causes the ridges at each n to slant away from the ν = 0 axis; closer examination reveals that the ridges have an S-curve shape that arises from the differential rotation, and also shows the ridge structure due to leakage, which will be discussed below in Section 2.3.

Figure 7:

Spectrum for l = 100 in the ν, m plane (top) and detail (bottom) of a single ridge (radial order) showing the curvature due to differential rotation and the multiple-ridge structure arising from spherical harmonic leakage.

Because modes of different m values sample different latitude ranges, with the sectoral (∣m∣ = l) modes confined to a belt around the equator and the zonal or m = 0 modes reaching to the poles, as illustrated in Figure 8, we can measure the rotation as a function of latitude.

Figure 8:

Spherical harmonic patterns for l = 50: left, m = 0; center, m = 45; right, m = 50.

A given (n, l) multiplet consists of 2l + 1 frequency measurements if each (l, m) spectrum is analyzed separately, though some fraction of these frequencies may be missing in any given data set. This amount of data was somewhat unwieldy in the early days of helioseismology. It is therefore common to express ν_nl(m) as a polynomial expansion, for example,

$${\nu _{nlm}} = {\nu _{nl}} + \sum\limits_{j = 1}^{{j_{\max}}} {{a_j}} (n,l){\cal P}_j^{(l)}(m),$$

(2)

where the basis functions are polynomials related to the Clebsch-Gordan coefficients $C_{j0lm}^{lm}$ by

$${\cal P}_j^{(l)}(m) = {{l\sqrt {(2l - j)!(2l + j + 1)!}} \over {(2l)!\sqrt {2l + 1}}}C_{j0lm}^{lm}$$

(3)

(Ritzwoller and Lavely, 1991). Indeed, in many analysis schemes coefficients of the expansion are derived by fitting directly to the acoustic spectrum and the individual frequencies are not measured. This approach can improve the stability of the fits, perhaps at the cost of imposing systematic errors. Early work used Legendre polynomials; however, most modern work uses either Clebsch-Gordan coefficients or the Ritzwoller-Lavely formulation, which come closer to being truly orthogonal for the solar rotation problem. Only the odd-order coefficients encode the rotational asymmetry, while the even-order coefficients contain information about the structural asphericity. Roughly speaking, the a₁ coefficient describes the rotation rate averaged over all latitudes, and the a₃ and higher coefficients describe the differential rotation.

2.3 Spherical harmonics and leakage

Spherical harmonic masks are used to separate the contributions from modes of different degree and azimuthal order into complex time series, which can then be transformed to acoustic Fourier spectra.

The radial component of the velocity at the solar surface from a mode with a given degree l, azimuthal order m, and radial order n is given by

$${V_{nlm}}(\phi, \theta, t) = {\mathop{\rm Re}\nolimits} [{a_{nlm}}(t)P_l^{\vert m\vert}(\cos \theta){e^{im\phi}})],$$

(4)

where Re[…] denotes the real part, ϕ is longitude, and θ is latitude (see, for example, Schou and Brown 1994a). The masks used separate the different spherical harmonics take the form

$${M_{l,m}} \propto {Y_{l,m}}(\theta, \phi)A(\rho),$$

(5)

where A is an apodization function and $\rho \equiv \sqrt {{{\cos}^2}\theta + {{\sin}^2}\theta {{\sin}^2}\phi}$ represents the fractional distance from disk center in the solar image. The line-of-sight projection factor is $V = \sqrt {1 - {\rho ^2}}$.

Because only part of the solar surface is visible at any time, the masks are not completely orthogonal and the modes “leak” into neighboring spectra. This leakage complicates the analysis and can cause systematic errors in the measured frequencies if it is not correctly taken into account. For a detailed discussion of the calculation of the so-called “leakage matrix,” see Schou and Brown (1994a) and Hill and Howe (1998). Briefly, the leakage matrix element s(l, m, l′, m′) for leakage from the l′, m′ mode to the l, m spectrum can be computed as

$$s(l,m,{l^{\prime}},{m^{\prime}}) = {1 \over \pi}\int_{- 1}^1 {\int_{- \pi/2}^{\pi/2} {P_l^m}} (x)P_{{l^{\prime}}}^{{m^{\prime}}}(x)\cos (m\phi)\cos ({m^{\prime}}\phi)V(\rho)A(\rho)dxd\phi.$$

(6)

Symmetry properties in this expression lead to some simple exclusion rules; for example, leaks with odd ∣δl + δm∣ (where δm ≡ m − m′ and δl ≡ l − l′) are not allowed.

One example of the importance of the leakage is in the contribution of the so-called m-leaks (δl = 0, δm = ±2) to the estimation of low-degree splittings. As pointed out, for example, by Howe and Thompson (1998), these leaks are strongest for small ∣m∣; they are also asymmetrical, especially for ∣m∣ = l, where the m = l peak has an m = l − 2 leak on one side and no counterbalancing m = l + 2 leak on the other. Especially for l = 1, this can introduce a serious systematic error into the estimate of the splitting if not properly accounted for.

Leakage also means that integrated-sunlight instruments (which effectively use the l = 0 mask) can detect modes of 0 ≤ l ≤ 5, though the sensitivity falls off rapidly for l > 1. All these modes appear in a single acoustic spectrum; the instruments are sensitive to odd-m modes for odd l and to even-m modes for even l, with the sectoral, or ∣m∣ = l, modes most strongly detected.

In general, the leakage has effects throughout the acoustic spectrum, but the most deleterious effects arise when the leaks cannot be resolved from the target peaks. This occurs for m-leaks at frequencies above about 2 mHz; for higher-degree modes the leakage between modes of adjacent l becomes a problem, as the ridges become both broader, and more closely spaced in frequency, at around l = 150. Beyond this point the peaks cannot be fitted independently, and some modeling of the leakage is essential in order to estimate the mode parameters.

2.4 Estimating rotation properties directly from coefficients

It is possible to make some inferences about the rotation profile without carrying out a full-scale inversion. Simple examination of the odd-order coefficients, sorted by the lower turning-point radius of the modes, reveals the existence of the near-surface shear, the differential rotation within the convection zone, and a discontinuity in the differential rotation at the base of the convection zone, as shown in Figure 9. More sophisticated analysis is also possible. For example, Wilson and Burtonclay (1995) gave approximate expressions for the rotation profile at different latitudes as sensed by a particular n, l multiplet, ${\bar \Omega ^{nl}}$, as follows:

$$\bar \Omega _0^{nl} \approx a_1^{nl} + a_3^{nl} + a_5^{nl},$$

(7)

$$\bar \Omega _{30}^{nl} \approx a_1^{nl} - {{a_3^{nl}} \over 4} - {{19a_5^{nl}} \over {16}},$$

(8)

$$\bar \Omega _{45}^{nl} \approx a_1^{nl} - {{3a_3^{nl}} \over 2} - {{3a_5^{nl}} \over 4},$$

(9)

$$\bar \Omega _{60}^{nl} \approx a_1^{nl} - {{11a_3^{nl}} \over 4} + {{37a_5^{nl}} \over {16}},$$

(10)

These estimates, where the subscripts on the LHS refer to the latitude in degrees, are noisy for individual multiplets, but Wilson and Burtonclay (1995) were able to build up a picture of the internal rotation from BBSO data by forming cumulative averages with the input data sorted in ascending order of ν/L.

Figure 9:

a₁ (top) and a₃ (bottom) coefficients for (left) 1986 BBSO observations, (middle) 108 days of GONG observations in 1996, (right) the mean of 35 consecutive 108-day periods of GONG observations from 1995–2005, plotted as a function of phase speed with the turning point radius marked on the upper axis.

3 Inversion Basics

Various inversion techniques are used to infer the internal rotation profile from the observed frequency splittings. The aim of the inversion procedure is to form linear combinations of the data that give well-localized inferences of the rotation at a particular location within the Sun. We will discuss only linear inversion methods, as non-linear approaches are not needed for the relatively low velocities involved in the global rotation.

3.1 The inversion problem

The basic 2-dimensional rotation inversion problem can be stated as follows: we have a number M of observations d_i, from which we wish to infer the rotation profile Ω(r, θ) where r is distance from the center of the Sun, and θ is (conventionally) colatitude. Each datum is a spatially weighted average of the rotation rate:

$${d_i} = \int_0^{{R_ \odot}} {\int_0^\pi {{K_i}}} (r,\theta)\Omega (r,\theta)drd\theta + {\epsilon _i},$$

(11)

where R_⊙ is the solar radius, the error term ε corresponds to the noise and measurement error in the data, and K is a model-dependent spatial weighting function known as the kernel (Hansen et al., 1977; Cuypers, 1980). For the two-dimensional rotation inversion, the radial part is related to the eigenfunction of the mode and the latitudinal part to the associated Legendre polynomial; Schou et al. (1994) give the expression for the kernel as

$$\matrix{{{K_{nlm}}(r,\theta) = {m \over {{I_{nl}}}}\left\{{{\xi _{nl}}(r)\left[ {{\xi _{nl}}(r) - {2 \over L}{\eta _{nl}}(r)} \right]P_l^m{{(x)}^2}} \right.} \cr {\left. {+ {{{\eta _{nl}}{{(r)}^2}} \over {{L^2}}}\left[ {{{\left({{{dP_l^m} \over {dx}}} \right)}^2}{{(1 - x)}^2} - 2P_l^m{{dP_l^m} \over {dx}}x + {{{m^2}} \over {1 - {x^2}}}P_l^m{{(x)}^2}} \right]} \right\}\rho (r)r\sin \theta,} \cr} $$

(13)

where

$${I_{nl}} = \int_0^{{R_ \odot}} {[{\xi _{nl}}{{(r)}^2} + {\eta _{nl}}{{(r)}^2}]} \rho (r){r^2}dr,$$

(13)

x = cos θ, L² = l(l + 1), ξ_nl is the radial displacement for the eigenfunction of the mode, L⁻¹η_nl is the horizontal displacement, and ρ(r) is the density (see Figure 10 for illustrations of sample kernels).

Figure 10:

Sections through rotation kernels for selected azimuthal orders for l = 3, n = 9 (top), and l = 20, n = 5 (bottom).

The aim of the inversion is to find

$$\bar \Omega ({r_0},{\theta _0}) = \sum\limits_{i = 1}^M {{c_i}} ({r_0},{\theta _0}){d_i},$$

(14)

where (r₀, θ₀) is the location at which the inferred rotation rate $\bar{\Omega}$ is to be found and the c_i are the coefficients to be used to weight the data; the inversion process can be thought of as the search for the best values for these coefficients.

3.2 Averaging kernels

By substituting Equation (11) into the RHS of Equation (14) we obtain

$$\bar \Omega ({r_0},{\theta _0}) = \int_0^{{R_ \odot}} {\int_0^\pi {{\cal K}({r_0},{\theta _0};r,\theta)} \Omega (r,\theta)drd\theta + {\varepsilon _i}},$$

(15)

where

$${\cal K}({r_0},{\theta _0};r,\theta) \equiv \sum\limits_{i = 1}^M {{c_i}} ({r_0},{\theta _0}){K_i}(r,\theta)$$

(16)

is the averaging kernel for the location (r₀, θ₀). The averaging kernels are independent of the values of the data. However, because the uncertainties in the data are used to weight the inversion calculation that generates the coefficients c_i, as described below in Sections 3.4 and 3.5, these do enter indirectly into the averaging kernels. The averaging kernels also depend on which modes are present in the input data set. They provide a useful tool for assessing the reliability of an inversion inference from a particular mode set (see, for example, Schou et al., 1992, 1994).

3.3 Inversion errors

If the errors on the input data are uncorrelated and properly described by a normal distribution whose width corresponds to the quoted uncertainty σ_i, the formal uncertainty on the inferred profile is given by

$${\sigma ^2}[\Omega ({r_0},{\theta _0})] = \sum\limits_i {{{[{c_i}({r_0},{\theta _0}){\sigma _i}]}^2}}.$$

(17)

In the (usually unrealistic) case where the errors on the input data are all equal, we can write

$${\sigma ^2}[\Omega ({r_0},{\theta _0})] = \Lambda ({r_0},{\theta _0})\sigma,$$

(18)

where the “error magnification” is given by

$$\Lambda ({r_0},{\theta _0}) = {\sum\limits_i {[{c_i}{{({r_0},{\theta _0})}^2}]} ^{1/2}}.$$

(19)

As discussed, for example, by Christensen-Dalsgaard et al. (1990), a quantitative choice of regularization parameters can then be made by finding the “knee” of a tradeoff curve where the error magnification is plotted against the width of the averaging kernel. However, in the two-dimensional case this does not always give a clear result, and this formulation of the error magnification is not very useful for modern data sets where the uncertainties on the parameters are anything but uniform. Instead, one can consider the uncertainty on the inferred quantity at a particular location.

Even when the errors on the input data are uncorrelated, the errors on the inferred profile will not be, as discussed by Howe and Thompson (1996). (As a simple way to understand this, consider the case where one measurement is significantly “off”, this will affect the inferred profile at every location where the inversion coefficient c_i for that datum is non-zero.) In the one-dimensional case, the correlation between the errors for two points r₀ and r₁ is given by

$$C({r_0},{r_1}) = {{\sum {{c_i}} ({r_0}){c_i}({r_1})\sigma _i^2} \over {{{[\sum {c_i^2} ({r_0})\sigma _i^2]}^{1/2}}{{[\sum {c_i^2} ({r_1})\sigma _i^2]}^{1/2}}}};$$

(20)

this can easily be generalized to the two-dimensional case. Howe and Thompson (1996) found that the spatial scale over which the inversion errors are significantly correlated is usually similar to that for the averaging kernels, though for some cases where the inversion parameters have been badly chosen the results can be correlated over long distances even when the averaging kernels appear well formed.

Error correlations by definition should not distort the inferred profile beyond the distribution predicted by the formal uncertainty on the inferences, provided always that the input uncertainties are correct. However, the finite width of averaging kernels also gives rise to a systematic error that can be much larger. Consider, for example, the case where a thin shear layer is not resolved; then all the estimated rotation rates on one side of the shear could be underestimated, and those on the other side overestimated, by several times the formal uncertainty. Such systematic errors and their relationship to the averaging kernels have been discussed, for example, by Christensen-Dalsgaard et al. (1990).

Gough et al. (1996) pointed out that it is not sufficient for the rotation rates at two locations to have non-overlapping errors as calculated in Equation (17), and described a method for increasing the error estimates on inversions to allow truly significant differences between the inferred rotation rate at different locations to be determined. This method, however, has not been widely used.

Because the input data are noisy and of finite resolution, the inversion problem does not have a unique solution; there will always be a tradeoff between noise and good localization. Two widely-used approaches to balancing these criteria are “regularized least squares” (RLS) and “optimally localized averaging” (OLA).

3.4 Regularized least squares

The RLS approach to the inversion problem is to find (essentially through a least-squares fit) the model profile that best fits the data, subject to a smoothness penalty term, or regularization. More regularization — a larger weighting for the penalty term — results in poorer spatial resolution (and potentially more systematic error) but smaller uncertainties. In one such implementation (Schou et al., 1994), we minimize

$$\sum\limits_i {{{{{[{d_i} - \int_0^R {\int_0^\pi {\bar \Omega}} (r,\theta){K_i}(r,\theta)drd\theta ]}^2}} \over {{{({\sigma _i}/\bar \sigma)}^2}}}} + \mu _r^2\int_0^R {{{\int_0^\pi {\left({{{{d^2}\bar \Omega} \over {d{r^2}}}} \right)}}^2}} drd\theta + \mu _\theta ^2\int_0^R {{{\int_0^\pi {\left({{{{d^2}\bar \Omega} \over {d{\theta ^2}}}} \right)}}^2}} drd\theta $$

(21)

with μ_r and μ_θ being the radial and latitudinal tradeoff parameters. The RLS inversion has the advantages of being computationally inexpensive and always (thanks to the second-derivative regularization, which amounts to an a priori assumption of smoothness) providing some kind of estimate of the quantity of interest even in locations that are not, strictly speaking, resolved by the data. In this method, the averaging kernels ${\cal K}$ can (but need not be) calculated from the coefficients in a separate step. They are not guaranteed to be well localized, though they are forced to have a center of mass at the specified location r₀,θ₀. Figure 11 illustrates typical averaging kernels for a 2dRLS inversion of an MDI data set.

Figure 11:

Averaging kernels for a typical RLS inversion of MDI data, for target latitudes 0 (a), 15 (b), 30 (c), 45 (d), 60 (e), and 75 (f) degrees as marked by the dashed radial lines, and target radii 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99 R_⊙ indicated by colors from blue to red as denoted by the dashed concentric circles. Contour intervals are 5% of the local maximum value closest to the target location, with dashed contours indicating negative values.

3.5 Optimally localized averaging

In the Subtractive OLA (SOLA) approach (Backus and Gilbert, 1968, 1970), the minimization is applied to the difference between the actual averaging kernels ${\cal K}$ and a target kernel ${\cal T}$, for example a 2-dimensional Gaussian or Lorentzian function. In this case (Pijpers and Thompson, 1992, 1994) the function minimized is

$$\int_0^R {\int_0^\pi {[{\cal T}({r_0},{\theta _0};r,\theta)} - {\cal K}} ({r_0},{\theta _0};r,\theta){]^2}rdrd\theta + \lambda {\sum\limits_{i = 1}^M {[{\sigma _i}{c_i}({r_0},{\theta _0})]} ^2}.$$

(22)

Both the tradeoff parameter λ and the radial and latitudinal resolution of the inversions must be chosen before running the inversion. If the choice of target kernel is poor — too narrow or too wide for the quantity and quality of the data — the reliability of the inversion will suffer. In OLA inversions, setting target locations outside the regions that can be resolved using the data will result in averaging kernels displaced from their targets, and this should be taken into account when interpreting the results. Figure 12 illustrates typical averaging kernels for a 2d SOLA inversion of an MDI data set.

Figure 12:

As Figure 11, for a SOLA inversion.

Another approach, older, and more computationally expensive, is the Multiplicative OLA (MOLA) described by Pijpers and Thompson (1992, 1994). Here, no target form is imposed on the averaging kernel, but it is multiplied by a term which penalizes large values away from the target location.

3.6 Other methods

Alternatives to full 2-dimensional inversions are the so-called “1.5-dimensional” approach, in which 1-dimensional radial inversions are carried out separately for each of the coefficients describing the latitudinal rotation variation, and “1d⊗1d” inversions in which the radial and latitudinal variations of the rotation rate are integrated separately. For details of many of these methods, please see Schou et al. (1998) and references therein.

3.7 Limitations

It is important to bear in mind the limitations of the inversion process when considering the results. The deepest and shallowest depths that can be resolved, for example, are limited by the deepest and shallowest turning-point radii of the available modes. The rotational splitting at a given m is to first order proportional to the rotation rate multiplied by ∣m∣; since the only mode whose latitudinal kernel reaches the pole is the m = 0 mode, which has no longitudinal structure and so can convey no rotational information, and the modes of small ∣m∣/l have only a few nodes around the equator and hence have low sensitivity to the rotation, the 2d inversion becomes progressively less reliable at high latitudes. Furthermore, since only modes of relatively low degree (l ≤ 20) penetrate into the radiative interior, the latitudinal resolution in this region is quite poor and becomes progressively worse with depth; radial resolution also becomes coarser in the interior. The practical effects of such limitations can be assessed by careful inspection of the averaging kernels, or by performing forward-calculation tests in which the averaging kernels are convolved with known test profiles.

Another point to bear in mind when considering inversion results is that the inversion can measure only the north — south symmetric part of the profile; any asymmetry between the hemispheres is averaged out. The inversions are also insensitive to meridional motions. Some information on hemispheric differences can be obtained using the techniques of local helioseismology, as reviewed by Birch and Gizon (2005), but these techniques, using high-degree modes, are mostly sensitive only to the outer layers of the Sun.

4 Observations: A Brief Historical Overview

Systematic helioseismic observations stretch back nearly 30 years, as illustrated in the schematic chart in Figure 13.

Figure 13:

Schematic time line of helioseismic observations in the last three solar cycles (top panel), with the filled part of each bar representing approximate duty cycle, plotted on the same temporal scale as the butterfly diagram (bottom panel) of the gross magnetic field strength from Kitt Peak observations.

Prior to the identification of global low-degree modes by Claverie et al. (1979), observing runs were usually short and carried out at a single site. However, the advantages of more extended observations (to obtain better frequency resolution), and of observations not modulated by the day-night cycle, were soon recognized. Grec et al. (1980), Duvall Jr and Harvey (1984), and Duvall Jr et al. (1984, 1986) carried out important observations at the South Pole during the Austral summer, but for long time series it is more practical to observe either from a network of sites spaced around the world, or from space.

Some of the first long-term sets of low-degree observations came from the Active Cavity Irradiance Monitor (ACRIM) experiment (Woodard and Noyes, 1985, 1988) aboard the Solar Maximum Mission spacecraft, which took helioseismic measurements in 1980 and 1984–1985, the Mark I instrument in Tenerife (Pallé et al., 1989), and the precursors of the Birmingham Integrated Solar Network (BiSON) (Elsworth et al., 1990a). Meanwhile, resolved-Sun observations were carried out at the South Pole by Duvall and collaborators, and by various other observers in the USA; these observations will be discussed in more detail later.

Libbrecht and Woodard (1990) observed the medium-degree modes from Big Bear Solar Observatory (BBSO) in the 1986–1990 rising phase of solar cycle 22. The first observations from widely separated sites were carried out by the Birmingham/Tenerife group in 1981 (Claverie et al., 1984), and by 1992 the six-station BiSON network was complete; it has been operating ever since. Another network of integrated-sunlight instruments, the French-based IRIS (Fossat, 1995), operated from 1989–2003.

The Global Oscillation Network Group (GONG) (Harvey et al., 1996) has been collecting continuous, high-duty-cycle observations of the medium-degree p modes since 1995, using a six-station worldwide network, and the Michelson Doppler Imager (MDI) instrument (Scherrer et al., 1995) aboard the SOHO spacecraft has been in operation since 1996, so that these two projects have essentially complete coverage of solar cycle 23. SOHO also carries instruments dedicated to the study of low-degree oscillations; LOI (Luminosity Oscillations Imager) (Fröhlich et al., 1995), and GOLF (Global Oscillations at Low Frequencies) (Gabriel et al., 1995). This wealth of high-quality data has given us the opportunity to study the solar interior rotation and its solar-cycle changes in more detail than ever before.

Also worth noting are the LOWL-ECHO project (Tomczyk et al., 1993), which made medium-degree observations from one or two sites from 1994 to 2004, and the high-degree Taiwanese Oscillations Network (Chou et al., 1995) deployed over the 1993–1996 period.

All these observations will be considered in more detail as we proceed to examine the results pertaining to the interior rotation.

5 The Core and Radiative Interior

5.1 The oblateness controversy

Interest in the rotation of the deep solar interior predates systematic helioseismic observation. One other possible diagnostic of the internal rotation is provided by the solar oblateness; because the Sun is not a solid body, both gravitational and rotational effects cause a very slight flattening. The lowest-order term in this effect is related to the quadrupole moment J₂; confusingly, the next-highest term, J₄, is sometimes called octopole and sometimes hexadecapole. According to Rozelot and Rösch (1997), who give a useful review of attempts to measure the solar oblateness, for a non-rotating Sun the oblateness Δr = r_eq − r_pol is given by

$${{\Delta r} \over {{r_0}}} = {3 \over 2}{J_2},$$

(23)

where r_eq and r_pol are the equatorial and polar radii, respectively, and r₀ is the radius of the best sphere passing through r_eq and r_pol. If there is an additional δr contribution from the surface rotation this expression becomes

$${{\Delta r - \delta r} \over {{r_0}}} = {3 \over 2}{J_2}.$$

(24)

The units of δr and Δr are conventionally arc ms.

Dicke (1964) noted that, if the Sun were oblate because of fast interior rotation, the effect on its gravitational potential might destroy the agreement between the predictions of General Relativity and the observations of the perihelia of the inner planets, (specifically Mercury, though Venus could in principle experience a smaller effect) potentially leaving room for alternate theories of gravitation. Dicke sets out to determine the solar oblateness from ground-based measurements — a challenging endeavor that produced controversial results. Models (e.g., Brandt 1966) suggested that the interior of the Sun could still be spinning at the rapid rate at which it originally formed, while the exterior had been slowed down by the torque of the solar wind. (As will be further discussed in Section 6, in the absence of direct observations of the solar interior the picture of solar interior dynamics was not at all clear, although the existence of something like what we now call the tachocline could be inferred from theoretical arguments.) Dicke and Goldenberg (1967b) reported finding a solar oblateness value of 5 × 10⁻⁵, which would be sufficient to create an 8% discrepancy between observations and the Einsteinian prediction for the precession of the perihelion of Mercury, and would imply a fast-rotating core.

The results, and the inferences Dicke and collaborators drew from them, raised a storm of controversy that may well have helped to stimulate interest in the Sun’s interior rotation profile. The criticisms and Dicke’s responses to them would fill a lengthy review article by themselves; we give only a few examples here. Roxburgh (1967) suggested that the result might be explained by the solar differential rotation, an idea rejected by Dicke and Goldenberg (1967a). Howard et al. (1967) concluded, on the basis of a variety of simple models of the solar “spin-down,” that the Sun should have reached a state of uniform rotation quite quickly after its initial formation. Sturrock and Gilvarry (1967) pointed out that the presence of magnetic field in the solar interior might well complicate the issue, and in an accompanying article Gilvarry and Sturrock (1967) suggested using a space probe in a highly eccentric orbit as a more direct test of general relativity — or, alternatively, that “more complete theoretical and observational knowledge of the visible layers and the interior of the Sun” was needed.

At least partly inspired by the controversy, Kraft (1967) studied the rotational velocities of young solar-type stars in the Pleiades and concluded that angular momentum was lost on a timescale of about half a billion years, but noted in his conclusion that “it is wrong to conclude that the present work in any way supports the Dicke result.” Goldreich and Schubert (1968) considered the stability of differentially rotating stars and concluded that it was possible but not likely that a radial rotation gradient such as that required by the Dicke and Goldenberg (1967b) result might exist.

H. Hill, a former colleague of Dicke who had helped to build the instrument with which the 1964 observations were made (Dicke, 1964), and collaborators, also attempted to measure the solar oblateness, using an instrument, SCLERA [Santa Catalina Observatory for Experimental Relativity by Astrometry], which was later to play a role in the early days of helioseismology. This measurement, carried out in 1973, (Hill and Stebbins, 1975), found a 9.6 × 10⁻⁶ value for the oblateness, much smaller than that of Dicke and Goldenberg (1967b). Hill et al. (1974) also pointed out a time-varying difference between the brightness of the solar limb and poles that might account for the anomalously high oblateness measurement.

Ulrich and Hawkins (1981a, b) made an early attempt to deduce what the J₂ and J₄ terms should be based on a simple differential rotation profile deduced from surface measurements, obtaining predicted values of between 1 and 1.5 × 10⁻⁷ for J₂ and between 2 and 5 × 10⁻⁹ for J₄ depending on the size of the convective envelope.

Dicke et al. (1986, 1987) repeated the 1966 measurements with an improved instrument, and obtained significantly smaller values for the oblateness, with some weak evidence for a solar-cycle variation. Lydon and Sofia (1996) made measurements using a balloon-based instrument and obtained values of 1.8 × 10⁻⁷ for J₂ and 9.8 × 10⁻⁷ for J₄. By this point, however, the focus in the solar oblateness studies had moved away from trying to infer the core rotation. Mecheri et al. (2004) used more realistic models of the internal rotation profile to suggest that the J₄ term should be particularly sensitive to the subsurface shear. Recent work on determining the oblateness from the shape of the solar limb has taken into account considerations of near-surface temperature or magnetic variations. Kuhn et al. (1998) and Emilio et al. (2007) used observations from MDI during rolls of the SOHO spacecraft and Fivian et al. (2008) used the RHESSI X-ray telescope. The work with SOHO revealed a temporal variation in the shape of the solar limb, with greater apparent oblateness at solar maximum, suggesting that hotter, brighter activity belts have greater apparent diameter. This poses an apparent contradiction to the results obtained from helioseismic inferences of the asphericity. Indeed, Fivian et al. (2008) suggest that all the temporally-varying, excess oblateness found in the observations can be corrected away by removing an ad-hoc term related to magnetic elements in the enhanced network.

Meanwhile, a much more flexible tool — helioseismology — had become available for probing the interior solar rotation.

5.2 Early low-degree helioseismic results

Around the early 1970s there were numerous attempts to search for global p-mode oscillations, with interest at first focusing on longer-period oscillations, the low-order, low-degree modes. Various theoretical predictions (Scuflaire et al., 1975; Iben Jr, 1976; Christensen-Dalsgaard and Gough, 1976) of the periods were available, offering the hope that global oscillations could be used to probe the rotation and structure deep inside the Sun. At first most of the results (Livingston et al., 1977; Musman and Nye, 1977; Grec and Fossat, 1977), were negative, except for the 160-minute period of Severnyi et al. (1976); Brookes et al. (1976), which was later (Elsworth et al., 1989) determined to be spurious and will not be further discussed here. The SCLERA group (Brown et al., 1978; Hill and Caudell, 1979; Caudell and Hill, 1980) found a variety of longer-period fluctuations in their solar-diameter data, but these results were not universally accepted; for example Fossat et al. (1981a, see also references therein) claimed that the SCLERA results were consistent with pure noise.

Low-degree helioseismology became a reality when the Birmingham group (Claverie et al., 1979) identified oscillations in the five-minute frequency band in integrated sunlight as low-degree global modes, using observations from Tenerife and Pic du Midi during the summers of 1976–1978; these initial data were adequate only to identify the spacing between modes of the same l and different n, without resolving separate l = 0 and l = 1 peaks.

A French-American team (Grec et al., 1980; Fossat et al., 1981b) obtained five days of continuous observations at the South Pole in the austral summer of 1979–1980, and were able to identify peaks of degree 0, 1, 2, and 3 and even a weak l = 4 peak by superposing sections of the acoustic spectrum with different radial order. These modes were identified as being of radial order around 12–30, as opposed to the very low-order modes that had been sought in the low-frequency spectrum; both the noise characteristics of the spectrum and the low amplitude of the lower-order modes mean that the fundamental (l = 0, n = 1) mode remains unobserved to this day, although some low-degree modes with single-digit n have been identified (Chaplin et al., 1996b).

Soon, the Birmingham team (Claverie et al., 1981), using 28 days of integrated-sunlight data from the Tenerife site and an analysis that involved “collapsing” segments of the acoustic spectrum so as to average together modes of the same degree and different radial order, reported finding three rotationally split components in the l = 1 modes and five in l = 2, with an average separation of 0.75 µHz. If correct, this would have implied a solar core rotation substantially faster than the surface. Isaak (1982) suggested that the excess component peaks (when two and three would be the expected number for l = 1 and l = 2 respectively) could be explained if the solar core were rotating on an oblique axis and had a very strong magnetic field; this idea, which was also mooted by Dicke (1983) to explain an oscillation of about half the solar rotation period seen in the oblateness data (Dicke, 1976), was rebutted in some detail by Gough (1982).

Fossat et al. (1981b) reported that initial results from 5 days of low-degree observations at the South Pole suggested quite short lifetimes, about 2 days; the l = 0 peaks appeared narrower than those of l = 1 and l = 2. Grec et al. (1983) later identified about 80 normal modes in the South Pole data, but did not confirm the Claverie et al. (1981) rotational splitting result, instead reporting that the l = 1 peak seemed too narrow to accommodate the reported splitting.

Claverie et al. (1982) reported a periodicity of approximately 13 days in the radial solar velocity, as measured using the resonant-scattering technique and the potassium D-line, and interpreted this as an effect of the solar core rotation; however, this effect was quickly explained away (Durrant and Schröter, 1983; Andersen and Maltby, 1983; Edmunds and Gough, 1983; Duvall Jr et al., 1983) as an artifact caused by the rotation of surface features — sunspots and plage — across the disk.

Meanwhile, the low-degree five-minute acoustic spectrum had also been observed using the Active Cavity Irradiance Monitor (ACRIM) aboard the Solar Maximum Mission (SMM) spacecraft (Woodard and Hudson, 1983b). Woodard and Hudson (1983a) agreed with Fossat et al. (1981b) in finding that the modes had lifetimes of about two days, too short for the rotational splitting reported by Claverie et al. (1981) to be real.

Later work (Libbrecht, 1988a; Elsworth et al., 1990b; Chaplin et al., 1997) revealed that the width of the peaks — inversely proportional to the mode lifetimes — was strongly dependent on frequency across the five-minute spectrum, with lifetimes of a few days in the middle of the five-minute band and weeks or months at low frequencies where, unfortunately, the amplitudes of the modes are also small. Reliable direct measurement of the low-degree splittings would have to wait for some years, while sufficiently long, high-quality time series of data accumulated.

5.3 Resolved-Sun measurements

In the meantime, resolved-Sun observations provided some information about the rotation in the radiative interior. Duvall Jr and Harvey (1984) reported observations at Kitt Peak, from 10–26 May 1983, for degrees 0 ≤ l ≤ 100. When plotted as a function of degree, the results show a slow decrease in the rotational splitting, from the highest degrees down to about l = 6, with an unexplained bump at l = 11, followed by an increase at lower degrees up to a value of 660 nHz for l = 1. These data, inverted by Duvall Jr et al. (1984), yielded a rotation profile with much of the radiative interior rotating at or below the surface rate, but with a modest increase in the interior. A similar pattern was found by Brown (1985), using 6 days of observations from the newly-developed Fourier Tachometer, a true 2-dimensional imaging instrument that gave access to all the azimuthal orders for degrees between 8 and 50; however, the coincidence of the l = 11 bump seems to have been merely a coincidence of noise, as it was not reproduced in the early observations from the Big Bear Solar Observatory (Libbrecht, 1986).

5.4 The SCLERA modes

Hill et al. (1982) derived splittings from the SCLERA low-frequency peaks, and from those inferred a core rotating at 6 times the surface rate; however, Woodard (1984) used ACRIM data to place an upper limit of 2.2 times the surface rate on the interior rotation rate, inconsistent with these splittings. Later, Hill (1985) identified low-degree rotational splittings in the five-minute band of the SCLERA acoustic spectrum, but Libbrecht (1986) and Brown et al. (1989) found that these results were inconsistent with the other evidence and were probably the result of misidentification of the modes. Given the complexity of the spectrum in question, whose derivation from measurements sampled at a few points on the solar limb made it difficult to separate out spectra of different degree, this seems a likely explanation.

5.5 Low-degree acoustic mode splittings 1988–2002

The next several years were active ones for low-degree helioseismology, with the development of the BiSON (Birmingham-based) and IRIS (based in Nice) networks. Together with the IPHIR instrument that rode the PHOBOS spacecraft on its cruise phase to Mars, and the ground tests of the LOI (Luminosity Oscillations Imager) instrument that would later be mounted on the SOHO spacecraft, these brought a succession of estimates of the low-degree splitting, as summarized in Table 1 and Figure 14. In addition to the MDI instrument for medium and high-degree observations, the SOHO spacecraft carried both LOI and GOLF (Global Oscillations at Low Frequencies) specifically for observing low-degree modes. Even though GOLF malfunctioned and could not be operated in its intended differential mode, instead being confined to making Doppler observations on one side of an absorption line, it provided some of the best available long-term, low-degree observations.

Figure 14:

l = 1 splitting estimates as a function of publication date.

Table 1: Summary of l = 1 splitting measurements, 1988–2002

The reported results show considerable variation, but apart from the early Tenerife result, which was based on much shorter and lower-duty-cycle observations than most of the others, they all cluster around the surface rotation rate, some (particularly the IRIS results) pointing to a core rotation faster than the surface rate and some (in particular the BiSON results) to one substantially below it, perhaps as low as zero. As we approach the present time and the observation and analysis improve, the values tend to converge on a splitting quite close to that which would correspond to the surface rate. Early in this period, there was room to speculate (e.g., Chaplin et al., 1996a) that the differences reflected a temporal variation, but this could not explain away all the discrepancies.

5.6 Pitfalls of low-degree splitting measurements

Unfortunately, all the measurements described in Section 5.5 suffer from similar problems, as summarized below:

• 1. The two components of the l = 1 mode are so close together (probably less than one microhertz apart) that they are resolved only for modes below about 2.2 mHz. This has implications for the measurements:

  • (a) Estimates of the splittings of unresolved components are highly prone to systematic error (Appourchaux et al., 2000a).

  • (b) The components that can be resolved have small amplitudes (Figure 15) and therefore require both observations over extended periods and high signal-to-noise ratios.

  • (c) On the other hand, these low-frequency modes have the advantage that they show very little frequency shift with the solar cycle, which simplifies the analysis of long time series.

• 2. Even though the low-degree modes penetrate deep into the solar interior, they spend most of their time in the outer layers of the Sun and are not very sensitive to the core; conversely, estimates of the core rotation are very sensitive to small errors in the splitting measurements.

• 3. In order to properly estimate the rotation profile in the deep interior it is necessary to combine the low-degree splittings with medium-degree ones in an inversion. However, because the low-degree modes are so few — a few dozen at most, compared to a couple of thousand medium-degree multiplets with tens of thousands of individual frequencies or coefficients — the need for extremely precise measurements is even more pressing. Also, combining data from different instruments with different systematic errors may cause problems, particularly if the observations were made at different epochs of the solar cycle.

Figure 15:

Power spectrum from 10 years of BiSON data, 1992–2002; the insets show the low-frequency end of the five-minute band (blue) and a single, rotationally split l = 1 peak (red).

Point 1 above was noted by Loudagh et al. (1993) and Elsworth et al. (1995), and point 2 by Loudagh et al. (1993) and Lazrek et al. (1996), who point out that “An accuracy of about 30 nHz, or (1 year)⁻¹ on the measurement of the l = 1 rotation splitting does not really permit, then, to discriminate between a solar core rotating twice as fast as the rest or not rotating at all!”. An approach to addressing point 3 was made by Tomczyk et al. (1995) with the newly-built LOWL instrument, an imaging instrument optimized for lower degrees. They obtained splittings for 1 ≤ l ≤ 100, and inferred a rotational profile down to 0.2 R_⊙, finding a rotation rate that barely varied with radius between 0.2 R_⊙ and 0.6 R_⊙, apart from a low-significance bump around 0.4 R_⊙.

Eff-Darwich and Korzennik (1998) further addressed point 3 when they combined results from several different instruments, including GONG, BiSON, MDI, and GOLF. They give a nice illustration of the tendency of higher-frequency low-degree mode splittings to be biased upward by the mode width, a point that was further illustrated by Chaplin et al. (2001), and conclude that with the then-available data it is not possible to rule out fast rotation in the core below 0.18R_⊙.

Charbonneau et al. (1998) used a genetic forward-modeling approach to analyze the LOWL data, with results favoring a rigidly-rotating core.

5.7 A new millennium for low-degree helioseismology

Starting around the turn of the century, there was a move towards more collaborations and comparisons between different projects in an effort to understand the systematic errors and better constrain the solar core dynamics. By this time, multi-year observations were available from GONG and the SOHO instruments, as well as good-quality observations from BiSON stretching back to 1991.

Chaplin et al. (1999) combined the LOWL higher-degree splittings with the very precise low-frequency BiSON splittings for the lowest-degree modes, and concluded that the data were consistent either with rigid rotation or with a slight downturn in the rotation rate in the core (the latter being at best a 1-σ result); on the other hand, Corbard et al. (1998b) had used a very similar analysis of GOLF and MDI data to deduce a slight increase in the rotation rate below 0.25 R_⊙, but García et al. (2003), also using MDI and GOLF data, obtained rather low splitting values from a 2243-day time series and tentatively concluded that they could rule out a high rotation rate in the core.

Eff-Darwich et al. (2002), following on from the work of Eff-Darwich and Korzennik (1998), again combined BiSON, GOLF, GONG and MDI data and found a very small downturn in rotation in the core, while Couvidat et al. (2003) found a flat rotation profile down to 0.2 R_⊙ using combined GOLF, MDI and LOWL data.

Fletcher et al. (2003) investigated the problem of fitting the poorly-resolved higher-frequency low-degree mode splittings to integrated-sunlight observations such as those from BiSON. Using genetic fitting algorithms, they were able to reduce, though not eliminate, the bias towards higher splittings for these fits. They also found, in common with previous work, a strong anticorrelation between the estimated splitting value and its formal error, which would tend to cause overestimated splittings to be more heavily weighted in inversions.

García et al. (2004) considered two years of “sun-as-a-star” observations from early in the solar cycle, obtained from GOLF, GONG, MDI, VIRGO and BiSON, and were able to extract not only sectoral splittings but also a₃ and a₅ coefficients from the data, suggesting that it may be possible to infer differential rotation even in stars from which we will never have resolved data.

Chaplin et al. (2004) used artificial data to address the question the detectability of a rotation-rate gradient in the core. They concluded that, based on the best available data from ten years of observations, the difference between the rotation rate at 0.1 R_⊙ and 0.35 R_⊙ would be detectable only if it exceeded 110 nHz.

Chaplin et al. (2006) carried out an exhaustive “hare-and-hounds” exercise, in which one participant (the “hare” supplies the same set of artificial data to the others, the “hounds,” who then apply their various fitting methods without knowing the “true” answer, and compare the results. They obtained good agreement between the different techniques for l = 1, but systematic differences for the l > 1 splittings, which are attributed to different assumptions about the relative heights and spacing of the non-sectoral (∣m∣ < l) components.

5.8 Summary of the acoustic-mode results

To summarize, the best evidence we have so far seems to imply that the rotation rate between about 0.2 R_⊙ and the base of the convection zone is most likely approximately constant with radius and spherically symmetric. It is not possible to rule out a different rotation rate for the inner core, but there is no evidence from p-mode observations to support such a difference. Between about 0.2 R_⊙ and the base of the tachocline, no significant departure from rigid-body rotation has been found. As discussed by Eff-Darwich et al. (2002), for example the available constraints already seem to rule out the simplest models of hydrodynamic spin-down, which would show a detectable increase in the rotation rate below 0.3 R_⊙. Understanding both of the relationship between p mode splittings and the interior rotation, and of the care needed to measure them, has greatly advanced since the early days of helioseismology, but the rotation rate of the innermost nuclear-burning core remains uncertain.

5.9 Gravity modes

One possible way to improve the constraints on the core rotation would be to use g modes, or gravity waves, instead of p modes. Because these modes have their greatest amplitude in the solar interior, they should be much more sensitive to the core properties. Unfortunately, they also have very small amplitudes at the surface. The history of helioseismology is littered with unconfirmed reports of g-mode identification; see, for example, Delache and Scherrer (1983), van der Raay (1988),Thomson et al. (1995), and the review by Hill et al. (1991b). The most promising recent work has been carried out using long time series from the GOLF instrument aboard SOHO. Appourchaux et al. (2000b) placed an upper limit of 10 mm/s on g-mode amplitudes based on two years of observations, and Gabriel et al. (2002) reduced this limit further, to 6 mm/s, using 5 years of data. Most recently, García et al. (2007) report finding a pattern of peaks with constant spacing in period corresponding to the model-predicted spacing for l = 2 g modes with δl = 0, δn = 1, and with a splitting that they interpret as corresponding to a core rotation rate of 3–5 times the surface rate; however, this is still a preliminary result in need of confirmation.

In a related paper, Mathur et al. (2007) point out that the current predictions for low-order g-mode frequencies are much more consistent than was the case a decade earlier, resulting in a period for the fundamental g-mode between 34–35 minutes. This finding does make one wonder about the usefulness of the g-mode observations for discriminating among models; on the other hand, it lends somewhat more credence to the current identification.

6 The Tachocline

6.1 Observations

While the bulk of radiative interior appears to rotate almost as a solid body, the base of the convection zone at 0.71 R_⊙ coincides with a region of strong radial shear, above which the convection zone exhibits a differential rotation pattern that depends strongly on latitude and only weakly on depth. This shear layer is known as the tachocline, a term introduced to the literature by Spiegel and Zahn (1992), who attribute to D.O. Gough the correction of the earlier term “tachycline” (Spiegel, 1972). As is evident from the date of the latter reference, the notion of a shear layer at the bottom of the convection zone had been present in models for some time prior to its observational discovery, though its exact location was somewhat uncertain.

The existence of a layer of radial shear around the base of the convection zone, with approximately solid-body rotation below it, was first demonstrated by Brown et al. (1989), using the data of Brown and Morrow (1987); however, the significance of their results was quite low and they were at pains to point out that other interpretations of the data were possible. Dziembowski et al. (1989) used BBSO data to improve the picture of rotation at the base of the convection zone, again finding that the low-latitude rotation rate increased, and the high-latitude rate decreased, towards a common value at the base of the convection zone. The position of the base of the convection zone was determined by Christensen-Dalsgaard et al. (1991) using sound-speed inversions of helioseismic frequencies from the work of Duvall Jr et al. (1988) and Libbrecht and Kaufman (1988); their value of 0.713 R_⊙, confirmed by Basu and Antia (1997), has been accepted ever since.

The discovery of this shear layer (as pointed out by Brown et al.) offered a solution to the puzzle of the apparent absence of a radial gradient of rotation in the convection zone that could drive a solar dynamo, leading to speculation that the dynamo must operate in the tachocline region instead of in the bulk of the convection zone.

The tachocline lies in the region where modes of l ≈ 20 have their lower turning points, and the resolution of the inversions is quite low — about 5–10% of the solar radius in the radial direction. The thickness of the shear layer is therefore likely not to be resolved in inversions, and some ingenuity (and forward modeling) is required to estimate it and account for the effect of the finite-width averaging kernels in smoothing out the inversion inferences. The results of various efforts to parameterize the tachocline shape at the equator are summarized in Table 2. They mostly concur in placing the centroid of the shear layer slightly below the seismically-determined base of the convection zone, and its thickness at around 0.05 R_⊙. The largest value for the thickness, that of Wilson et al. (1996b), was obtained using forward calculation and direct combination of splitting coefficients rather than a true inversion, while the very low value of Corbard et al. (1999) was obtained using an inversion technique specifically designed to allow a discontinuous step in the rotation rate at the tachocline. The analysis of Elliott and Gough (1999) was somewhat different from the others, in that it involved calibrating a particular model of the tachocline against the inferred sound-speed rather than against a rotation profile.

Table 2: Tachocline radius r and width Γ.

Antia et al. (1998) and Corbard et al. (1999) found no significant evidence for a variation in the position or thickness of the tachocline with latitude, but Charbonneau et al. (1999) found a significant prolateness, with the tachocline (0.024 ± 0.004) R_⊙ shallower at latitude 60° than at the equator. Basu and Antia (2003) also found a slightly thicker and shallower tachocline at high latitudes, and speculated that the tachocline location might be discontinuous at the latitude (around 30°) where the shear vanishes and changes sign.

6.2 Models and the tachocline

Even the most generous estimates for the observed tachocline thickness are small enough to pose an interesting theoretical question: what prevents the shear from spreading further into the radiative interior, destroying the observed uniform rotation? The literature on tachocline modeling is extensive, far beyond the scope of this review. In brief, three main candidate mechanisms have been proposed: turbulent flows (Spiegel and Zahn, 1992); “fossil” magnetic fields (e.g., Gough and McIntyre 1998); and gravity waves, known to observational helioseismologists as g modes (e.g., Zahn et al. 1997), but all these scenarios have problems when considered as the sole mechanism. Recent advances in computing have made possible detailed three-dimensional simulations to explore these issues, but these models have not yet been able to reproduce a self-sustaining tachocline. For a review from a modeler’s perspective, see Miesch (2005). Also, a variety of discussions of tachocline models are collected in the book edited by Hughes et al. (2007).

7 Rotation in the Bulk of the Convection Zone

7.1 Observational history

The surface differential rotation, with the equator rotating faster than the poles, was known from, for example, sunspot tracking, long before helioseismology opened up the solar interior. Most models in the pre-helioseismology era predicted or assumed a rotation rate constant on cylinders parallel to the axis of rotation. This is a consequence of the so-called Taylor-Proudman constraint, a well-known result in fluid dynamics.

Duvall Jr and Harvey (1984) and Duvall Jr et al. (1984) observed from the South Pole, using only sectoral modes; their instrument used intensity images in a Calcium absorption line, scanning rather than imaging the whole Sun at once. Their main conclusion was that: “Most of the Sun’s volume rotates at a rate close to that of the surface”.

Brown (1985) had a different instrument, the Fourier Tachometer, which produced 100 × 100 pixel velocity images. Brown’s initial crude analysis of five days of data used cross-correlation, and expanded the multiplet frequencies using low-order polynomial fits; the results showed little sign of any depth variation in the differential rotation.

Duvall Jr et al. (1986), again using data from South Pole observations but now covering the full range of azimuthal orders, found values of the a₃ coefficient (the first-order measure of differential rotation) consistent with the surface rotation and rather larger than was consistent with the results of Brown (1985).

Brown and Morrow (1987), with 15 days (not all consecutive) of Fourier Tachometer data, could not distinguish between rotation on cylinders and latitude-dependence, but found that there was definitely less differential rotation in the radiative interior below the convection zone; their a₃ values were now closer to those of Duvall Jr et al. (1986), and they declared the previous ones erroneous. Brown et al. (1989) carried out a much more detailed analysis of the Brown and Morrow (1987) data, strengthening the evidence for mostly depth-independent rotation in the convection zone, as shown in Figure 16.

Figure 16:

Summary of rotation inferences from Brown et al. (1989) (reproduced by permission of the AAS).

Both the South Pole observations and those of Brown and collaborators were relatively noisy and of poor resolution; although they strongly hinted at a picture with little radial differential rotation in the convection zone and little differential rotation at all below it, other interpretations were possible.

Libbrecht (1989) published splittings from 100 days of BBSO observations in summer 1986, broadly confirming the results of Brown et al. (1989) with substantially smaller uncertainties. Dziembowski et al. (1989) inverted these data, and inferred a sharp radial gradient at the base of the convection zone and roughly constant rotation at each latitude above that. They also found a bump in the rotation rate in the middle of the convection zone, to which we will return below. Other inversions of the same data set were presented by Christensen-Dalsgaard and Schou (1988) and Libbrecht (1988b), with similar results, though not all the early inversions (cf. Korzennik et al. 1988; Sekii 1991) produce such recognizable results; this may be an example of the difficulty of using OLA-type techniques for data with insufficient higher-degree modes. Another (2dOLA) inversion of these data, shown in Figure 17, was carried out by Schou et al. (1992), who illustrated their averaging kernels; these were rather broad, but adequate to rule out a rotation-on-cylinders model. This paper was also the first to make the important point that the so-called “polar” rotation rate inferred from inversions is actually localized somewhat away from the pole.

Figure 17:

Rotation profile based on analysis of BBSO splittings, (Schou et al., 1992) (reproduced by permission of the AAS).

Gough et al. (1993) continued to challenge the observers to completely exclude rotation on cylinders, pointing out that it was possible to construct a cylindrical model that satisfied the constraint of the BBSO data, but Schou and Brown (1994b) showed that such a model could not be made consistent with both the Fourier Tachometer data and the gravitational stability of the rotating Sun.

Bachmann et al. (1993) analyzed Fourier Tachometer observations from 1989 and pointed out a “wiggle” in the splitting coefficients at ν/L ≃ 40µHz, (corresponding to a turning-point radius of about 0.85 R_⊙); attributed to daily modulation of the observations, this now well-known effect accounts for the “feature” seen in the middle of the convection zone in many inversions of single-site data.

Better data, with long time series free from daily modulation, were obviously needed before much more progress could be made, and with the advent of the GONG network in 1995 and the MDI instrument aboard SOHO in 1996 such data became available. Preliminary rotation profiles were presented by Thompson et al. (1996) for GONG and by Kosovichev et al. (1997) for MDI, both showing the now familiar pattern of almost-constant rotation in the convection zone, with shear layers both at the base of the convection zone and below the surface.

Schou et al. (1998) carried out a comprehensive analysis of the rotation profile based on the first 144 days of observations from MDI, using and comparing several different rotation inversion techniques with an input data set consisting of coefficients up to a₃₆ for p modes up to l = 194 and f modes up to l = 250. They were able to obtain consistent and robust results from the surface to about 0.5 R_⊙ at low latitudes; at higher latitudes the domain of reliability was shallower. Roughly speaking, the inversions could not be well localized within about 0.2 R_⊙ of the rotation axis. The results (Figure 18) showed that the rotation in the bulk of the convection zone, below 0.95 R_⊙, had a slow increase with radius at most latitudes, but was definitely incompatible with rotation on cylinders.

Figure 18:

Four inferred rotation profiles from the first 144 days of MDI observations (Schou et al., 1998): (a) 2dRLS, (b) 2dSOLA, (c) 1d×1dSOLA, (d) 1.5dRLS, from Schou et al. (1998) (reproduced by permission of the AAS).

7.2 The “polar jet”

In addition to the other general features described here, Schou et al. found some evidence for a “jet” of faster rotation at about 75° latitude and 0.95 R_⊙; although this was more obvious in some inversions than in others, it did seem to have a signature in the coefficients themselves (see also Howe et al. 1998). However, this feature was not reproduced in inversions of GONG data (e.g., Howe et al., 2000b), or even in inversions of MDI data analyzed with the GONG pipeline (Schou et al., 2002), and it is now believed to be an artifact related to the MDI data analysis.

7.3 GONG/MDI comparison

Once both GONG and MDI had been running for a few years, it became evident that the two projects were producing inferences of the interior rotation profile that were different in some significant details, particularly at high latitudes within the convection zone. Schou et al. (2002) carried out a careful comparison, taking data from three epochs at different phases of the solar cycle from each project and deriving rotational splittings or splitting coefficients from each, both with the usual algorithms and with those regularly used for the other project’s data, before using both RLS and OLA inversions. The results clearly showed that most of the discrepancies arose from the analysis pipelines rather from the data themselves. The “CA” peak-fitting algorithm used for the MDI data was able to extract modes from the GONG data to somewhat higher degrees and lower frequencies than the “AZ” algorithm could manage with either GONG or MDI input data. However, for both MDI and GONG data, the “CA” algorithm introduced an anomaly in the splitting coefficients centered at around 3.3 mHz, which in turn caused the inversion inferences to show a higher rotation rate deep in the convection zone at higher latitudes. Excluding these data brought the GONG and MDI data (analyzed with the “AZ” and “CA” pipelines respectively) into much better agreement, at the cost of somewhat degraded resolution. Restricting both data sets to the common mode set below 3 mHz reduced the discrepancies even farther, but did not remove the “jet” in the MDI data. Since the “jet” feature was only seen in the MDI data analyzed with the CA pipeline, however, the authors concluded that this feature was probably spurious.

7.4 Slanted contours

Although much of the debate in the early 1990s centered on discriminating between rotation constant on cylinders and rotation constant along radial lines, neither picture gave a complete description of the data. Gilman and Howe (2003) and Howe et al. (2005) pointed out that the differential rotation in the bulk of the convection zone, at least at low- to mid-latitudes, could be quite well described by saying that the contours of constant rotation lay at about a 25° angle to the rotation axis, as illustrated in Figure 19.

Figure 19:

Mean rotation profile from GONG data; contours of constant rotation (left), showing lines at 25° to the rotation axis as dashed lines, after Howe et al. (2005), and cuts at constant latitude as a function of radius (right), after Howe et al. (2000b).

Figure 20 compares idealized rotation profiles for the cylindrical, radial, and slanted-contour configurations.

Figure 20:

Idealized rotation profiles for rotation constant on cylinders (left), radial lines (middle), and lines at 25° to the rotation axis (right). The top row shows contours of constant rotation, while the lower row shows rotation rate as a function of radius at constant latitude for latitudes at 15° intervals from the equator (top) to 75° (bottom). The rotation rate is matched to the GONG inferences at 0.99 R_⊙ and smoothed to simulate the broadening effect of inversion resolution on the tachocline; the near-surface shear was not included.

7.5 Polar rotation

Another interesting feature revealed by the early GONG and MDI observations (Schou et al., 1998; Birch and Kosovichev, 1998) was that, while the surface rotation rate was mostly well described by the usual three-term expansion in the cosine of the colatitude θ, Ω(θ) = A + B cos² θ + C cos⁴ θ, (e.g., Snodgrass 1984) the rotation rate close to the poles was significantly slower than that. The authors speculated that this might be a result of drag from the solar wind, and that the effect might therefore disappear or become less marked at epochs of higher activity. In fact, though the inferred high-latitude rate did speed up during solar maximum — as seen, for example, in Howe et al. (2005) and in Figure 26 — it remained at all times lower than the extrapolation of the three-term fit.

7.6 Models and rotation in the convection zone

The interior rotation is only one part of the complex system that drives the solar cycle, but it is perhaps still the easiest part to measure in the solar interior; the meridional circulation can be directly measured only in the shallower subsurface layers, and buried magnetic fields can at best only be inferred indirectly. The differential rotation in the convection zone must arise from the interaction of convection cells and Coriolis forces, with the meridional motions playing an important part.

Early depictions of the solar dynamo (see, for example, Köhler 1974; Durney 1975) required a rotation rate increasing inward, and a meridional flow rising at the poles and sinking at the equator, in order to drive the solar cycle migration of the activity belts in the observed sense. This picture, taken together with rotation on cylinders, would have meant that the observed surface differential rotation was a superficial phenomenon, with the dynamo operating in the unobservable deeper layers. At this stage, there does not seem to have been a clear distinction made between the direction of the meridional circulation at the surface and the direction of migration of the magnetic activity belts during the solar cycle, which are of course now understood to operate in opposite directions; the poleward meridional flow at the surface was first measured by Duvall Jr (1979).

The models of Glatzmaier (1985) and Gilman and Miller (1986), which were among the first numerical simulations of solar rotation and the dynamo, have been cited, for example by Wilson (1992) as dating from “Prior to the advent of helioseismology,” but this is not quite correct. In fact, both these papers refer to the Duvall and Harvey data, and Gilman and Miller (1986) also mentions the observations of Brown (1985), suggesting that the model results could be consistent with the helioseismic observations if there were a layer of inward-increasing velocity below the surface and above the domain of the simulation. The simulations in both cases, like their precursors over the previous several years such as that described by Gilman and Miller (1981), produced rotation approximately constant on cylinders and increasing outward, which would result in a dynamo wave propagating poleward if the dynamo were operating in the bulk of the convection zone. The main message that modelers in the late 1980s seem to have taken from the observations was that the rotation rate was increasing outward, in agreement with the simulations of Gilman and Miller (1986) but in disagreement with the α-effect dynamo picture, which required a rotation rate increasing inward; see Parker (1987) for a review representing a theorist’s perspective on the state of play at this stage. This led Gilman and Miller (1986) to suggest (not for the first time; see also, for example, Galloway and Weiss 1981) that the dynamo might be operating in a thin layer at the bottom of the convection zone; this speculation was further reinforced by the later helioseismic inferences that clearly showed this shear layer, or tachocline (see Section 6) and the approximately radial configuration of the rotation in the convection zone.

Even quite recent global simulations of convection (Brun et al., 2004, for example,) still show some tendency towards rotation on cylinders, but the higher-resolution calculation of Miesch et al. (2008) mostly eliminates the cylindrical effect and produces a rotation pattern, based on giant convection cells, that after suitable temporal averaging looks quite solar-like, as illustrated in Figure 21.

Figure 21:

Three temporally averaged rotation profiles from the spherical-shell simulations of (a) Brun et al. (2004), (b) Browning et al. (2006), and (c) Miesch et al. (2008) (reproduced by permission of the AAS).

8 The Near-Surface Shear

One persistent puzzle in the measurements of rotation at the photosphere had been that direct Doppler measurements consistently gave somewhat slower rotation rates than the measurements made by tracing surface features. For example, Brown et al. (1989) summarized the results of Snodgrass (1983, 1984) as

$${{{\Omega _m}} \over {2\pi}} = 462 - 74{\mu ^2} - 53{\mu ^4}{\rm{nHz}}$$

(25)

for magnetic features and

$${{{\Omega _p}} \over {2\pi}} = 452 - 49{\mu ^2} - 84{\mu ^4}{\rm{nHz}}$$

(26)

for the surface plasma, respectively, where μ is the sine of the latitude. For an overview of such measurements, see Beck (2000). The usual explanation for the discrepancy is that while the Doppler techniques measure the velocity at the surface, the tracers such as sunspots are anchored in a faster-rotating layer deeper down. For example, Gilman and Foukal () noted that the observations implied a subsurface shear layer and suggested that this might arise from angular momentum conservation in the supergranular layer.

An extremely early attempt to measure the subsurface rotation was made by Rhodes Jr et al. (1979), when the identification of the 5-minute oscillations with p modes was still a relatively recent discovery. These authors used high-degree modes, probing about the upper 20 Mm (0.03 R_⊙) of the convection zone, and detected an inwards-increasing gradient. If these measurements are reliable, they represent the first detection of the subsurface shear. However, most of the early helioseismic measurements of the internal rotation profile were restricted to a degree range that did not allow the near-surface shear to be resolved in inversions. Rhodes Jr et al. (1990), attempting to measure the rotation in the bulk of the convection zone, also saw hints of a gradient, opposite to that seen at the base of the convection zone, below the surface, and Wilson (1992) used forward calculation techniques on the data of Brown and Morrow (1987) and Libbrecht (1989) to deduce that the rotation rate must increase inward immediately below the surface. We should remember, however, that at this time the picture of the internal rotation profile was not as clear as it is today, and it is not always obvious whether interpretation of the observations as gradients of rotation refers to the near-surface shear, the shear at the base of the convection zone, or some unresolved amalgamation of the two. Wilson, for example, was not arguing for a near-surface shear layer but against the model with rotation constant along radii.

With the advent of GONG and MDI, measuring modes to higher degrees than had previously been possible, the near-surface shear could be seen in global inversions; it is visible in the early results presented by Thompson et al. (1996) for GONG and by Kosovichev et al. (1997) for MDI, in both cases apparently changing sign at higher latitudes.

Schou et al. (1998) found clear evidence of the near-surface shear in inversions of MDI data. All the inversion methods agreed well on the shear at low latitudes, but at high latitudes the picture was complicated by the proximity of the submerged “jet” feature and the methods agreed less well. The disagreement may have been partly due to systematic errors in the splitting coefficients. In the comparisons of MDI and GONG data and analysis carried out by Schou et al. (2002), the high-latitude reversal of the shear is seen only in data analyzed with the “CA” pipeline; this may be partly because the “AZ” pipeline mostly fails to recover the splittings of the (narrow, low-amplitude) f-mode peaks, but the reversal persists in the MDI data even for the restricted common mode set.

The near-surface shear (down to about 15 Mm) was studied in detail by Corbard and Thompson (2002), using f modes from MDI data. They measured the slope of the rotation rate, close to the surface at low latitudes, as about −400 nHz/R_⊙, decreasing to a very small value by about 30° latitude and possibly reversing in sign at higher latitudes (though this result, seen in only the outer 5 Mm, was dependent on only the highest-degree modes, those with l ≥ 250). The low-latitude rotation rate was found to vary almost linearly with depth in the subsurface region, while if angular momentum was conserved in parcels of fluid moving with respect to the rotation axis, it would be expected to vary with the inverse square of the distance from the axis.

The near-surface shear is also accessible to the methods of local helioseismology, at least for latitudes below 50–60°. Basu et al. (1999) and Howe et al. (2006a) compared results from local ring-diagram analysis and global inversions and found, at latitudes ≤ 30°, quite good agreement between the dΩ/dr values obtained from local and global inversion results. However, although the slope from local measurements does show some variation with latitude (Figure 22), it by no means vanishes at 52.5°, the highest latitude at which the measurement is made. The ring-diagram results allow us to consider the northern and southern hemispheres separately, but Basu et al. (1999) found very little difference in the shear between the two hemispheres.

Figure 22:

Radial variation of the mean rotation rate after subtraction of the tracking rate, for global inversions (blue) and north-south averaged local inversions of MDI (green) and GONG (red) data at latitudes 0° (a), 15° (b), 30° (c), and 45° (d); similar to Howe et al. (2006a).

Some attempts have been made to use the near-surface shear to drive or at least contribute to a solar dynamo, for example by Brandenburg (2005), but Dikpati et al. (2002) showed that any dynamo contribution from the shear of the outer layers could only provide a fraction of the effect needed to power the solar cycle.

9 The Torsional Oscillation

9.1 The torsional oscillation before helioseismology

The so-called Torsional Oscillation is a pattern of migrating bands of faster- and slower-than-average zonal (i.e., parallel to the equator) flow associated with the equatorward drift of the activity belts during the solar cycle. It was first described by Howard and LaBonte (1980), who used 12 years (1966–1978) of full-disk velocity observations from the 150-foot tower at the Mount Wilson observatory and found evidence of a pattern of flow bands migrating towards the equator; the greatest concentration of active regions is associated with the poleward edge of the main equatorward-moving band. They initially interpreted the high-latitude variations as consisting of bands of faster rotation starting at the poles and taking a full 22-year Hale cycle to drift to the equator. Scherrer and Wilcox (1980a) and Scherrer et al. (1980), observing at the Stanford Solar Observatory, found no evidence of changes in the equatorial rotation rate for data from 1976–1979, but as this period was close to a solar minimum, and the resolution of the Stanford instrument was not high, this is neither surprising nor inconsistent with the results of Howard and LaBonte. LaBonte and Howard (1982) note that Scherrer and Wilcox (1980b) (at a AAS meeting), had “confirmed the existence of the global velocity field,” though this is not apparent from the latter’s published abstract.

A somewhat different pattern of velocity variations is seen when magnetic features rather than Doppler measurements are used to determine the surface rotation rate, as described for example by Komm et al. (1993a), who found that the pattern derived from magnetograms lay equatorward of that from Doppler measurements, with the slower-than-average bands coinciding with the zones of greater magnetic flux.

Mount Wilson Doppler observations since 1986, clearly showing the pattern of migrating zonal-flow bands, were presented by Ulrich (1998, 2001); see also Howe et al. (2006a) for updated results. The bands extend over about 10° in latitude, and have zonal velocities a few meters per second faster or slower than the surrounding material, corresponding to excess angular velocity of less than 0.5% of the overall rotation, or a few nanohertz.

9.2 Early helioseismic measurements

The first hints of the signature of the migrating flow bands in helioseismic data can be seen in the BBSO data (Woodard and Libbrecht, 1993), as was pointed out by Howe et al. (2000c), but these measurements do not give much information on the radial extent of the flows. Kosovichev and Schou (1997) found evidence of the flows, a few meters per second faster than the general rotation profile, in f-mode measurements from early MDI data; Giles et al. (1998) found a similar pattern using the time-distance technique of local helioseismology, while Schou (1998) and Schou (1999) clearly showed that these flows were migrating in a manner consistent with the Mount Wilson Doppler observations. The first radially-resolved evidence of zonal flow migration was reported by Howe et al. (2000c) for GONG and by Toomre et al. (2000) for MDI, while Howe et al. (2000a) combined MDI and GONG data and concluded that the equatorward-migrating part of the flow pattern (at latitudes below about 40°) penetrated to at least 0.92 R_⊙ (56 Mm below the surface). Antia and Basu (2000) also reported similar findings. Antia and Basu (2001) studied the evolution of the variations poleward of 50°, which had much higher amplitudes than the equatorward-moving flows and which showed signs of propagating poleward over time. The larger amplitude of the high-latitude signal may be related to the smaller angular momentum closer to the rotation axis.

9.3 Recent results

As more data accumulated, the signature of the torsional oscillation pattern in the helioseismic observations became clearer. Vorontsov et al. (2002) studied the evolution of the flows in MDI data from 1996 through 2001. They concluded that at least the high-latitude region of changing rotation involves the whole depth of the convection zone. The results on the radial extent of the flows at lower latitudes were less clear, with evidence that the bands of slower rotation might penetrate close to the base of the convection zone, while the bands of faster rotation appeared to reach about 0.9 R_⊙ but no deeper. Another interesting feature of that paper was the introduction of the use of 11-year sinusoids to characterize the variation of the rotation rate at any given location. This innovation had the useful effect of clarifying the pattern, making obvious the poleward propagation of the high-latitude flows even with data from little more than half a cycle. The existence of a weak third-harmonic component to the 11-year cycle, however, was not confirmed in later work.

Basu and Antia (2003) found similar results in MDI and GONG data up to 2002, as seen in Figure 23. These results also hint at another subtlety; at low latitudes, the phase of the flow pattern is not constant along radial lines. In fact, the variation in the lower part of the convection zone appears to lead that close to the surface by a year or two, with the low-latitude band of faster rotation following roughly the same 25° slant as the rotation contours. This tendency was further studied by Howe et al. (2005, 2006b), who compared inversions of MDI and GONG data with forward-modeled profiles based on different flow configurations, including some derived from dynamo models. Although some detail was lost and distorted due to the resolution and uncertainties in the inversions, the authors were able to conclude that the low-latitude branch probably penetrates through much of the convection zone, but is sufficiently displaced in phase at greater depths that the correlation between the surface pattern and that deeper down almost vanishes. In this work, the 11-year sinusoid analysis showed evidence of a second-harmonic component rather than the third harmonic reported by Vorontsov et al. (2002).

Figure 23:

Contour diagrams of constant rotation velocity residuals at 0.98R_⊙, obtained using two dimensional RLS inversion of the GONG data, from Basu and Antia (2003) (reproduced by permission of the AAS).

Figure 24:

Zonal flow pattern derived from MDI f-mode measurements, with smooth profile subtracted. Based on a figure from Schou (1999), updated and used by kind permission of J. Schou (2008, private communication.)

Figures 25, 26, and 27 show the variations in rotation rate, based on the results and figures in Howe et al. (2005, 2006b), but brought up to date with the most recent GONG and MDI observations available at the time of writing. The plots were prepared using the same 2dRLS inversion codes for both MDI and GONG medium-degree data, and 2dSOLA for MDI, that were used for the work of Howe et al. (2000a) and the other related papers. Figure 28 shows the phase and amplitude profiles for 11-year sine functions fitted to the rotation variations.

Figure 25:

Rotation rate after subtraction of a temporal mean at each location, as a function of latitude and time at selected depths, for OLA (top) and RLS (middle) inversions of MDI data, and for RLS inversions of GONG data (bottom).

Figure 26:

Rotation rate after subtraction of a temporal mean at each location, as a function of depth and time at selected latitudes. Latitudes are 0, 15, 30, 45, 60° from left to right; inversions are MDI OLA (top), MDI RLS (middle), and GONG RLS (bottom).

Figure 27:

Rotation rates at selected latitudes and depths as a function of time, after subtraction of a temporal mean. The results are from GONG RLS (black), MDI RLS (red), and MDI OLA (blue) inversions.

Figure 28:

Phase (left) and amplitude (right) of 11-year sine functions fitted to temporal variation of the rotation rate for OLA (top) and RLS (middle) inversions of around 11 years of MDI observations and for RLS inversions of GONG data (bottom).

9.4 Local helioseismology and the torsional oscillation

The torsional oscillation pattern, at least at lower latitudes and closer to the surface, is also suitable for measurements using the techniques of local helioseismology, in which short-wavelength, short-lived waves are used to infer the structure and dynamics of localized areas of the Sun. Because these waves do not penetrate very far below the surface, such techniques are restricted to the outer few megameters of the solar envelope, but this region can be studied in much greater detail and with shorter averaging times than is possible with global helioseismology.

Basu and Antia (2000) detected the zonal flow migration using MDI data and the ring-diagram technique (Hill, 1988), in which the displacement of three-dimensional acoustic power spectra derived from small areas of the solar disk is used to infer horizontal flows in both the zonal and meridional directions. Later, Haber et al. (2002) measured both the zonal flows and a corresponding modulation of the meridional flow pattern, as seen in Figure 29 (left). Beck et al. (2002), using the time-distance technique, which considers the correlations between oscillations at spatially separated locations, also found bands of meridional flow away from the activity belts associated with the zonal flow bands. Chou and Dai (2001) and Chou and Ladenkov (2005), using data from the Taiwan Oscillations Network [TON]), also found diverging meridional flows associated with the activity belts. Zhao and Kosovichev (2004) measured the zonal (Figure 29 right) and meridional flows with the time-distance technique, and reported meridional flow converging on the activity belts above a depth of 12 Mm, with diverging flows below 18 Mm, forming circulation cells around the activity belts. The presence of inflows into the activity belts was also observed at the surface by Komm et al. (1993b) and Komm (1994). Komm et al. (2005) studied the flows in about a year of high-resolution GONG (“GONG+”) data, and concluded that the overall flow pattern existed whether or not active regions were included in the analysis; in other words, the zonal flow bands and their associated converging/diverging meridional flows appear to exist independently of the flows in the immediate vicinity of strong active regions.

Figure 29:

Local helioseismic inferences of zonal flows close to the surface, from Haber (2002) (left) and Zhao and Kosovichev (2004) (right) (reproduced by permission of the AAS).

Figure 30:

Zonal flows since 1986, from Mount Wilson Doppler measurements (top), global helioseismic measurements from BBSO and MDI (middle) and MDI ring-diagram analysis (bottom). The color scale is in nHz.

Howe et al. (2006a) compared the results from ring-diagram analysis of the MDI data, global analysis of MDI and GONG data, and the Mount Wilson Doppler observations. They found very similar results for the north-south symmetrized flow pattern close to the surface in all three observations. Both the global and local helioseismic data indicated that the strength of the flow pattern did not fall off steeply below the surface.

It should be noted that the local helioseismic observations are somewhat prone to systematic errors, some of which follow the changing B₀ angle, or tilt of the solar rotation axis relative to the observer, as shown for example by Zaatri et al. (2006). This can result, for example, in a pronounced and almost certainly non-solar north-south variation of the zonal flow measurements, which is generally corrected for by subtracting suitable averages.

Some further features of the torsional oscillation pattern as we know it from a full cycle of observations from GONG and MDI (and nearly two cycles of surface Doppler observations) are worth noting:

• 1. The exact appearance of the pattern is quite sensitive to the background term that is subtracted. For example, compare the f-mode results shown in Figure 24, which were plotted as the difference from a smooth 3-term expansion of the rotation rate, with the plots in Figure 25, which were plotted by subtracting the temporal mean at each location.

• 2. Although the pattern repeats — of course not precisely — with each (approximately) 11-year activity cycle, each equatorward-migrating flow band exists for about eighteen years, emerging at mid-latitudes soon after the maximum of one cycle and finally disappearing at the equator a couple of years after the minimum of the following cycle; thus, the band of faster rotation associated with the activity of cycle 22 was still visible at the beginning of GONG and MDI observations in early cycle 23, and the band that is expected to accompany cycle 24 became visible around 2002 (if we look at the mean-subtracted residuals), or 2005–2006 (if we use the smooth-function subtraction). On the other hand, each poleward-moving branch seems to last only about nine years, appearing a year or so after solar minimum and moving to the pole before the next minimum.

• 3. Although the equatorward-migrating bands of faster rotation are clearly associated with the migrating activity belts of the magnetic butterfly diagram, the relationship is not completely straightforward. The new equatorward-propagating branch is clearly visible some years before noticeable new cycle active regions begin to erupt, and the phase/latitude profiles of the magnetic index and the velocity are very different. Also, as was noted by LaBonte and Howard (1982) and by Howe et al. (2006a), the strength of the torsional oscillation signal has not shown much change over the last few solar cycles, while the level of magnetic activity varies much more from one cycle to another.

• 4. Although the equatorward branch of the zonal flow migration pattern shows some relationship to the pattern of enhanced activity in the Fe xiv corona going back to 1973 (Altrock, 1997), the “extended solar cycle” seen in these observations starts at a much higher latitude, apparently about 70°, before migrating to the equator over about eighteen years; thus even the equatorward edge of these coronal activity bands seems to be at higher latitude than the observed new branch in the zonal flows that starts at about the same time.

• 5. Finally, we note that because the angular velocity changes associated with the torsional oscillation signal are relatively small compared to the difference in angular velocity between the surface and the bottom of the near-surface shear layer, while the amplitude of the signal does not decrease rapidly with depth, the magnitude of the shear at a given location varies by only a fraction of its value during the solar cycle. However, the fractional change in the shear is much greater than the fractional change in the rotation rate.

9.5 Models of the torsional oscillation

While observers, for example Howard and LaBonte (1980) and Ulrich (2001), have speculated that the torsional oscillation pattern might itself be part of the driving mechanism for the solar cycle, perhaps generating activity by shearing magnetic loops, modelers have generally seen it rather as a side-effect of the magnetic fields.

Schüssler (1981) and Yoshimura (1981) modeled the torsional oscillation as a result of the Lorentz force due to dynamo waves; according to the latter paper, the phenomenon would be important only close to the surface, and would have only equatorward, not poleward, moving bands. LaBonte and Howard (1982) objected to the Yoshimura model on the grounds that it would predict a strong correlation between the strength of the surface magnetic field and that of the velocity signal, which did not seem to be the case in the observations.

Küker et al. (1996) used a different mechanism to generate the torsional oscillation signal in their model, considering it as the response of the Reynolds stress on the time-dependent dynamo magnetic field rather than a direct effect of the large-scale Lorentz force. This model gave a very weak poleward branch for the torsional oscillation signal.

Once the flows had been shown observationally to penetrate well below the surface, Durney (2000) suggested that, “the pattern of torsional oscillations appear to have the potential of critically discriminating between different dynamo models as, e.g., the Babcock-Leighton and interface models.”

Covas et al. (2000) used a model in which the observed rotation profile was imposed and the rotation variations arises from the action of the Lorentz force of the dynamo-generated magnetic field on the angular velocity. They were able to simulate approximately solar-like patterns of zonal flow bands and magnetic activity. In subsequent papers they focused on the the possibility of so-called “spatio-temporal fragmentation” allowing cycles of different periods in different regions, and in calculations with no density stratification in the convection zone they found this to be feasible (Covas et al., 2001a). The effect was not too sensitive to uncertainties in the rotation law (Covas et al., 2001b, 2002), and somewhat sensitive to the boundary conditions at the outer surface (Tavakol et al., 2002). Adding density stratification (Covas et al., 2004) did not substantially change the results, though the amplitude of the oscillations in the deeper layers of the convection zone did decrease as the density gradient increased. However, they did find that introducing quite a small amount of α-quenching (magnetic feedback on turbulent convection) would suppress the torsional oscillation effect.

Spruit (2003) modeled the torsional oscillation pattern as a “geostrophic flow” driven by temperature variations near the surface associated with magnetic activity, and therefore having its greatest amplitude at the surface and falling to 1/3 of its surface value at 0.92 R_⊙. This model also accounts for the observed inflows into the activity belts. There are some problems in reconciling this model with the observations; it is difficult to see how the observed depth-dependent phase pattern could arise from a surface-originated cause, and the existence of the flows even at epochs where there are no active regions is also hard to explain, though Spruit suggested that the flows might be produced by unobserved small-scale and short-lived magnetic regions.

Rempel (2007) used a mean-field flux-transport dynamo model, with a model-derived differential rotation profile and meridional flow, to investigate the effects of various driving mechanisms for the torsional oscillation. The author concluded that the poleward-propagating branch of the pattern could be explained by a periodic forcing at mid-latitudes without any underlying migration of buried polar field. On the other hand, in this type of model the observed equatorward-propagating branch could not be reproduced without adding a thermal forcing after the manner of the Spruit (2003) model. Howe et al. (2006b) compared such a model with the observations, and found it not to be completely consistent with the observed interior behavior of the flows at lower latitudes.

10 Tachocline Variations

10.1 The 1.3 year signal

Howe et al. (2000b) reported finding variations of the equatorial rotation rate close to the tachocline with a 1.3 year period during the early years (1995–1999) of GONG and MDI observations. The strongest signal was seen at 0.72 R_⊙, with a weaker anticorrelated signal below the tachocline at 0.63 R_⊙. At higher latitudes, there was also an apparent 1-year periodicity. The signal was more clearly seen in the GONG data, and due to the different temporal sample of the MDI data it was difficult to make a quantitative comparison, but the visual appearance of similar variations in both data sets was quite persuasive. Figure 31 extends the data up to the present for the equatorial locations just above and below the tachocline.

Figure 31:

Rotation-rate residuals at the equator at 0.72 R_⊙ (top) and 0.63 R_⊙ (bottom), for RLS (filled) and OLA (open) inversions of MDI (red triangles) and GONG (black circles) data.

Because of the role of the tachocline region in the dynamo, as well as the coincidence of the period with that seen in some heliospheric and geomagnetic observations (Silverman and Shapiro, 1983; Richardson et al., 1994; Paularena et al., 1995), this claim attracted considerable interest, inspiring modelers such as Covas et al. (2001a) to try to build models in which different periods could exist at the top and bottom of the convection zone. However, Antia and Basu (2000) and Basu and Antia (2001), with a slightly different analysis of the same MDI and GONG data, reported finding no significant variations (Basu and Antia (2001) did see a signal somewhat similar to that reported by Howe et al. (2000b) but did not consider it significant).

Moreover, the periodic signal disappears in the post-2001 data even in the original authors’ analysis (Toomre et al., 2003; Howe et al., 2007), as shown in Figure 32, and it seems likely that the high-latitude 1-year period was an artifact. Intermittency in short-period variations is a known phenomenon in the geomagnetic-index data (Silverman and Shapiro, 1983) and does not in itself imply that the phenomenon was not real. It will be interesting to see whether the oscillation will reappear in the new solar cycle.

Figure 32:

Sine-wave power in the rotation rate residuals from RLS inversions of GONG data, at 0.72 R_⊙, 0°, plotted as a function of frequency for a) 1995–2000, b) 1995–2003, c) 1995–2005, d) 2000–2005.

10.2 Tachocline jets

Christensen-Dalsgaard et al. (2004) searched for evidence of jets close to the tachocline, which are predicted, for example, by the model of Dikpati et al. (2004). Using GONG data they reported finding possible evidence of a jet at the tachocline, migrating equator-wards by about 30 degrees in two years but not at the same latitude as the surface activity belts. The significance and meaning of this finding remain unestablished.

10.3 Angular momentum variations

Given estimates of both density and rotation as functions of depth and latitude, one can calculate the solar angular momentum locally or globally. Of course, such calculations will reflect, and in some cases enhance, any errors in the input data, and should therefore be approached with caution.

Komm et al. (2003) investigated the angular momentum variation based on the inversions of GONG and MDI data used by Howe et al. (2000b, a) and found variations reflecting the torsional oscillation well into the convection zone and 1.3 year variations close to the tachocline. Because the density increases steeply with decreasing radius, variations at greater depths will be more strongly seen in the angular momentum than in the rotation rate, but it should be remembered that no new information has been added to the data.

Lanza (2007) approached the problem from the other direction, considering the role of angular momentum transport in the modeling of the torsional oscillation.

Antia et al. (2008) investigated temporal variations of the solar kinetic energy, angular momentum and higher-order gravitational multipole moments as derived from helioseismic inferences of the internal rotation rate; they found variations on the time scale of the solar cycle (but not the 1.3 year cycle), with some discrepancies between MDI and GONG results. They also speculate that the kinetic-energy changes might contribute to the observed irradiance variations during the solar cycle; however, it is not clear that such a contribution is needed, as the usual view is that the solar-cycle variation in irradiance can be modeled simply from the effects of sunspots and plage on the surface, as discussed, for example, by Jones et al. (2008).

11 Summary and Discussion

Since the 1970s, helioseismology has provided several insights into the interior solar rotation: the approximately-rigid rotation of the radiative interior; the differential rotation throughout the convection zone; the thin tachocline; the extension of the surface torsional oscillation throughout the convection zone. More than once, these discoveries have overturned theoretical expectations, inspiring modelers to improve their calculations in an effort to reproduce the observed behavior. Because of the surprising nature of many of the findings, it has been important to have more than one source of observations, so that it is possible to distinguish between real solar features — especially the unexpected ones — and systematic error.

It may be that in the future solar cycle 23, with MDI and GONG operating in parallel, will be seen as a golden age of helioseismology. At the time of writing, we eagerly anticipate the launch of the Solar Dynamics Observatory [SDO] with its Helioseismic and Magnetic Imager [HMI], a successor to MDI that will provide near-continuous helioseismic observations at higher resolutions than ever before and may help in unraveling the relationships between active region flows, magnetic fields, and geoeffective solar activity as well as providing a continued watch on the longer-term variations in the solar velocity fields. Sadly, however, current plans call for both GONG and MDI to cease to collect data soon after the successful launch of SDO, which would leave HMI without any independent cross-checks, while on the low-degree front the BiSON network has recently lost its funding and there are no new dedicated low-degree space-based instruments currently scheduled.

There are still areas — such as the strength of the near-surface shear at high latitudes, the rotation of the inner core, and any inhomogeneities and changes in the tachocline — that remain unclear. Furthermore, a complete numerical model of the solar dynamo — vital for any long-term predictive capability — is still lacking, and helioseismic observations still have an important part to play in constraining such models as they develop.