Abstract

To investigate whether the Sun's global meridional circulation is primarily mechanically or thermally driven, we solve a hydrodynamical model to determine what meridional flow pattern is most likely to occur for the Sun's relatively well-observed differential rotation. We include Coriolis forces and turbulent Reynolds stresses, and hold the inclusion of thermodynamics for a future work. We find that the steady-state solution from this model is sensitive to the magnitude of the density decline with radius, yielding a single flow cell with poleward surface flow for density declines of less than about a factor of 4 within the convection zone, but two cells, with a reversed circulation in high latitudes, for all higher density declines, including that for the adiabatically stratified solar convection zone. This two-celled pattern is virtually independent of the magnitude of the decline of the turbulent viscosity with depth within the convection zone. For a solar-like density increase with depth, two cells are found for a wide variety of differential rotations, including that closest to solar observations. We compare our meridional circulation with axisymmetric convective rolls that occur in 3D global convective simulations, which tend to align with rotation axis. For solar-like turbulent viscosity, the meridional flow speeds are much larger than observed, implying that an additional physical mechanism is needed that works against the meridional flow. Candidates include negative buoyancy forces, anisotropic small-scale turbulent diffusion of momentum or heat and latitude-dependent radial heat flux by solar convection influenced by rotation. The role of these additional mechanisms in determining flow speed and in producing latitudinal and radial structures in the meridional circulation will be investigated in a future paper.

INTRODUCTION

Analysis of long-term meridional flow data available at Mount Wilson Observatory (MWO) shows a reversed meridional flow cell at high latitudes, poleward of 60° or so, during solar cycle 19, 20, 21 and 22 (Ulrich 2010). The second cell disappeared during most of cycle 23, but has reappeared during the ascending phase of the current cycle 24. Surface Doppler measurements are more able to measure meridional flow at the highest latitudes than are helioseismic measures, since latitudinal flow is nearly parallel to the line of sight there. Accounting for solar B0 angle variations can interfere with helioseismic determinations of polar meridional flow (Komm et al. 2013). Measurements of meridional flow by tracking magnetic features has now been extended to very high latitudes (Rightmire-Upton, Hathaway & Kosak 2012). These show that for one year from 2010 April 2010 to 2011 March, there is no evidence of a reversed cell up to 85° latitude (Rightmire-Upton et al. 2012), but earlier measures using similar techniques (Hathaway & Rightmire 2011) do show such reversals for extended periods in both hemispheres.

It is not our intent here to argue the relative merits of the various methods of measuring meridional circulation, or what measurement most accurately captures the flow of the plasma. It is clear that there is no consensus at the present time about what the meridional circulation is in polar latitudes at the surface of the Sun. From the point of view of solar dynamo theory, and understanding how polar magnetic fields on the Sun evolve, it is very important to determine what the circulation pattern is normally in the steady state and as a function of time.

The latitudinal extent as well as the speed of the primary meridional flow cell have enormous influence in determining the duration of a cycle in a class of dynamo, called the flux-transport dynamo. Meridional circulation plays the role of a conveyor belt in such models, carrying the dynamo-generated magnetic flux from the surface to the pole and then to the bottom of the convection zone. The unusually long duration of cycle 23 followed by a long minimum at its end has been explained as being due to the primary cell of plasma flowing all the way to the pole during most of cycle 23 (Dikpati et al. 2010).

Observational evidence of meridional circulation in the Earth's atmosphere has existed for at least a century. Such flows, called Hadley and Ferrel cells, were postulated for explaining the mean zonal winds at least two centuries before that (Lorenz 1967). Reliable observations of meridional flow in the Sun have been available only for the past two decades, so models that incorporate them are of relatively recent origin.

There have been many axisymmetric ‘mean-field’ type calculations of solar differential rotation and meridional flow beginning in the 1960's. This early work is reviewed in Rüdiger (1989) and Stix (1989). These early models are of basically two types: those that are purely mechanical, which rely on anisotropic turbulent viscosity, such as Köhler (1970), and those that rely upon latitude-dependent radial convective heat transport, such as Durney & Roxburgh (1971) and Belvedere, Paterno & Stix (1980). Both models are able to simulate the surface latitudinal rotation gradient of the Sun, but the meridional circulation produced by both is towards the equator near the outer boundary, rather than towards the poles as observed. In Köhler (1970), if the sense of anisotropy is reversed, the surface meridional flow is poleward, but then the poles rotate faster than the equator. In Durney & Roxburgh (1971) and Belvedere et al. (1980), the surface temperature difference between equator and pole is substantially larger than solar observations indicate.

A reason these early models need equatorward flow near the surface is that, in order to transport angular momentum towards the equator to maintain the high rotation there, the angular momentum per unit mass of the fluid rings moving towards the equator must be higher than that of the poleward moving fluid at the same latitude. This is achieved by equatorward flow being at a larger radial distance from the centre of the Sun, so the fluid there has a larger moment arm relative to the solar rotation axis. A later model (Pidatella et al. 1986) does get poleward surface flow with a reversed latitudinal temperature gradient, apparently by relying upon anisotropic eddy viscosity to maintain the equatorial acceleration. Both meridional flow speed and pole–equator temperature differences are within observational limits. However, this model usually gives two cells stacked in radius, each of which reaches from equator to pole. Very recent observations from Solar Dynamics Observatory/Helioseismic and Magnetic Imager data are indicating such cells stacked in depth (Zhao et al. 2013); therefore, the theoretical reasons for such a pattern need to be reinvestigated in the context of our improved knowledge of solar interior.

A rather different approach to mean-field theory for stellar differential rotation and meridional circulation was taken by Rüdiger et al. (1998), in which they calculated differential rotation arising from turbulent Reynolds stresses while ignoring the dynamical effect of Coriolis forces from meridional circulation. They found that for the turbulent viscosity needed to give the observed differential rotation, the meridional circulation that would result due to Coriolis forces is orders of magnitude larger than solar observations allow. Conversely, if the turbulent viscosity is chosen to limit the meridional flow to solar values, its value would exceed solar estimates by two orders of magnitude. They conclude that an additional force is needed that opposes the large meridional flow speed to bring it within observational limits when the model correctly determines the amplitude of differential rotation.

