# Evidence for Large-Scale Subsurface Convection in the Sun

###### Abstract

A helioseismic statistical waveform analysis of subsurface flow was performed on two 720-day time series of SOHO/MDI medium- spherical-harmonic coefficients. The time series coincide with epochs of high and low solar activity. Time-dependent coupling-strength coefficients of modes of the same radial order and degree , but different azimuthal order , were inferred from the waveform analysis. These coefficients are sensitive to flows and general aspherical structure. For odd values of , the coefficient measures an average over depth of the amplitude of one spherical-harmonic component of the toroidal flow velocity field. The depth-dependent weighting function defining the average velocity is the fractional kinetic energy density in radius of modes of the multiplet. A mean-square -dependent flow velocity was inferred from the -coefficients for in the range through for each and in the respective ranges through and through for the epochs of high and low activity. A further averaging, over , yielded a root-mean-square flow velocity as a function of for each epoch, which average increases from about at to at . The inferred velocities are consistent with (though perhaps do not demand) a cellular pattern of flow extending over the vertical range of mode sensitivity, estimated to be about four percent of the solar radius below the photosphere.

###### keywords:

Sun:oscillations – Sun:helioseismology – stars:interiors – convection^{†}

^{†}pubyear: 2016

^{†}

^{†}pagerange: Evidence for Large-Scale Subsurface Convection in the Sun–Evidence for Large-Scale Subsurface Convection in the Sun

## 1 Introduction

Turbulent motions in a sun-like star are believed to play a leading role in convecting the star’s luminosity through the outer portions of its envelope. Turbulent convection is also believed to participate, along with meridional circulation, in distributing fluid angular momentum within stellar convection zones, thereby setting up the differential rotation (or ‘angular velocity’) within the star (e.g. Hansen et al., 2004; Stix, 2004). A great deal of effort has been directed toward understanding the interactions of the Sun’s mass flows, which are observed over a great range of spatial scales (Kitchatinov & Rüdiger, 1999; Durney, 2000; De Rosa et al., 2002; Küker & Rüdiger, 2005; Miesch, 2007). The most stringent observational constraints on the internal angular velocity are provided by whole-Sun oscillation-frequency-splitting measurements (Thompson et al., 2003; Howe, 2009). For other large-scale subsurface flows, the methods of local helioseismology are used (Gizon & Birch, 2005). Stellar structure and evolution calculations have traditionally relied on the so-called mixing-length theory of convection (Böhm-Vitense, 1958). In the last few decades, however, numerical simulations of convection-zone-scale flows have begun to provide a detailed theoretical picture of the dynamic and magnetic interiors of stars in general and in particular of the Sun (e.g. Miesch, 2005).

Helioseismic measurements are beginning to reveal deep-seated large-scale turbulent motion in the Sun. Unlike the differential rotation and the near-surface meridional flow (Giles et al., 1997; Braun & Fan, 1998; González Hernández et al., 1999; Giles, 2000; Haber et al., 2002; Hughes & Thompson, 2003; Zhao & Kosovichev, 2004; Chou & Ladenkov, 2005), the largest turbulent scales have yet to be well, or even consistently, characterized. Large-scale turbulent motions are seen in the first Mm below the photosphere (Hathaway et al., 2000; Featherstone et al., 2006). Although there does not appear to be an excess of velocity power at “giant-cell” scales, i.e. horizontal scales comparable to the Mm depth of the convection zone, there is evidence of elongated and persistent large-scale flow structures. Existing helioseismic observations of convection at greater depths are hard to reconcile with one another. Hanasoge et al. (2012), using time-distance analysis, report only upper limits on flow speeds at and , which limits are at least an order of magnitude smaller than the speeds seen in simulations (which more or less agree with mixing-length theory). The disparity between theory and observation becomes apparent on horizontal scales somewhat larger than the Mm supergranulation scale (corresponding to angular wavenumber ) and increases dramatically with increasing horizontal scale. The upper limits therefore challenge the conventional picture of convective energy transport in stars. However, the time-distance results seem to contradict the more recent ring-diagram measurements of Greer et al. (2015), which are in closer agreement with simulations and in fact show the amplitude of turbulent motions increasing with depth at depths approaching . Supergranulation tracking measurements (Hathaway et al., 2013) suggest large-scale convective motions as deep as Mm, with amplitudes similar to those seen by Greer et al.

