Solar Interior Rotation and its Variation

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.


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  (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.

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 = 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.
As was first discovered by Deubner (1975), the velocity or intensity variations at the solar surface have a spectrum in − or − 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 " modes" can be classified by their radial order , spherical harmonic degree , and azimuthal order ; as discussed, for example, in Section 2.2 of Birch and Gizon (2005), the radial displacement of a fluid element at time , latitude and longitude can be written in the form where is the radial eigenfunction of the mode with frequency and ( , ) is a spherical harmonic. As seen in Figure 4, the power in the spectrum falls along distinct "ridges" in the − plane, each ridge corresponding to one radial order. The modes making up the = 0 ridge are the so-called modes, which are surface gravity waves. The modes, so called because their restoring force is pressure, are excited at the surface and have their largest amplitudes there. Another class Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 of modes, the 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). The longer the horizontal wavelength -and the lower the degree -the more deeply the modes penetrate, with the radial = 0 mode going all the way to the core of the Sun (but providing no rotational information), while modes with ≥ 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 used for local helioseismology. The lower turning point radius, , is a monotonic function of the phase speed / , where = √︀ ( + 1) ≈ + 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.
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 ( ≤ 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.

Differential rotation and rotational splitting
The Sun's rotation lifts the degeneracy between modes of the same and different , 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 , ≡ − , − + , is proportional to the rotation rate multiplied by . Figure 7 shows a typical − acoustic spectrum of GONG data at = 100. The effect of the rotation causes the ridges at each 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.
Because modes of different values sample different latitude ranges, with the sectoral (| | = ) modes confined to a belt around the equator and the zonal or = 0 modes reaching to the poles, as illustrated in Figure 8, we can measure the rotation as a function of latitude.
A given ( , ) multiplet consists of 2 + 1 frequency measurements if each ( , ) 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 ( ) as a polynomial expansion, for example, where the basis functions are polynomials related to the Clebsch-Gordan coefficients 0 by (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 1 coefficient describes the rotation rate averaged over all latitudes, and the 3 and higher coefficients describe the differential rotation.

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 , azimuthal order , and radial order is given by 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 where is an apodization function and ≡ √︀ cos 2 + sin 2 sin 2 represents the fractional distance from disk center in the solar image. The line-of-sight projection factor is = √︀ 1 − 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 ( , , ′ , ′ ) for leakage from the ′ , ′ mode to the , spectrum can be computed as Symmetry properties in this expression lead to some simple exclusion rules; for example, leaks with odd | + | (where ≡ − ′ and ≡ − ′ ) are not allowed. One example of the importance of the leakage is in the contribution of the so-called -leaks ( = 0, = ±2) to the estimation of low-degree splittings. As pointed out, for example, by , these leaks are strongest for small | |; they are also asymmetrical, especially for | | = , where the = peak has an = − 2 leak on one side and no counterbalancing = + 2 leak on the other. Especially for = 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 = 0 mask) can detect modes of 0 ≤ ≤ 5, though the sensitivity falls off rapidly for > 1. All these modes appear in a single acoustic spectrum; the instruments are sensitive to odd-modes for odd and to even-modes for even , with the sectoral, or | | = , 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 -leaks at frequencies above about 2 mHz; for higher-degree modes the leakage between modes of adjacent becomes a problem, as the ridges become both broader, and more closely spaced in frequency, at around = 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.

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 , multiplet,Ω , as follows: 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 / .

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.

The inversion problem
The basic 2-dimensional rotation inversion problem can be stated as follows: we have a number of observations , from which we wish to infer the rotation profile Ω( , ) where is distance from the center of the Sun, and is (conventionally) colatitude. Each datum is a spatially weighted average of the rotation rate: where ⊙ is the solar radius, the error term corresponds to the noise and measurement error in the data, and 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 = cos , 2 = ( + 1), is the radial displacement for the eigenfunction of the mode, −1 is the horizontal displacement, and ( ) is the density (see Figure 10 for illustrations of sample kernels).
The aim of the inversion is to find where ( 0 , 0 ) is the location at which the inferred rotation rateΩ is to be found and the 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.

Averaging kernels
By substituting Equation (11) into the RHS of Equation (14) we obtain Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 where is the averaging kernel for the location ( 0 , 0 ). 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 , 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., 1992Schou et al., , 1994.

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 , the formal uncertainty on the inferred profile is given by In the (usually unrealistic) case where the errors on the input data are all equal, we can write where the "error magnification" is given by 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 for that datum is non-zero.) In the one-dimensional case, the correlation between the errors for two points 0 and 1 is given by 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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 widelyused approaches to balancing these criteria are "regularized least squares" (RLS) and "optimally localized averaging" (OLA).

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 with 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 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 0 , 0 . Figure 11 illustrates typical averaging kernels for a 2dRLS inversion of an MDI data set.