The most recent results are found in Kitchatinov & Olemskoy (2011). Their primary focus has been on global solutions for differential rotation driven by the ‘Λ effect’ or non-diffusive Reynolds stresses, and meridional circulation driven by Coriolis forces and latitudinal entropy gradients, all working against turbulent diffusion represented by the diffusive parts of the velocity correlation stress tensor. In Kitchatinov & Olemskoy (2011) it is shown that these combined effects lead to a single celled meridional flow with poleward flow all the way to the poles of about the correct magnitude, coupled with a differential rotation that conforms to helioseismic observations that include radial shear layers at the bottom and top of the convection zone and the characteristic tilt of interior angular velocity contours of ∼25° from the outward radial direction. Almost no evidence of a second reversed cell near the poles has been found in this or most previous model studies.

The mean-field models just described can be thought of as part of a hierarchy of models, the most complex of which are 3D global convection models, such as the ASH codes (Miesch 2005). Then come the mean-field models, intermediate in complexity, which solve for meridional circulation and differential rotation that are driven by parametrized turbulent Reynolds stresses and/or turbulent heat transport. In this paper, we build an axisymmetric mean-field model to examine the relative roles played by all forces acting in meridional planes to determine the meridional circulation.

Still another approach is to estimate the meridional circulation from profiles of the convergence or divergence of the turbulent Reynolds stresses found in the azimuthal equation that describes angular momentum balance. This is the so-called gyroscopic pumping mechanism. Conversely, it is possible, in principle, to estimate the turbulent Reynolds stresses that must have driven differential rotation, by invoking the gyroscopic pumping mechanism in reverse. In the end, the combination of meridional circulation and differential rotation must satisfy balance equations in both meridional plane and the azimuthal direction.

In contrast to all previous models, our model here takes advantage of the fact that from helioseismic measurements the differential rotation of the Sun as a function of latitude and radius within the solar convection zone is rather well known. We model the meridional flow in global spherical geometry for a deep convection zone, including mechanical forces to determine the form and amplitude of meridional circulation in the Sun from all the physics known from observations. With this model we address the question of what form of meridional circulation is most likely to occur. Will there be just one primary cell between equator and pole, or multiple cells stacked in latitude or depth? We hold for the next paper the study of the effect of buoyancy forces arising from thermodynamic effects.

MODEL AND EQUATIONS

Since the bulk of the solar convection zone is very close to adiabatic in stratification (∼10−5), we take it to be adiabatic and therefore remove thermodynamics from the problem. This is very different than for global dynamics of the Earth's lower atmosphere, in which the stratification is perhaps 20–30 per cent below adiabatic stratification, implying that thermodynamics must be taken into account. As a result of this choice, we can solve a purely mechanical system, consisting of the three equations of motion and the mass continuity equation. Within this mechanical system, we can experiment with different density and turbulent viscosity profiles with radius to determine what meridional circulation pattern is favoured for different profiles of these quantities and differential rotation.

In limiting our model to be purely mechanical, neglecting thermal effects, we are guided by the theoretical discussion and results in Rempel (2005a), where the most extensive analysis of the role of departures from the adiabatic gradient has been carried out. Rempel (2005a) showed that, in a model including thermodynamics in the tachocline or convection zone or both, the resulting differential rotation is quite sensitive to the details of the departures from adiabatic – where they occur, their magnitude and whether they are positive or negative (super or subadiabatic). There is very little guidance from either observations or convection theory to constrain these departures from an adiabatic stratification. But by specifying the differential rotation, we avoid this uncertainty, taking advantage of the relatively greater observational knowledge of the differential rotation throughout the convection zone and tachocline. Each differential rotation used would be itself driven by a specific profile of non-viscous turbulent Reynolds stress, working against the turbulent viscous part of this stress.

Furthermore, as pointed out by Rempel (2005a), in a mean-field model, it is highly desirable to avoid exciting axisymmetric convective rings that would be all that remain of full 3D global convection. Isolated from the full 3D explicitly calculated global convection zone, these rings would be a very unrealistic representation of the convection, and would only obscure the patterns of meridional circulation driven by non-convective effects. Rempel (2005a) avoided them by always making sure that the Rayleigh number for his mean-field model was subcritical so that axisymmetric convective rings were not excited. Here we achieve the same goal by initially leaving out the thermodynamics altogether. When thermodynamics are included, buoyancy forces are most likely to act to brake the meridional circulation, not drive it.

There are two relatively distinct effects of density variations: one is associated with buoyancy effects and the other is inertial effects, due to the basic and large radial density variation through the convection zone. Here, we study the inertial effect of density variation. It weights the Coriolis forces from differential rotation in substantially different ways than happens in the incompressible constant density case, and we can learn from that, separately from the effects of buoyancy.

Once the differential rotation is specified, we can solve for the meridional circulation that is consistent with all forces in the meridional plane for each differential rotation we choose. Related ideas were expressed in Durney (1996). In each case, as also argued by Rempel (2005a), the different profiles of Coriolis force from the differential rotation, particularly the radial component, are not quite in balance with other forces in the meridional plane [turbulent viscous, pressure gradient, and buoyancy (when there are departures from adiabatic stratification)], and so each will drive a different meridional circulation. Just as with departures from adiabatic stratification and latitudinal entropy gradients, the buoyancy forces are hard to constrain by observations or convective theory. Therefore, there is value in determining the meridional circulation that results only from forces that we can constrain reasonably well from observations and theoretical considerations. We return to the question of buoyancy forces in the discussion of results.

We should expect to find a meridional circulation for most differential rotation profiles, due to the combination of Coriolis and viscous effects, as influenced by the presence or absence of a radial density gradient. The case of a cylindrical rotation profile is special; we discuss this case further below, and in Dikpati (2013b).

In addition to neglecting thermodynamic effects, we omit magnetic fields. It is known that the solar differential rotation varies very little with phase of the magnetic cycle (by less than 0.5 per cent, in the form of torsional oscillations). These properties imply that the feedback of magnetic fields induced by the solar dynamo on the flow (particularly the differential rotation) is relatively small. To address magnetic effects on meridional circulation requires a model that couples the global circulation with the workings of the dynamo, which is beyond the scope of this paper.