This paper describes a helioseismic analysis of large-scale solar subsurface turbulence based on years of medium- Doppler images from the Michelson Doppler Imager (MDI) on the SOHO spacecraft (Scherrer et al., 1995). The analysis of these data takes a global-mode, statistical waveform approach, in which time-dependent mode-coupling strengths are first inferred from spectral-domain covariance data. The coupling strengths are then used to estimate the power in subsurface turbulence over the range of angular wavenumber.

The data sensitivity model used in the analysis is detailed in Section 2 in the context of a more general sensitivity model. The analysis of MDI time series resulting in flow velocity measurements is presented in Section 3, together with analogous measurements of near-surface magnetic activity. Section 4 summarizes the findings and discusses their implications and the prospects for detecting yet deeper flows.

## 2 The Forward Model

### 2.1 A general data sensitivity model for helioseismic analysis of flows

In helioseismic statistical waveform analysis, simple products of the observed solar oscillation signal are inverted for the subsurface mass flow velocity and aspherical structure. In waveform analysis and in other seismic analysis, rigorous forward modeling is crucial to the accurate retrieval of subsurface conditions. The data analysis described in the following section used oscillation signal in the form of coefficients in an expansion of the observed photospheric Doppler velocity in scalar spherical harmonic functions, , of heliographic colatitude, , and longitude, , and in sinusoidal functions, , of time, . The analysis uses covariance data , where ‘’ means conjugation, for many combinations of and . The forward model specifies the dependence of on the time-varying interior flow velocity, , with denoting spherical polar coordinates , , and and denoting statistical expectation. In view of the weakness of the flows being investigated, the expectation is taken to be the sum of a zeroth-order term , describing the effect of a spherical reference Sun, and a perturbation describing the effect of large-scale turbulent flows.

To model the effect of deep-seated physical perturbations on waves observed in the photosphere one needs to consider the behavior of the wave field throughout a substantial volume of the Sun. For the present analysis the wave field is represented by frequency-domain mode amplitudes, , defined by the expansion

(1) |

where is the component at frequency and at position in the solar interior of the Lagrangian wave displacement and are the displacement eigenfunctions of the normal modes of the spherical reference model. The index refers to the usual () indices of global oscillation modes. (In this paper the convention is used for the Fourier component of a time series of finite duration.)

The present data analysis uses approximate expressions for the zeroth-order (unperturbed) wave-field covariances and the flow-induced covariance perturbations that are analogous to expressions derived for local helioseismic analysis (Woodard, 2006). In zeroth order, and for modes that are independently excited, the expectation is zero unless and and the mode-amplitude spectrum has the lorentzian profile

(2) |

The function

(3) |

represents the complex response of a simple damped oscillator to a unit harmonic driving force of frequency . The parameters and are the resonant frequency and damping rate of mode and is a normalization factor. For the perturbation, one has

(4) |

where the complex-valued coupling coefficient, , describes the effect of subsurface flows and aspherical structure. The expression

(5) |

(e.g. Woodard, 2014), in which denotes an element of mass and the integration is carried out over the entire Sun, describes the first-order effect of the flow velocity. Mode-coupling effects of various structural asphericities have been quantified in Lavely & Ritzwoller (1992). Magnetic activity affects stellar oscillations (Gizon et al., 2008) and would be expected to contribute to the coupling coefficients. The Hermitian property was assumed in obtaining the above expression for . This should be a good approximation, provided that the divergence of the turbulent mass flux can be neglected in the solar interior and that the vertical flux can be ignored at the outer turning points of the observed waves.

Because present-day seismic observations sample only the Earth-facing hemisphere of the Sun, the mode amplitudes are not independently observable. More precisely, the response of the observed signal to the amplitudes is given approximately by

(6) |

where the leakage matrix is described, for instance, in Schou & Brown (1994) and is a background signal. Expressions for and follow straightforwardly from the above leakage relation and from Equations (2) and (4).

As in Woodard (2014, hereafter W14), the mode couplings can be expanded

(7) |

where, apart from sign factors, the are Clebsch-Gordon coefficients and . The usefulness of this expansion stems from the fact that is sensitive to the vector spherical-harmonic component of degree and azimuthal order of the flow velocity field at frequency .

### 2.2 Modeling flow-dependent couplings of modes within multiplets