Optimally localized averaging
In the Subtractive OLA (SOLA) approach Gilbert, 1968, 1970), the minimization is applied to the difference between the actual averaging kernels and a target kernel , for example a 2-dimensional Gaussian or Lorentzian function. In this case Thompson, 1992, 1994) the function minimized is 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.
Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 Another approach, older, and more computationally expensive, is the Multiplicative OLA (MOLA) described by 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.

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.

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 is to first order proportional to the rotation rate multiplied by | |; since the only mode whose latitudinal kernel reaches the pole is the = 0 mode, which has no longitudinal structure and so can convey no rotational information, and the modes of small | |/ 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 ( ≤ 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.

Observations: A Brief Historical Overview
Systematic helioseismic observations stretch back nearly 30 years, as illustrated in the schematic chart in Figure 13. 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, , 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 andNoyes, 1985, 1988) aboard the Solar Maximum Mission spacecraft, which took helioseismic measurements in 1980 and 1984 -1985, the Mark I instrument in Tenerife , 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 , operated from 1989 -2003.
The Global Oscillation Network Group (GONG)  has been collecting continuous, high-duty-cycle observations of the medium-degree modes since 1995, using a sixstation 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 highquality 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 mediumdegree 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.

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 2 ; confusingly, the nexthighest term, 4 , 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 Δ = eq − pol is given by where eq and pol are the equatorial and polar radii, respectively, and 0 is the radius of the best sphere passing through eq and pol . If there is an additional contribution from the surface rotation this expression becomes The units of and Δ 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 measurementsa 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 −5 , 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) Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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 −6 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 2 and 4 terms should be based on a simple differential rotation profile deduced from surface measurements, obtaining predicted values of between 1 and 1.5 × 10 −7 for 2 and between 2 and 5 × 10 −9 for 4 depending on the size of the convective envelope. Dicke et al. (1986Dicke et al. ( , 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 −7 for 2 and 9.8 × 10 −7 for 4 . 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 4 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.

Early low-degree helioseismic results
Around the early 1970s there were numerous attempts to search for global -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 160minute 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 longerperiod 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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 and different , without resolving separate = 0 and = 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 = 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 ( = 0, = 1) mode remains unobserved to this day, although some low-degree modes with single-digit 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 = 1 modes and five in = 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 = 1 and = 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 = 0 peaks appeared narrower than those of = 1 and = 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 = 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 fiveminute 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.

Resolved-Sun measurements
In the meantime, resolved-Sun observations provided some information about the rotation in the radiative interior. Duvall Jr and  reported observations at Kitt Peak, from 10 -26 May 1983, for degrees 0 ≤ ≤ 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 = 6, with an unexplained bump at = 11, followed by an increase at lower degrees up to a value of 660 nHz Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 for = 1. These data, inverted , 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 = 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). 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.

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.
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.

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 = 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: Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1   (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 coefficientsthe 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.
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) -1 on the measurement of the = 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 LOWL instrument, an imaging instrument optimized for lower degrees. They obtained splittings for 1 ≤ ≤ 100, and inferred a rotational profile down to 0.2 ⊙ , finding a rotation rate that barely varied with radius between 0.2 ⊙ and 0.6 ⊙ , apart from a low-significance bump around 0.4 ⊙ .
Eff-Darwich and  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.18 ⊙ . Charbonneau et al. (1998) used a genetic forward-modeling approach to analyze the LOWL data, with results favoring a rigidly-rotating core.

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 lowfrequency 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 ⊙ , 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 , 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 ⊙ 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 3 and 5 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 rotationrate 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 ⊙ and 0.35 ⊙ 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 = 1, but systematic differences for the > 1 splittings, which are attributed to different assumptions about the relative heights and spacing of the non-sectoral (| | < ) components.

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 ⊙ 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 -mode observations to support such a difference. Between about 0.2 ⊙ 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 ⊙ . Understanding both of the relationship between 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.

Gravity modes
One possible way to improve the constraints on the core rotation would be to use modes, or gravity waves, instead of 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 ofmode 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 -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 = 2 modes with = 0, = 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 -mode frequencies are much more consistent than was the case a decade earlier, resulting in a period for the fundamental -mode between 34 -35 minutes. This finding does make one wonder about the usefulness of the -mode observations for discriminating among models; on the other hand, it lends somewhat more credence to the current identification.

Observations
While the bulk of radiative interior appears to rotate almost as a solid body, the base of the convection zone at 0.71 ⊙ 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 ⊙ , 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 ≈ 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 finitewidth 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 ⊙ . 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. 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) ⊙ 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.

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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 Reference Kosovichev (1996) 0  (Spiegel and Zahn, 1992); "fossil" magnetic fields (e.g., Gough and McIntyre 1998); and gravity waves, known to observational helioseismologists as 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).
Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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  and  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 , again using data from South Pole observations but now covering the full range of azimuthal orders, found values of the 3 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