There are various possible starting points for developing the axisymmetric equations we solve, including the appendix of Miesch (2005), chapter 5 of Rüdiger (1989) and Rempel (2005a). We first give the three equations for axisymmetric flow, denoting u as the linear rotational velocity, v the velocity in colatitude, w the radial velocity, Ω the rotation rate of the coordinate system, ρo the reference state density, p the pressure gradient relative to a spherically symmetric hydrostatic pressure and Fr, Fθ, Fϕ the force components arising from the Reynolds stress tensor. We linearize these equations of motions, so that terms involving products of the axisymmetric flow components w, v, u are removed. The resulting equations are
(1)
(2)
(3)
in which Fr, Fθ, Fϕ are related to the Reynolds stress tensor Ri, k by
(4)
(5)
(6)
The Reynolds stress tensor is given by
(7)
in which vi are the turbulent velocity components and Eik is the deformation tensor, given in spherical coordinates as
(8)
(9)
(10)
(11)
(12)
and
(13)
Λik represents the non-diffusive parts of the Reynolds stress, which are responsible for driving the differential rotation. In principle, such non-diffusive stresses exist for driving the meridional circulation too, but virtually nothing is known about these, so in mean-field models they are usually set equal to zero (Rüdiger 1989). We will also do that in this initial study, but will revisit this question in a later paper.
Since all the velocities will be small compared to the sound speed, they are constrained by the mass continuity equation
(14)

In contrast to Rempel (2005a), we do use the anelastic approximation, so the time derivative of a perturbation density, found in equation 1 of Rempel (2005a), is not present in equation (14). Since we have also omitted thermodynamics and buoyancy forces from the problem in this first study, the perturbation density does not appear anywhere in the equations.

Equation (14) can be satisfied exactly if we define a streamfunction χ for the meridional velocities w, v such that
(15)
Substituting into equations (4)–(6) for the Reynolds stress tensor, using the deformation tensor components defined in equations (8)–(13), we can write the three force components Fr, Fθ, Fϕ as,
and
Then we can eliminate the perturbation pressure from equations (1) and (2) by multiplying equation (1) by −ρo, equation (2) by rρo, differentiating the modified equation (1) with respect to r, modified equation (2) with respect to θ, and adding the result, followed by substitution for the Reynolds stress and deformation tensors using equations (4)–(13). Then the equations for u and χ can be written as
(16)
in which the function Hϕ includes the non-diffusive Reynolds stresses contained in Λik, and
(17)
in which the Lχ is defined by
(18)

The non-viscous Reynolds Stress terms acting on the meridional flow are represented by |${\cal H}_{\theta }$| and |${\cal H}_r$|⁠, which we set to zero in the present calculation; a forthcoming paper will study the effect of inclusion of |${\cal H}_{\theta }$| and |${\cal H}_r$| deriving from a convection simulation (Guerrero et al. 2013).

Contained within the right-hand sides of equations (16) and (17) are the force balances referred to in Section 1. In particular, the Coriolis terms (first two terms on the right-hand side of equation 16), when balanced against the convergence of the turbulent Reynolds stress that would be contained in |${\cal H}_{\phi }$|⁠, constitute the gyroscopic pumping balance equation. In this balance there is no viscosity, so all of the terms involving the viscosity in (16) would be zero. If the viscous terms are retained, then the meridional circulation would be different than what would be found from the simpler gyroscopic pumping balance equation.

Since the observed solar differential rotation varies so little with time, here we will specify it with guidance from a variety of solar observations. We will solve for the steady-state meridional circulation for several idealized differential rotation profiles, as well as one that corresponds most closely to current observations, to illustrate the range of meridional circulations possible.

A variety of choices of turbulent viscosity are possible. Here we use the simplest form, namely that for isotropic eddy viscosity. In most results, we take this viscosity νt to be independent of latitude, but with a simple depth dependence of the form,
(19)
in which r0, νc, νs and dν are parameters and erf is the error function. We take νc = 7 × 1011 cm2 s−1 and dν = 0.1R and vary νs over a wide range.
For an adiabatic stratification in the convection zone that contains an ideal gas with adiabatic exponent γ = 5/3, a gravitational acceleration of 1/r2 type, we take
(20)
with m = 1.5 and ρb = 0.2 gm cm−3. Note that in order to avoid a singularity at the surface, we have taken (R/r − 0.97) instead of (R/r − 1) on the right-hand side of equation (20). This ensures that the density does not vanish within the computational domain, 0.65RrR and 0° ≤ θ ≤ 90°. Rüdiger (1989, chapter 6) argues that if the turbulent stresses acting on meridional flow are very small scale and relatively uninfluenced by rotation, then their effects can be represented by an isotropic eddy viscosity. If we include no additional forces, then equation (17) can be written as
(21)

We solve equation (21) numerically to seek steady-state solutions for the meridional circulation for various prescribed differential rotations for the Sun. Rempel (2005b) has pointed out that a small change in differential rotation can result in a large change in meridional circulation, at least temporarily, because the small imbalance among large forces in the meridional plane temporarily becomes much larger. Examination of these effects using a time-dependent differential rotation in our model is deferred to a later paper.

SOLUTION METHOD AND BOUNDARY CONDITIONS

The steady-state equation to be solved is
(22)
We employ a relaxation method to solve equation (22), using a multigrid algorithm as implemented in mudpack distributed from the link http://www2.cisl.ucar.edu/resources/legacy/mudpack.
Our computational domain extends from r = rb to r = R in the radial direction, and the North Pole (θ = 0) to the equator (θ = π/2) in θ direction. The following four boundary conditions are applied to solve equation (22): at the pole (θ = 0), χ = 0 for rbrR. At the equator (θ = π/2), χ = 0 for rbrR. At r = R (top boundary) and at r = rb (bottom), for 0 ≤ θ ≤ π/2, we, respectively, apply,
and

RESULTS

We solve equation (22) after making it dimensionless using 1.09 × 1010 cm as the dimensionless length and 1.1 × 108 s as the dimensionless time, the same as used in Dikpati (2011). We present our results in dimensional form, as it is easier for comparing our model output with observations.

For a solar convection zone that includes a fluid density decrease with radius, we compute meridional circulations driven by a wide range of possible differential rotations, some used in well-known solar dynamo calculations, others as instructive limiting cases, and still others that closely follow observations of solar differential rotation. In particular, we calculate meridional circulation for the following four differential rotations: (i) latitudinal differential rotation in the main bulk and a radial shear layer (tachocline) at the base of the convection zone (Brown et al. 1989; Tomczyk, Schou & Thompson 1995; Charbonneau et al. 1997; Dikpati & Charbonneau 1999); (ii) purely latitudinal differential rotation; (iii) a purely radial differential rotation; (iv) helioseismically observed solar-like differential rotation with a tachocline and a near-surface shear layer as well as with a ∼25° slanted rotation contours inside the bulk of the convection zone (Corbard & Thompson 2002).