The presence of the resonant factors and in the expression for wave covariance sensitivity (Equations (4), (2), and (3)) implies that the data can be particularly sensitive to the flow velocity when and are close to and , respectively, for specified and . Therefore, in view of Equation (5), these covariance data should be especially sensitive to turbulent frequencies of order . Giant-cell patterns seen in observations (Hathaway et al., 2015) and in numerical simulations of solar convection (Miesch et al., 2008) have lifetimes of weeks to months. These patterns are thus expected to produce relatively strong couplings between oscillation modes of frequency differing by less than , in particular between nearly- degenerate modes of the same () multiplet.

The present data analysis is based on measurements of (near-resonant) couplings of modes of the same and . Replacing the indices , , and by (), (), and , respectively, in Equation (4) and using Equation (7) one obtains an explicit expression for the wave-covariance sensitivity to the parameters, for any pair of oscillation modes. Replacing by and by in the resulting expression and suppressing the multiplet indices , one obtains the sensitivity relation

(8) |

for modes of one multiplet, where the sensitivity kernel is given by

(9) |

The present analysis uses the approximate expressions given in W14 for the sensitivity of the -coefficients to the subsurface velocity appropriate in the limit . In the low- approximation, the -coefficients are sensitive to only the odd- toroidal components of the flow velocity. For either odd or even , the toroidal component of spherical indices () has the form

(10) |

(Equation 25 of Lavely & Ritzwoller, 1992), where (, , ) are the unit vectors of the spherical-polar coordinates () and is the corresponding scalar spherical harmonic function. As in W14, it is convenient to use the scaled radial function

(11) |

As discussed in that paper, the scaled function contributes to the mean square flow velocity at radius .

Manipulation of Equations 35 through 47 of W14 gives

(12) |

for odd , where is an average over of , defined above, and . The -dependent weighting function used to obtain is the fractional kinetic energy per unit of modes of the () multiplet.

## 3 Data Analysis and Results

Although the covariance data can be inverted directly for the flow velocity, the approach taken here was more indirect in that -coefficients were first extracted from the data, by a least-squares fitting procedure, then subjected to further analysis. The -coefficients isolate vector spherical-harmonic and frequency contributions of the velocity field, thus the analysis was similar in spirit to a multi-channel deconvolution (MCD, Jensen et al., 1998).

The signal components were computed for two -day time series of MDI medium- spherical-harmonic coefficients, for in the range . The midpoints of the time series occur in mid 1997 and late 2000, during epochs of high and low solar magnetic activity. By Equation (6) and due to the narrowness of the global oscillation frequency profiles, can be used as a proxy for the mode amplitude for close to the mode frequency . More precisely, and the sensitivity of to the -coefficients is obtained by substituting ‘’ for ‘’ and ‘’ for ‘’, where , in Equation (8).

The form of the sensitivity kernel depends on parameters, such as the frequency and line width, defining the frequency profiles of global oscillation modes. The mode parameters used for this study were obtained from a recent re-analysis of MDI spherical-harmonic power spectra (Larson & Schou, 2015). The -dependence of the mode frequencies reflects the Sun’s latitude- and radius-dependent angular velocity referred to the Sun-orbiting frame of the SOHO spacecraft. The frequencies, , of turbulent components also apply to this frame.

The perturbation is the expectation of ‘residual’ covariance data and the -parameters were obtained from a simple linear least-squares fit of the forward model to these data. The zeroth-order expectation of is obtained from Equations (2) and (6). In practice, due to the finite line width of the modes and the fact that the mode frequencies depend only weakly on , a given mode signal leaks into more than one component of the observed signal. A similar leakage, of the photospheric granulation signal and perhaps other signals of solar origin, is expected to occur (via the background term of Equation (6)). The observed background power is important at low frequency. As a consequence, differs substantially from zero when is a small integer. Therefore, in the interest of simplicity, the fitting procedure used only covariances for which . This restriction precludes the detection of flow components of degree less than .

The -parameters for a given and were obtained one at a time from the covariance data according to

(13) |