Both the South Pole observations and those of Brown and collaborators were relatively noisy
Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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.  (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.
Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 Bachmann et al. (1993) analyzed Fourier Tachometer observations from 1989 and pointed out a "wiggle" in the splitting coefficients at / ≃ 40 Hz, (corresponding to a turning-point radius of about 0.85 ⊙ ); 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  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 36 for modes up to = 194 and modes up to = 250. They were able to obtain consistent and robust results from the surface to about 0.5 ⊙ 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 ⊙ of the rotation axis. The results ( Figure 18) showed that the rotation in the bulk of the convection zone, below 0.95 ⊙ , had a slow increase with radius at most latitudes, but was definitely incompatible with rotation on cylinders.

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 ⊙ ; although this was more obvious in some inversions than in others, it did seem to have a signature in the coefficients themselves (see also ). 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 , and it is now believed to be an artifact related to the MDI data analysis.

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.

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 20 compares idealized rotation profiles for the cylindrical, radial, and slanted-contour configurations.

Polar rotation
Another interesting feature revealed by the early GONG and MDI observations 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 , Ω( ) = + cos 2 + cos 4 , (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.

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.

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 (1983Snodgrass ( , 1984 as for magnetic features and Ω 2 = 452 − 49 2 − 84 4 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 (1979) 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 modes was still a relatively recent discovery. These authors used high-degree modes, probing about the upper 20 Mm (0.03 ⊙ ) 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  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) -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 modes from MDI data. They measured the slope of the rotation rate, close to the surface at low latitudes, as about −400 nHz/ ⊙ , 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 5 Mm, was dependent on only the highest-degree modes, those with ≥ 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 Ω/ 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. 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.
Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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-thanaverage 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 zonalflow bands, were presented by Ulrich (1998Ulrich ( , 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.

Early helioseismic measurements
The first hints of the signature of the migrating flow bands in helioseismic data can be seen in the BBSO data , as was pointed out by Howe et al. (2000c), but these measurements do not give much information on the radial extent of the flows.  found evidence of the flows, a few meters per second faster than the general rotation profile, in -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 ⊙ (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 highlatitude signal may be related to the smaller angular momentum closer to the rotation axis.

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 ⊙ 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. (2005Howe et al. ( , 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 Figure 24: Zonal flow pattern derived from MDI -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.) 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). Figures 25,26,and 27 show the variations in rotation rate, based on the results and figures in Howe et al. (2005Howe et al. ( , 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.

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, shortlived 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1    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. 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 0 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 -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) 11year 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, Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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.

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 socalled "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, and somewhat sensitive to the boundary conditions at the outer surface . 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 ⊙ . 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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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.  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 ⊙ , with a weaker anticorrelated signal below the tachocline at 0.63 ⊙ . 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.

Living Reviews in Solar Physics
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 Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 finding no significant variations  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. Living Reviews in Solar Physics http://www.livingreviews.org/lrsp-2009-1 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.

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).

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 featuresespecially 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.

Acknowledgements
This review has made use of NASA's Astrophysics Data System. This work utilizes data obtained by the Global Oscillation Network Group (GONG) program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. The Solar Oscillations Investigation (SOI) involving MDI is supported by NASA grant NAG5-13261 to Stanford University. SOHO is a mission of international cooperation between ESA and NASA. NSO/Kitt Peak magnetic data used here are produced cooperatively by NSF/NOAO, NASA/GSFC and NOAA/SEL. The Mt. Wilson observations have been supported over several decades by a series of grants from NASA, NSF and ONR and are currently supported by NSF/ATM. The Mt. Wilson observatory is managed by the Mt. Wilson Institute under agreement with the Observatories of the Carnegie Institution of Washington.
BiSON has been funded by the U.K. Particle Physics and Astronomy Research Council.