We omit the calculation of meridional circulation for a cylindrical differential rotation, because it will not produce any meridional circulation for an incompressible fluid, and will produce meridional circulation for a compressible fluid in baroclinic state – the details can be found in Dikpati (2013b).

Latitudinal differential rotation in the main bulk and a radial shear layer at the base of the convection zone

The parametrized form for this differential rotation was derived from helioseismic observations in Charbonneau et al. (1997). This differential rotation has been used in many flux-transport dynamo simulations starting from Dikpati & Charbonneau (1999). We reproduce the mathematical formula as used in Dikpati & Charbonneau (1999) below:
(23)
in which, Ωs(θ) = ΩEq + a2cos 2θ + a4cos 4θ is the surface latitudinal differential rotation, Ωc/2π = 432.8, ΩEq/2π = 460.7, a2/2π = −62.69 and a4/2π = −67.13 nHz.
We do our calculations in this paper in the rotating frame of reference (with respect to the core rotation), for which expression (23) becomes
(24)
All differential rotation profiles below will be given in the rotating frame.

Fig. 1(a) shows the differential rotation profile given by expression (24). For a constant viscosity throughout the convection zone, the meridional circulation produced for this differential rotation in our hydrodynamical model is heavily influenced by the density stratification. Panels (b)–(d) of Fig. 1 show the steady-state meridional circulation solutions for increasing density stratification.

Panel (a) shows colour-filled contours of solar differential rotation parametrized from helioseismic observations by Charbonneau et al. (1997), steady-state meridional circulation patterns produced from our model for density function exponent m = 0.2, 0.5 and 1.5 are presented in panels (b), (c) and (d).
Figure 1.

Panel (a) shows colour-filled contours of solar differential rotation parametrized from helioseismic observations by Charbonneau et al. (1997), steady-state meridional circulation patterns produced from our model for density function exponent m = 0.2, 0.5 and 1.5 are presented in panels (b), (c) and (d).

We use equation (20) to compute the density decrease with radius. By fixing the scaling factor ρb, we pick the values of the exponent m as m = 0.2, 0.5 and 1.5. The value m = 0 would correspond to the incompressible case, and m = 1.5 implies an adiabatically stratified convection zone. The density stratification obtained from the expression (20) for the values m = 0.2 and 0.5 correspond to intermediate cases, which help us see the effect of increasing density stratification on the generation of meridional circulation.

We see in Fig. 1(b) that for m = 0.2, i.e. for the density decrease of only a factor of 2 between bottom and top, we still get a single primary cell with poleward flow near the top, just the same as in the incompressible case [m = 0; see Dikpati (2013a)]. But for m = 0.5, i.e. for a density decrease of a factor of 4, a second, reversed circulation cell appears, for the same Coriolis force per unit mass (see Fig. 1c). For still higher density decreases, this second cell expands in latitude and strengthens relative to the circulation in the primary cell. For the solar-like density decrease for an adiabatic convection zone, we can see in Fig. 1(d) that the boundary between the two circulations has migrated to 55° or so.

How is this evolution possible, if the differential rotation and therefore the Coriolis force driving meridional circulation is the same for all cases? It is possible because it is only the Coriolis force per unit mass that is the same. The Coriolis force per unit volume changes with the change in density profile. In equation (22) for the meridional circulation streamfunction χ, it is the force per unit volume that is doing the driving. In the incompressible case, the effect of Coriolis force is distributed throughout the volume, but with a strong density decrease with radius, that effect is concentrated much more near the bottom of the domain. On the other hand, the radial derivative of the latitudinal component of the Coriolis force per unit volume peaks near the top at high latitudes, because the radial density gradient is very large there. We will examine the details of the Coriolis forcing function in Section 5 below.

Having established the effect of the radial density gradient in determining the meridional flow, we now examine for a solar density profile with radius, how the result is influenced by different profiles of turbulent viscosity with radius. These results are shown in Fig. 2.

Panel (a) shows meridional circulation streamlines for density exponent m = 1.5 for different turbulent viscosity profiles with radius. Panel (b) shows streamlines for the viscosity profile in panel (a) with smallest decrease with depth; and panel (c) shows them for the viscosity with the largest decrease with depth.
Figure 2.

Panel (a) shows meridional circulation streamlines for density exponent m = 1.5 for different turbulent viscosity profiles with radius. Panel (b) shows streamlines for the viscosity profile in panel (a) with smallest decrease with depth; and panel (c) shows them for the viscosity with the largest decrease with depth.

We see in Fig. 2 that the meridional circulation for the viscosity that declines most with depth is almost the same as that for a very small decline, or none at all (see panel d of Fig. 1). Thus, in contrast to the effect of the density decrease with radius, the pattern of meridional circulation is quite insensitive to the viscosity profile. This is fortunate, because we do not know very well how the turbulent viscosity changes with depth within the solar convection zone.

The turbulent viscosity νs at the outer boundary of the domain is a parameter of the problem that can be set. We find that the structure of the meridional flow with latitude and radius is independent of νs; only the amplitude varies. In other words, the streamline pattern in panels (b) and (c) of Fig. 2 apply to all values of νs. Fig. 3(a) shows the peak amplitude of the latitudinal flow as a function of νs; panel (b) displays a family of profiles of meridional flow speed for selected νs in the range shown in panel (a). The thin horizontal line at 25 ms−1 in panel (a) depicts a peak flow speed typical of the solar surface observations. In panel (b), profiles of flow speed for even lower νs than shown have the same shape but larger amplitude.

Panel (a): maximum latitudinal flow speed (negative, poleward) as function of turbulent viscosity νs at outer boundary. Panel (b): latitudinal flow near outer boundary as function of colatitude θ.
Figure 3.

Panel (a): maximum latitudinal flow speed (negative, poleward) as function of turbulent viscosity νs at outer boundary. Panel (b): latitudinal flow near outer boundary as function of colatitude θ.

Two features of Fig. 3 are clearly evident. First, to get peak amplitudes of meridional flow that conform to observations requires we must take a turbulent viscosity near the outer boundary close to 1015 cm2 s−1, at least two orders of magnitude larger than is typically used for the Sun. This is also what Rüdiger et al. (1998) found. Secondly, the polar reversed cell has almost the same peak amplitude as does the primary cell, and the boundary between them is near 45° latitude, a higher relative amplitude and lower latitudinal boundary than observed.

But the main issue with the results in Fig. 3 is the amplitude. It implies that in the Sun there must be some other force, not included in our model, that opposes the meridional flow and thereby reduces its amplitude. We return to this point in Section 5.