for ranging between and , subject to the condition , and for lying within one full width of the centroids, , of the (lorentzian) mode profiles. is the sensitivity kernel described earlier in this Section. The above least-squares estimator is analogous to the one given by Equation 80 of Woodard (2006). Parameter estimates were obtained for through and for through at each . These ranges are optimized for sensitivity to near-sectoral-harmonic components of giant-cell-scale flows. The were estimated for within about of the expected (rotational-advection) frequency of flow-velocity features of azimuthal order . The sampled frequency range was intended to be sufficiently wide to capture the time-scales of rotationally-advected giant-cell-scale flow patterns. For both the low- and high-activity epochs the analysis was performed on modes of through and through . These modes are collectively sensitive to the subsurface flow velocity to a depth of about of the solar radius below the photosphere (e.g., Figure 5 of Christensen-Dalsgaard, 2002).

The sensitivity relation (12) provides a basis for inverting the odd- -coefficients for the toroidal velocity profiles . But as it was not clear at the outset whether large-scale convective motions would even be detectable, the next step in the analysis was simply to establish whether or not a flow signature is actually present in the measured -coefficients. To this end, unweighted averages, over at fixed , of the estimated signed coefficients were computed, on the assumption that their expected values vary slowly with . Given that the sensitivity functions of the signed -coefficients are, by Equation (12), slowly varying functions of , this assumption is reasonable provided that the profiles do not vary too rapidly with in the layers explored by the modes analyzed. Denoting the -averaged signed -coefficients by , and in view of the relatively narrow range of used in the analysis, one obtains, from Equation (12), the expression

(14) |

where is a typical value of in the measured range.

While the averaging procedure does suppress measurement noise, the scatter of the individual samples suggests that the signal-to-noise ratio of the individual is rather low. Since there is no obvious reason to expect correlation in the amplitudes and phases of different flow components , further averaging of the signed -coefficients seemed unlikely to improve the statistics of the flow velocity measurement. Instead, estimates of the power (squared velocity) of the flow were made, as in the analysis of Hanasoge et al. (2012). The approach exploits the fact that, by Equation (14), makes a positive, though small, contribution to each and therefore the statistics of the measurements can be improved by summing or averaging the . The noise of the measurements themselves is the main contributor to and must be known to determine the flow power contribution. (Wave realization noise (e.g., Gizon & Birch, 2004) is thought to be the main source of noise in helioseismic measurements.) To estimate the noise power, the -coefficients themselves were averaged over in the same way as the signed coefficients, yielding analogous to . But because the sensitivity of to alternates rapidly with , according to Equation (12), the are expected to have far less sensitivity to the flow velocity than the . Accordingly, was taken to represent the noise contribution to .

For the high- and low-activity epochs the and were averaged over and the squared moduli of the resulting averages were then summed over and averaged over for the observed values of , , and , yielding and . Figure 1 shows and as a function of , for odd and even separately and for the two epochs. The errors in quantities derived from the measured -coefficients are standard deviations based on scatter, where it was assumed that the are statistically independent. For the even- case, the measured power, , displays a striking excess over , the estimated noise power, near solar maximum. The excess is much smaller near solar minimum, consistent with the even-s -coefficients being sensitive to magnetic activity. The analogous excesses for the odd- case, though not as striking as the even- excesses, are statistically significant. Moreover, the approximate equality of the odd- excess power for the high- and low-activity periods does suggest a turbulent, rather than magnetic, origin for the excess power.

The -dependence of the power excesses is of great interest since it should reflect the depth dependence of flow power and solar activity. Summing and over and averaging the results over observed yields and and the -dependent power excesses . The -dependent excesses are not as statistically significant as those produced by combining modes of different . To enhance the statistics of the measurement, the excesses were summed over , at each , yielding . Inspection of Figure 1 suggests that the signal-to-noise ratio of the excess power diminishes for greater than about , so values exceeding were excluded from the sums. Figure 2 shows the -dependence of the summed even- and odd- excesses for the two epochs. The rapid increase in the even- excess power with increasing is consistent with solar activity close to the photosphere (Libbrecht & Woodard, 1990). The much gentler trend in the odd- excess seems to be consistent with large-scale turbulent motions over a substantial range of depths.

The depth-averaged flow velocity can be expressed in terms of using Equation (14). Summing the resulting expression for over and gives

(15) |

for the mean-square depth-averaged flow velocity at degree . The factor appears in the above expression because the quantity , whose -dependence was ignored for simplicity, had to be summed, rather than averaged, over . The factor , rather than , appears in the preceeding equation to ensure that the right hand expression provides an unbiased measurement of .

The measured were averaged over at each , with a weighting inversely proportional to the estimated measurement variance. The averages were multiplied by to provide estimates, , of the power in the depth-averaged velocity up to angular wavenumber . (It was assumed that can be smoothly interpolated to the unobserved even values of and values less than .) Figure 3 shows the measured root-mean-square depth-averaged linear velocity for the two epochs and the average velocity of the different epochs.

## 4 Summary and Outlook

The flow velocities seen in Figure 3 are roughly consistent with large-scale velocity amplitudes observed in the photosphere and with the seismic results of Greer et al. (2015). However, they appear to be at odds with the claims of Hanasoge et al. (2012), although the present measurements do not go as deep as the time-distance measurements. (The lower turning radius of the most deeply penetrating, , modes used for the present analysis, is estimated to be about , corresponding to a depth of about Mm.) The modes sample the flow velocity from the photosphere down to their inner turning radii, with a weighting function that peaks just below the photosphere. Therefore the modest but apparently significant increase in r.m.s. velocity with increasing radial order shown in Figure 3 suggests an even greater increase in velocity with physical depth. Although an inversion for the depth-dependence of the flow velocity was not carried out as part of this analysis, it seems possible that the present results imply a depth dependence consistent with the flow velocities at depths approaching Mm, as shown in Figure 4 of Greer et al.. Note that this figure is not equivalent to Figure 3 of the current paper, because the present analysis is limited to angular wavenumbers while the former includes supergranulation scales. The supergranulation-tracking analysis of Hathaway et al. (2013) also suggests an increase in turbulent power with depth.

In extrapolating the measured flow power, based on near-sectoral (i.e., ) flow components, to all , strong anisotropy-inducing effects of solar rotation (seen in numerical convection simulations) were ignored. Also, the angular wavenumber () dependence of velocity power was ignored in computing the overall power. In addition, only the toroidal flow velocity was measured. For these reasons, the overall level of measured turbulent power could easily be in error by a factor of or more. However, unless the angular anisotropy of large-scale flow patterns changes rapidly between and the photosphere, the increase in the r.m.s. velocity with mode depth would seem to indicate a real increase in the amplitude of large-scale turbulent motion with increasing depth over the sampled depth range, as Greer et al. found.

That the signal-to-noise ratio of the measured flows is significant even at the lowest wavenumbers observed suggests the possibility of detecting turbulent power at angular wavenumbers less than , corresponding to the largest angular scales, and at greater depths than the present analysis permits. To measure larger angular scales a better treatment of spatial leakage is desirable. To probe deeper into the Sun the analysis would need to be extended to oscillation modes of lower . The signal-to-noise of the measurement can be improved, as only a fraction of the existing and projected database of solar oscillations has been analyzed for this work. Similarly, only the near-sectoral-harmonic components of the flow have been utilized in the analysis. The present analysis was based on the dynamical couplings of modes of the same and , which are sensitive mainly to toroidal flow. To detect large-scale poloidal flows the couplings of modes of different and/or would have to be measured.

## Acknowledgements

I thank Curt Cutler and David Hathaway for useful discussions and Tim Larson for help in accessing MDI spherical harmonic time series. I also thank an anonymous referee for constructive comments. This research was supported by NASA grant NNX14AH84G to NWRA/CoRA. This article has been accepted for publication in Monthly Notices of the Royal Astronomical Society, Published by Oxford University Press on behalf of the Royal Astronomical Society

## References