Even though the meridional flow speeds for solar values of turbulent viscosity are much too high, for a given profile of viscosity with depth, the structure of the flow (how many cells, the latitude of the boundary between them, velocity profile with depth and latitude) is independent of the turbulent viscosity amplitude. It also varies little with changes in viscosity profile with depth in the lower part of the domain. Therefore, this flow cell structure is a very robust result. This is unlike the Cartesian case studied in Dikpati & Gilman (2012), in which, as the viscosity was reduced, the number of high latitude cells increased. This shows how important it is to include the full spherical geometry and density increase with depth, both of which were absent from the Cartesian model.

Leaving aside issues of flow speed, as opposed to flow profile, we conclude that the differential rotation profile used by Dikpati & Charbonneau (1999) is consistent with a two-celled meridional flow in the convection zone.

Purely latitudinal differential rotation

Another limiting case of interest is that of purely latitudinal differential rotation, set to solar values. This is the same as the Dikpati & Charbonneau (1999) profile, except that there is no tachocline. The form is given in the expression (25) below:
(25)

This differential rotation and the steady-state solution of meridional circulation are shown in Fig. 4. Here we see that we get two cells again. The node between cells is near 45° rather than 55° for the differential rotation with tachocline (see Fig. 1d), and the secondary cell is somewhat stronger than in that case. But qualitatively the meridional circulations are similar. Thus, we infer that the addition of a tachocline changes the meridional circulation quantitatively, but not its basic form. It is tempting to infer that the presence or absence of a tachocline is of secondary importance in determining the number of meridional circulation cells, at least mechanically. From results in Figs 1, 2 and 4, we can infer that the presence or absence of a tachocline is much less important than the density decrease with radius in determining whether there are two circulation cells or just one.

Panel (a) shows colour-filled contours of purely latitudinal differential rotation (same as in Fig. 1(a), but without tachocline). Steady-state meridional circulation pattern produced is presented in panel (b).
Figure 4.

Panel (a) shows colour-filled contours of purely latitudinal differential rotation (same as in Fig. 1(a), but without tachocline). Steady-state meridional circulation pattern produced is presented in panel (b).

Purely radial differential rotation

Another idealized case of interest is that of a purely radial differential rotation. Fig. 5 shows one such case, which was computed by setting a2 and a4 equal to zero in Ωs and taking d1 in equation (24) to be large enough that the differential rotation is linear in r. We see that the resulting meridional circulation is a single cell with poleward flow near the outer boundary. This is what we should expect, since the radial component of Coriolis force is outward at all latitudes, but strongest in low latitudes. Since a purely latitudinal differential rotation always yields two meridional circulation cells, clearly if the differential rotation can be thought of as a weighted combination of latitudinal and radial differential rotation, the resulting meridional circulation will have one cell if the radial gradient dominates, and two cells if the latitude gradient dominates.

Panel (a) shows colour-filled contours of purely radial differential rotation. Steady-state meridional circulation pattern produced is presented in panel (b).
Figure 5.

Panel (a) shows colour-filled contours of purely radial differential rotation. Steady-state meridional circulation pattern produced is presented in panel (b).

Differential rotation with surface shear layer, tilted rotation contours and tachocline

Since the 1970s, there have been indication of a near-surface shear layer inferred from various observations. A detailed differential rotation pattern was derived by Corbard & Thompson (2002) from helioseismic observations, which contains (i) a near-surface radial shear layer, (ii) ∼25° tilted rotation contours in the main bulk of the convection zone and (iii) a tachocline at the bottom. We use the mathematical formula for this differential rotation following Dikpati et al. (2002), which is (in rotating frame) given by
(26)
in which
and
Here ΩEq, Ωc, a2 and a4 have slightly different values compared to that noted in Section 3.1; ΩEq/2π = 452.5 nHz, Ωc/2π = 435 nHz, a2/2π = −61 nHz and a4/2π = −73.5 nHz. What the other parameters denote and their values are as follows: the location at which the slope of the rotation contours starts, rcz = 0.71R; the equatorial rotation rate at rcz, Ωcz/2π = Ωca2/5 − 9a4/105 = 453.5 nHz; the width, ωcz = 0.05R; the position of the tachocline, rc = 0.69R; width of the tachocline, ωc = 0.05R; the position of the maximum of near-surface shear, rs = 0.97R; the width, ωs = 0.05R, and the near-surface shear is defined by β0/2π = 437 nHz, β3/2π = −214 nHz and β6/2π = −503 nHz. This produces the profile shown in fig. 6 of Corbard & Thompson (2002). Based on helioseismic observations, this is the most realistic case for the Sun.
Panel (a): differential rotation contours including near-surface shear layer, as shown in the expression (26). Panel (b): streamfunction for meridional flow driven by Coriolis forces associated with differential rotation profile shown in panel (a). Panel (c): latitudinal flow speed near outer boundary for streamfunction shown in panel (b) (solid curve) and for case with near-surface shear layer removed (dotted curve).
Figure 6.

Panel (a): differential rotation contours including near-surface shear layer, as shown in the expression (26). Panel (b): streamfunction for meridional flow driven by Coriolis forces associated with differential rotation profile shown in panel (a). Panel (c): latitudinal flow speed near outer boundary for streamfunction shown in panel (b) (solid curve) and for case with near-surface shear layer removed (dotted curve).

The results for this case are shown in Fig. 6. Panel (a) gives the differential rotation contours; the near-surface shear layer is clearly visible. Panel (b) shows the resulting meridional streamlines driven by Coriolis forces from this differential rotation and turbulent Reynolds stresses. Panel (c) displays the latitudinal flow at the outer surface as a function of colatitude, for two cases, one with and one without the near-surface shear layer.

We see in Fig. 6 that the two-cell pattern prevails, both with and without the near-surface shear layer. The main difference between the two cases is that with the near-surface shear, the primary cell with poleward flow at the outer surface is somewhat stronger than the counter cell at high latitudes, compared to the case without near-surface shear. This is a quantitative improvement when compared to observations of surface meridional flow.

We conclude that meridional circulation driven by Coriolis forces from the most realistic solar differential rotation profile known favours a two-celled circulation, with a stronger primary cell in low and mid-latitudes, and a somewhat weaker counter cell at high latitudes, with the boundary between cells near 55°. To deviate strongly from this pattern would require additional physical processes that could compete with the Coriolis forces.

Meridional circulation from antisolar differential rotation profiles