- Böhm-Vitense (1958) Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
- Braun & Fan (1998) Braun D. C., Fan Y., 1998, ApJ, 508, L105
- Chou & Ladenkov (2005) Chou D.-Y., Ladenkov O., 2005, ApJ, 630, 1206
- Christensen-Dalsgaard (2002) Christensen-Dalsgaard J., 2002, Reviews of Modern Physics, 74, 1073
- De Rosa et al. (2002) De Rosa M. L., Gilman P. A., Toomre J., 2002, ApJ, 581, 1356
- Durney (2000) Durney B. R., 2000, ApJ, 528, 486
- Featherstone et al. (2006) Featherstone N. A., Haber D. A., Hindman B. W., Toomre J., 2006, in Proceedings of SOHO 18/GONG 2006/HELAS I, Beyond the spherical Sun.
- Giles (2000) Giles P. M., 2000, PhD thesis, STANFORD UNIVERSITY
- Giles et al. (1997) Giles P. M., Duvall T. L., Scherrer P. H., Bogart R. S., 1997, Nature, 390, 52
- Gizon & Birch (2004) Gizon L., Birch A. C., 2004, ApJ, 614, 472
- Gizon & Birch (2005) Gizon L., Birch A. C., 2005, Living Reviews in Solar Physics, 2, 6
- Gizon et al. (2008) Gizon L., Cally P., Leibacher J., 2008, Sol. Phys., 251, 1
- González Hernández et al. (1999) González Hernández I., Patrón J., Bogart R. S., The SOI Ring Diagram Team 1999, ApJ, 510, L153
- Greer et al. (2015) Greer B. J., Hindman B. W., Featherstone N. A., Toomre J., 2015, ApJ, 803, L17
- Haber et al. (2002) Haber D. A., Hindman B. W., Toomre J., Bogart R. S., Larsen R. M., Hill F., 2002, ApJ, 570, 855
- Hanasoge et al. (2012) Hanasoge S. M., Duvall T. L., Sreenivasan K. R., 2012, Proceedings of the National Academy of Science, 109, 11928
- Hansen et al. (2004) Hansen C. J., Kawaler S. D., Trimble V., 2004, Stellar interiors : physical principles, structure, and evolution, 2nd edn. Springer–Verlag, New York
- Hathaway et al. (2000) Hathaway D. H., Beck J. G., Bogart R. S., Bachmann K. T., Khatri G., Petitto J. M., Han S., Raymond J., 2000, Sol. Phys., 193, 299
- Hathaway et al. (2013) Hathaway D. H., Upton L., Colegrove O., 2013, Science, 342, 1217
- Hathaway et al. (2015) Hathaway D. H., Teil T., Norton A. A., Kitiashvili I., 2015, ApJ, 811, 105
- Howe (2009) Howe R., 2009, Living Reviews in Solar Physics, 6, 1
- Hughes & Thompson (2003) Hughes S. J., Thompson M. J., 2003, in H. Sawaya-Lacoste ed., ESA Special Publication Vol. 517, GONG+ 2002. Local and Global Helioseismology: the Present and Future. pp 307–310
- Jensen et al. (1998) Jensen J. M., Jacobsen B. H., Christensen–Dalsgaard J., 1998, in Korzennik S., ed., ESA Special Publication Vol. 418, Structure and Dynamics of the Interior of the Sun and Sun-like Stars. p. 635
- Kitchatinov & Rüdiger (1999) Kitchatinov L. L., Rüdiger G., 1999, A&A, 344, 911
- Küker & Rüdiger (2005) Küker M., Rüdiger G., 2005, Astronomische Nachrichten, 326, 265
- Larson & Schou (2015) Larson T. P., Schou J., 2015, Sol. Phys., 290, 3221
- Lavely & Ritzwoller (1992) Lavely E. M., Ritzwoller M. H., 1992, Royal Society of London Philosophical Transactions Series A, 339, 431
- Libbrecht & Woodard (1990) Libbrecht K. G., Woodard M. F., 1990, Nature, 345, 779
- Miesch (2005) Miesch M. S., 2005, Living Reviews in Solar Physics, 2, 1
- Miesch (2007) Miesch M. S., 2007, Astronomische Nachrichten, 328, 998
- Miesch et al. (2008) Miesch M. S., Brun A. S., DeRosa M. L., Toomre J., 2008, ApJ, 673, 557
- Scherrer et al. (1995) Scherrer P. H., et al., 1995, Sol. Phys., 162, 129
- Schou & Brown (1994) Schou J., Brown T. M., 1994, A&AS, 107, 541
- Stix (2004) Stix M., 2004, The sun : an introduction, 2nd edn. Springer–Verlag, Berlin
- Thompson et al. (2003) Thompson M. J., Christensen-Dalsgaard J., Miesch M. S., Toomre J., 2003, ARA&A, 41, 599
- Woodard (2006) Woodard M. F., 2006, ApJ, 649, 1140
- Woodard (2014) Woodard M., 2014, Sol. Phys., 289, 1085
- Zhao & Kosovichev (2004) Zhao J., Kosovichev A. G., 2004, ApJ, 603, 776