For every differential rotation profile used in Sections 4.1–4.4, there is a corresponding antisolar differential rotation profile which would also drive meridional circulation. Leaving aside the near-surface shear layer for the moment, we reason that antisolar differential rotation profiles are specified by keeping the rotation of the core fixed, and changing the sign of the coefficients that represent all departures from that core rate. Without the near-surface shear layer, from equation (22) it is clear by inspection that for antisolar differential rotation, because the forcing function changes sign, the meridional circulation streamfunction χ simply changes sign, implying the opposite sense of circulation. Thus, all primary cells would have equatorward flow near the outer surface. Fig. 7 gives an example of an antisolar differential rotation and the meridional circulation produced in this case.

Panel (a): antisolar differential rotation profile produced from profile in Fig. 1(a), found by reversing sign of all rotation parameters in Ωs in equation (24). Panel (b): streamlines of meridional flow driven by Coriolis forces associated with differential rotation profile shown in panel (a); note the pattern is the same as in Fig. 1(a), but with the flow direction reversed.
Figure 7.

Panel (a): antisolar differential rotation profile produced from profile in Fig. 1(a), found by reversing sign of all rotation parameters in Ωs in equation (24). Panel (b): streamlines of meridional flow driven by Coriolis forces associated with differential rotation profile shown in panel (a); note the pattern is the same as in Fig. 1(a), but with the flow direction reversed.

The surface shear layer is different, in that it is formed by the influence of the overall rotation of the system, which has not changed sign, so one should take that profile as it was for the solar case. Then the resulting meridional circulation will not be a simple change of sign. We have not displayed that case here. There is of course an even more general family of antisolar differential rotation profiles, for stars that rotate in the opposite direction to the Sun. In that case, the near-surface shear layer would be reversed in sign, just as would all the other components of the differential rotation.

INTERPRETATION AND DISCUSSION

It is clear from all of the above results that, when Coriolis forces associated with differential rotation are important, the meridional circulation that is produced by these forces is sensitive quantitatively to the profile of differential rotation. When there is a substantial decrease in fluid density with radius, such as in the solar convection zone, a two-celled meridional circulation pattern, with reversed flow in high latitudes, virtually always results. It follows that when we know the profile of differential rotation rather well, as we do for the Sun, the meridional circulation pattern to be expected is fairly strongly constrained. Unless forces not included in our model act to change the number of meridional circulation cells with latitude, we should expect the two-celled pattern to be the norm. It is of interest to note that Tassoul & Tassoul (1995) also found two meridional circulation cells in the radiative envelopes of young stars, in that case driven in part by thermal effects but strongly influenced by Coriolis forces.

One way to understand the above results is to look more closely at the forcing function that comes from the Coriolis forces in the latitudinal and radial equations of motion. It is important to note that these forcing functions do involve the radial gradient of density, since it is the Coriolis force per unitvolume that is being differentiated to form the forcing function. In Fig. 8, we display the two parts of this forcing function for the differential rotation profile shown in expression (24), for the incompressible case (panels a and b) and the solar-like m = 1.5 compressible case (panels c and d). Panels (a) and (c) give the part that involves the radial derivative of the differential rotation, which comes from the latitudinal component of the Coriolis force; panels (b) and (d) show the part arising from the latitudinal derivative of differential rotation, which comes from the radial component of Coriolis force. The colour scales are the same for panels (a) and (b), and separately also the same for panels (c) and (d), so they can be directly compared in determining the total forcing function in incompressible and compressible cases. Blue and green shaded areas give rise to counterclockwise circulation, while red and orange areas lead to clockwise circulation. These circulation tendencies represent the action of the Coriolis forces in the meridian plane to create vorticity in the longitudinal direction. For a given differential rotation profile, these forcing function patterns are independent of the turbulent viscosity used.

Components of forcing function due to Coriolis forces for differential rotation profile shown in Fig. 1a. Panels (a) and (b) are for the incompressible case, panels (c) and (d) for the compressible case. Panels (a) and (c): component from latitudinal Coriolis force (term differentiated with respect to r); Panels (b) and (d): component from radial Coriolis force (term differentiated with respect to θ).
Figure 8.

Components of forcing function due to Coriolis forces for differential rotation profile shown in Fig. 1a. Panels (a) and (b) are for the incompressible case, panels (c) and (d) for the compressible case. Panels (a) and (c): component from latitudinal Coriolis force (term differentiated with respect to r); Panels (b) and (d): component from radial Coriolis force (term differentiated with respect to θ).

We see from Fig. 8 that the radial component of Coriolis force dominates in low latitudes for both cases, and is always trying to produce the primary cell with poleward flow near the outer boundary, so it is no surprise that the primary flow cell has poleward flow near the outer boundary no matter what the density decrease with radius is.

By contrast, the two cases differ substantially at high latitudes. For the incompressible case, in high latitudes the latitude component of Coriolis force is still trying to produce an equatorward flow near the bottom (the tightly packed blue contours there) which in effect extends the primary cell all the way to the pole. Near the top at high latitudes the two Coriolis force components almost cancel each other out, so the strong forcing near the bottom dominates at all depths. The net result of these forcings is that in the incompressible case the primary cell extends all the way to the pole.

In the compressible case, at high latitudes the situation is nearly reversed: red contours dominate except very near the bottom, where red and blue cancel out. Here the Coriolis forces in the upper part of the shell are dominating in determining the net circulation. The radial gradient of density, large near the top, is coupled with the latitudinal component of Coriolis force to produce a longitudinal component of vorticity with the opposite swirl to that of the primary cell. In other words a reversed cell is formed, with equatorward flow near the outer boundary. This effect comes from the ‘deformation’ of fluid elements, because in this mechanical model the density and pressure surfaces do not coincide [see equation 3.13 of Rüdiger & Hollerbach (2004) and adjacent text].

These results quantify the reasoning that went into fig. 2 of Dikpati & Gilman (2012) as a way to picture why we should get a primary and a secondary reversed cell in the Sun. Figures similar to Fig. 8 can be constructed for all other differential rotation profiles we have studied, to explain the resulting meridional flow in each case. We conclude that these forcing functions from the Coriolis forces per unit volume in the meridian plane control the structure of the meridional flow, even when other physical effects are included.

The results presented above differ significantly from those found previously by Dikpati & Gilman (2012), in which polar regions were ‘unwound’ into the form of a Cartesian channel, in which the fluid density was constant with depth and the meridional circulation was forced from the equatorward side of the channel. There much higher turbulent viscosities were needed to obtain a single cell in the channel, or even a primary cell with just one secondary reversed cell near the polar side of the channel. The reason is that when the polar domain is greatly stretched to produce a straight channel, the effect of a particular viscosity value becomes much smaller – the longitude dimension becomes much bigger. This allows meridional flow with multiple nodes across the channel to be excited, and also creates the phenomenon of node merging. In spherical geometry, this finer flow structure is not favoured, because in spherical geometry the physical distances get shorter and shorter as the pole is approached, and the damping effects of viscosity are correspondingly larger. This effect minimizes the number of nodes possible in the meridional circulation, and keeps the single node that occurs in most cases from getting very close to the poles. A similar effect should be at work in the polar regions of the solar convection zone.

Since in this mechanical model the fluid pressure in the solutions varies with both latitude and radius, while the density is only a function of radius, when we add thermodynamics to the problem we should expect ‘baroclinic’ effects, because the pressure, density, temperature and specific entropy surfaces would not coincide. The structure of such surfaces found, when thermodynamics is included, could be compared with those from 3D global convection simulations that generate solar-like differential rotation and meridional circulation.

Similarly, profiles of turbulent Reynolds stresses that transport angular momentum that are consistent with the observed differential rotation and calculated meridional circulation could be compared with stresses generated by the 3D global convection models. In this way, it can be determined how well the balance equations for the meridional plane and the zonal direction described earlier agree with the outputs from 3D global convection simulations. In other words, these balance equations can be used as diagnostics tools for analysing the much more complex 3D model simulations.

It is very important to distinguish the meridional circulation patterns driven by Coriolis forces we have found here from axisymmetric convective rings that could be excited in full 3D global convection models and even axisymmetric mean-field models. Rempel (2005a) and we here are able to prevent the axisymmetric convective rings from being excited, and are therefore able to isolate the meridional circulation from such effects. But no full 3D global convection model, e.g., Gilman & Miller (1986), is able to do this, so the total meridional flow found in these models is likely to contain more cells in latitude, and perhaps stacked in depth, particularly in lower latitudes. Some influence of rotation on the orientation of these rings in the meridional plane is likely. Schematically, the different patterns are depicted in Fig. 9. The turbulent Reynolds stresses and Coriolis force from differential rotation provide a relatively steady driving for meridional circulation, while convective driving, by virtue of its nature, will be very unsteady. We infer that it will be hard to find persistent one- or two-celled global meridional circulation for even a few months, if the Sun's meridional circulation is primarily convectively driven.

Schematic diagram showing distinction between forms of global meridional circulation (red streamlines) and axisymmetric convective rings (blue). Global meridional circulation is driven by Coriolis forces from differential rotation, Reynolds stresses and buoyancy forces arising from thermodynamics, whereas axisymmetric convective rings are driven by buoyancy forces. Orientation of axisymmetric convective rings is influenced by Ω; they will generate Coriolis forces that will modify differential rotation.
Figure 9.

Schematic diagram showing distinction between forms of global meridional circulation (red streamlines) and axisymmetric convective rings (blue). Global meridional circulation is driven by Coriolis forces from differential rotation, Reynolds stresses and buoyancy forces arising from thermodynamics, whereas axisymmetric convective rings are driven by buoyancy forces. Orientation of axisymmetric convective rings is influenced by Ω; they will generate Coriolis forces that will modify differential rotation.

The overall energetics of the system for meridional flow we have presented here follows closely the energy flow diagram in fig. 4 of Rempel (2005b), except that here we have specified the differential rotation, so its driving by the Λ effect and braking by viscous loss is implicit, not explicit. Rempel (2005b) energy diagram also says that the buoyancy force must brake the meridional circulation that is driven by the Coriolis forces from differential rotation. Viscous loss by itself is not enough to keep the meridional flow speeds within the range of observations.

We can also argue strongly against a convective or even thermodynamic origin for the meridional circulation from the fact that we have shown that, if Coriolis forces from differential rotation are opposed only by turbulent viscous forces, the resulting meridional circulation is larger than observed by a factor of 10–100, even though the profile of meridional flow compares well with surface observations. This means that a buoyancy force, which is the only way thermodynamics can directly come into the equations of motion, must always oppose the meridional circulation, rather than drive it. In the presence of thermodynamics, in whatever form taken, the Coriolis force from the differential rotation is still present. This point is well recognized in the energy flow diagram shown in Rempel (2005b). Since without an opposing buoyancy force, the meridional circulation is 10–100 times too large, the buoyancy force must come quite close to balancing the Coriolis force everywhere in the domain. This leads to much lower meridional flow speeds when solar-like turbulent viscosities are used in the model. But the balance cannot be exact, or there would be no meridional flow at all. In that case, the differential rotation would be a ‘thermal wind’. This closeness can be used to estimate the latitude-dependent density structure of a nearly adiabatically stratified convection zone for a particular differential rotation, from which the buoyancy force that opposes meridional flow can be calculated. But this is the subject of future papers.

SUMMARY

We have shown that our purely mechanical model yields a two-celled circulation for a solar-like differential rotation and turbulent Reynolds stresses. This result is independent of the value of turbulent viscosity used, and not very dependent on its profile with depth. Instead, it is due largely to the density decline with radius in the convection zone.

The meridional circulation amplitude we find is much too large for the Sun when solar-type turbulent viscosities are used, so there must be another process working against the meridional flow that is driven by the turbulent Reynolds stresses and the Coriolis forces from the differential rotation. There are at least three possibilities, all processes that have been studied separately in the past, as reviewed in Section 1. Any process that itself tends to drive meridional flow that is towards the equator at the outer surface is a candidate. These include anisotropy in the small-scale turbulent momentum transport, such as in Köhler (1970), latitudinal gradients in radial heat flux, such as in Durney & Roxburgh (1971) and negative radial buoyancy forces, as in Rempel (2005a). These additional forces may therefore increase the probability of the persistence of a secondary reverse cell in high latitudes.

The problem is that, unlike the Coriolis forces associated with differential rotation, the magnitude of each of these other processes is relatively poorly known. Nevertheless, starting from the present model, each of these effects could be included, and the amplitude of each needed to reduce the meridional circulation to observed values could be determined. This may be the only feasible way to estimate what their amplitudes are. This effort would require new calculations beyond what we have done here; the approach taken here provides a quantitative modelling framework within which to carry out such studies in the future. More detailed models would be needed to include all the relevant processes. However, assuming a detailed differential rotation profile in the solar problem implies a particular profile of Reynolds stresses and thermal effects, so they are included implicitly.

A remaining open question is what was different in the Sun during cycle 23 that caused the poleward meridional flow to reach all the way to the poles, while for all other recent cycles the primary cell reached only to about 60°, at least in the MWO data. One possibility is that there are domains of slightly subadiabatic stratification that the meridional circulation advects to give a slightly stronger or weaker negative buoyancy force, perhaps in high latitudes, resulting in more or less suppression of the secondary cell. Recent simulations by Masada, Yamada & Kageyama (2013) show that including a subadiabatic layer below the base of the convection zone changes the form of meridional circulation from predominantly one cell in latitude to two cells. This could also happen with weakly subadiabatic stratification within the lower half of the convection zone itself. However, we note that in figs 6(c) and (d) of Masada et al. (2013) we may be seeing a combination of axisymmetric convective rolls and meridional circulation driven by the Coriolis forces from the differential rotation, because these calculations, like Gilman & Miller (1986), are also done for substantially supercritical Rayleigh numbers.

Another possibility is that varying magnetic fields induced by the solar dynamo modulate the meridional flow in a different way in one cycle compared to another. To assess this possibility requires an magnetohydrodynamic theory for solar meridional circulation, which does not yet exist. Future efforts are needed to address the role of magnetic fields in determining the meridional flow and its variations in the Sun.

We thank Peter Gilman for many helpful discussions on the theory of meridional circulation. Our thanks go to Gustavo Guerrero for fruitful discussions and for his suggestion of computing the antisolar case. This work is partially supported by NASA's LWS grant with award number NNX08AQ34G. The National Center for Atmospheric Research is sponsored by the National Science Foundation.

REFERENCES

Belvedere
G.
Paterno
L.
Stix
M.
Geophys. Astrophys. Fluid Dyn.
,
1980
, vol.
14
pg.
209
Brown
T. M.
Christensen-Dalsgaard
J.
Dziembowski
W. A.
Goode
P. R.
Gough
D. O.
Morrow
C. A.
ApJ
,
1989
, vol.
343
pg.
526
Charbonneau
P.
Christensen-Dalsgaard
J.
Henning
R.
Schou
J.
Thompson
M. J.
Tomczyk
S.
Provost
J.
Schmider
F.-X.
Proc. IAU Symp. 181, Sounding Solar and Stellar Interiors. Kluwer, Dordrecht
,
1997
pg.
161
Corbard
T.
Thompson
M. J.
Sol. Phys.
,
2002
, vol.
205
pg.
211
Dikpati
M.
ApJ
,
2011
, vol.
733
pg.
90
Dikpati
M.
Jain
K.
Tripathy
S. C.
Hill
F.
Leibacher
J. W.
Pevtsov
A. A.
ASP Conf. Ser. Vol. 478
,
2013a
San Fransisco
Fifty Years of Seismology of the Sun and Stars. Astron. Soc. Pac.
pg.
277
Dikpati
M.
Geophys. Astrophys. Fluid Dyn.
,
2013b
, vol.
478
pg.
277
Dikpati
M.
Charbonneau
P.
ApJ
,
1999
, vol.
518
pg.
508
Dikpati
M.
Gilman
P. A.
ApJ
,
2012
, vol.
746
pg.
65
Dikpati
M.
Corbard
T.
Thompson
M. J.
Gilman
P. A.
ApJ
,
2002
, vol.
575
pg.
L41
Dikpati
M.
Gilman
P. A.
de Toma
G.
Ulrich
R. K.
Geophys. Res. Lett.
,
2010
, vol.
37
pg.
L14107
Durney
B. R.
Sol. Phys
,
1996
, vol.
169
pg.
1
Durney
B. R.
Roxburgh
I. W.
Sol. Phys
,
1971
, vol.
16
pg.
3
Gilman
P. A.
Miller
J.
ApJS
,
1986
, vol.
61
pg.
585
Guerrero
G.
Smolarkiewicz
P. K.
Kosovichev
A. G.
Mansour
N. N.
ApJ
,
2013
, vol.
779
pg.
176
Hathaway
D. H.
Rightmire
L.
ApJ
,
2011
, vol.
729
pg.
80
Kitchatinov
L. L.
Olemskoy
S. V.
MNRAS
,
2011
, vol.
411
pg.
1059
Köhler
H.
Sol. Phys.
,
1970
, vol.
13
pg.
3
Komm
R.
Gonzalez-Hernandez
I.
Hill
F.
Bogart
R.
Rabello-Soares
M. C.
Haber
D.
Sol. Phys.
,
2013
, vol.
287
pg.
85
Lorenz
E. N.
The Nature and Theory of the General Circulation of the Atmosphere. World Meteorological Organization, Geneva
,
1967
pg.
115
Masada
Y.
Yamada
K.
Kageyama
A.
ApJ
,
2013
, vol.
778
pg.
1
Miesch
M. S.
Living Rev. Sol. Phys.
,
2005
, vol.
2
pg.
1
Pidatella
R. M.
Stix
M.
Belvedere
G.
Paterno
L.
A&A
,
1986
, vol.
156
pg.
22
Rempel
M.
ApJ
,
2005a
, vol.
622
pg.
1320
Rempel
M.
ApJ
,
2005b
, vol.
631
pg.
1286
Rightmire-Upton
L.
Hathaway
D. H.
Kosak
K.
ApJ
,
2012
, vol.
761
pg.
L14
Rüdiger
G.
Differential Rotation and Stellar Convection
,
1989
Berlin
Akademie-Verlag
pg.
328
Rüdiger
G.
Hollerbach
R.
The Magnetic Universe
,
2004
Weinheim
Wiley
pg.
332
Rüdiger
G.
von Rekowski
B.
Donahue
R. A.
Baliunas
S. L.
ApJ
,
1998
, vol.
494
pg.
691
Stix
M.
Rev. Mod. Astron.
,
1989
, vol.
2
pg.
248
Tassoul
M.
Tassoul
J.-L.
ApJ
,
1995
, vol.
440
pg.
789
Tomczyk
S.
Schou
J.
Thompson
M. J.
ApJ
,
1995
, vol.
448
pg.
L57
Ulrich
R. K.
ApJ
,
2010
, vol.
725
pg.
658
Zhao
J.
Bogart
R. S.
Kosovichev
A. G.
Duvall
T. L.
Hartlep
T.
ApJ
,
2013
, vol.
774
pg.
L29