Abstract
Convection in the Sun’s outer envelope generates turbulence and drives differential rotation, meridional circulation, and the global magnetic cycle. We develop a greater understanding of these processes by contrasting observations with simulations of global convection. These comparisons also enhance our comprehension of the physics of distant Sun-like stars. Here, we infer toroidal flow power as a function of wave number, frequency, and depth in the solar interior through helioseismic analyses of space-based observations. The inferred flows grow with spatial wave number and temporal frequency and are confined to low latitudes, supporting the argument that rotation induces systematic differences between the poles and equator. In contrast, the simulations used here show the opposite trends—power diminishing with increasing wave number and frequency while flow amplitudes become weakest at low latitudes. These differences highlight gaps in our understanding of solar convection and point to challenges ahead.
INTRODUCTION
Despite several decades of space- and ground-based high-resolution observations and substantial theoretical and computational effort, fundamental questions on solar dynamics, such as the emergence of differential rotation, meridional circulation, and solar magnetism, remain unanswered. Convection is thought to be centrally important to these phenomena but is among the least understood aspects of solar dynamics due to the extreme parameter ranges in which it manifests. Gaps in our theoretical appreciation of convection persist. The mixing-length theory (1), a widely adopted prescription for modeling stellar convection and evolution (2), is a coarse description and does not predict flow details. Numerical simulations and laboratory experiments are also unable to access the relevant parameter regime (3, 4). Under the circumstances, a fruitful avenue to pursue consists of contrasting specific details between simulations and observations, highlighting disagreements and resolving them where possible. Although comparisons are stymied because it is often the case that simulations and observations do not pertain to same quantities, they represent important steps toward an improved understanding. The first such attempt using helioseismology concluded that convective velocity amplitudes were small in comparison to numerical simulations (5). This result was subsequently challenged (6) by an analysis using a different helioseismic technique, which inferred amplitudes to be consistent with numerical simulations. The disagreement remains unresolved, as a comparison between a multitude of constraints, helioseismic or otherwise, and different numerical simulations (7) reveals. A major source of this conflict is the lack of consistency in various parameters in the comparisons: differences in the extents of temporal and spatial averaging, the type of flow (poloidal or toroidal), etc. These efforts, however, were productive, focusing attention on this important question and have moved the field forward, leading to the development of concrete models and new results on the physics of convection (8).
Solar convection is a multiscale phenomenon, occurring over broad ranges in spatial and temporal frequencies. The large scales are very interesting since they likely traverse layers deep in the convection zone (3, 7) and are therefore very challenging to model, given that they encounter vast gradients in density and adiabatic index. A major puzzle concerns flows on the largest scales, which are observed at the surface to contain the least power. In contrast, classical turbulence theory suggests that most of the turbulent power lies in spatial scales comparable to the system size. By this definition, the Sun is not a classical turbulent system, since velocity power smoothly decays with decreasing wave number. Supergranulation, another puzzling feature of convection, spans the frequency range from 3 to 11 μHz and spherical harmonic degrees s dominantly ranging from 80 to 140, with a power that goes down to very low wave numbers as well. Therefore, the very largest scales of convection are likely to evolve at even lower frequencies, i.e., 3 μHz, and wave numbers of 100. On the basis of these estimates, in the present analysis, we infer toroidal flows within the temporal frequency range of [0,3.05] μHz, for wave numbers of [1,50].
MATERIALS AND METHODS
Helioseismology provides a window into the internal dynamics of the Sun. Vigorous convective turbulence at the surface generates small-amplitude acoustic waves, which propagate in the solar interior and re-emerge at the surface, forming resonances in cavities (9). The modulation of the solar surface by acoustic waves, measured as Doppler shifts in absorption lines, is subjected to seismic analysis. In contrast to other helioseismic techniques, measuring normal-mode coupling from observations is straightforward (10–13). The additional benefit is that interpretation is tightly connected to well-posed semianalytical models; the price is that these models are nontrivial to evaluate and rely on tedious algebra. Normal modes also provide information about the poloidal-toroidal structure of fluid motions.
The method involves first constructing eigenfunctions for a reference model of the Sun, e.g., model S (9), which, for the linear solar wave operator, form an orthonormal and complete basis (14). Since the Sun contains convective flows and the reference model does not, the “true” eigenfunctions depart from the reference set. Specifically, convection is treated perturbatively since velocity amplitudes on large scales are substantially smaller than the local sound speed (3). Completeness of the reference basis allows us to express the true eigenfunctions as a weighted linear superposition of reference eigenfunctions. Modes of the Sun are thus in a “coupled” state in relation to the reference. The weights, termed as coupling coefficients, may be measured from observations and related to model perturbations through a straightforward inverse problem. The improvement in inferential accuracy obtained by including differential rotation in the background model, which is known to cause changes to the reference basis (10), is at most on the order of a few percent for the modes studied here. However, since the computational and technical overhead required to expand around the rotationally distorted eigenfunctions is significantly greater, we expand around the degenerate eigenfunctions of model S.
A useful aspect of this methodology, compared with other techniques of helioseismology, is that it allows us access to time variations of nonaxisymmetric perturbations on a global scale, an appropriate formalism for the application to dynamically evolving phenomena such as convection and Rossby waves. In addition, it provides a rigorous formalism with which to analyze vector flows, enabling a proper comparison with simulations, a notable improvement over prior work. In particular, we apply the Chandrasekhar-Kendall relation, with which vector flows may be decomposed into poloidal (first two terms of Eq. 1 below) and toroidal quantities (last term) (15)
A snapshot of computed convection (A) and simulated and observed differential rotation (B and C) (38). Although Taylor-Proudman columns appear at low latitudes in the simulations, the amplitudes and contours compare favorably. (D) Numerical simulations of convection (4) show that toroidal
Data analysis and inversion
The Helioseismic and Magnetic Imager (HMI) (18), onboard the Solar Dynamics Observatory, has been observing the Sun since 2010. We analyze 8 years of HMI global-mode data in the ℓ range of 50 to 192 (ℓ is harmonic degree) and select modes with resonant frequencies in the range of 1800 μHz ≤ ωnℓ/(2π) ≤ 3600 μHz, where n denotes the radial order (19). The upper frequency limit of 3600 μHz is placed to exclude modes with known systematics (19). From the lower end, we exclude modes with very small linewidths (that low-frequency modes have), since characterizing the frequency dependence of convection requires sampling modes off resonance. Small linewidths imply that measuring the mode even slightly off resonance can lead to background power seeping into inferences. The global time series data are further divided into eight segments of 360-day-long sets.
These time series are transformed to the temporal Fourier domain and denoted by
Correlations thus computed contain several independent parameters ω, σ, t, m, and ℓ and are also noisy, rendering interpretation difficult. Consequently, we introduce B coefficients, which condense and average the information present in raw correlations. The coefficient
The function
Modes in the Sun are stochastically driven, and the excitation process is assumed to be statistically independent for each wave number and temporal frequency channel because, at any given time, a very large number of isotropically distributed (roughly) granules across the surface of the Sun are each exciting waves with random phases and amplitudes. This implies that, for reference modes,
The frequencies and other mode parameter fits vary from one 360-day period to another, which is why the noise background power shows dispersion. The theoretically computed noise background follows the observed B-coefficient power closely except at low s − ∣t∣, where the observed B-coefficient power is significantly larger than the model. This additionally gives us confidence that our models for the noise are reasonable.
Last, we comment on the trustworthiness of the result. We have experimented with different sets of modes; specifically, we studied the results using the entire set of modes, i.e., those that have been successfully identified by the HMI pipeline from analyses of 360-day datasets (26). We found that the inferred amplitudes of prograde and retrograde flows are somewhat sensitive to the high-frequency modes (>3600 μHz), which are known to be problematic (19). However, one feature that is robust to all choices of mode sets is the variation of toroidal flow power with s − ∣t∣, i.e., sectoral modes are dominant and there is very little power in tesseral and zonal modes (for the spatiotemporal band under consideration here).
Simulation setup and importance of toroidal flows
At the outset, it is important to note that this simulation is merely one of innumerable possibilities, with many possible choices for turbulence modeling, boundary conditions, resolutions, numerical schemes, etc. The comparisons with simulations here are therefore solely a commentary on the setups chosen for this work.
There is simply no means by which simulations will be able to reach solar convection parameter regimes in the foreseeable future. However, it is within the realm of possibility that there exists an effective theory of convection, i.e., where all the ingredients required to properly recover large-scale properties of these flows are captured in the calculations. An instance of this is the simulation of surface granulation, which is able to reproduce line formation, granule sizes, shapes, correlation functions, etc. [e.g., (27)] despite the very large range of scales prevalent in near-surface turbulence. Thus, the question is as follows: Are simulations of global convection at a stage where we can trust them to accurately reproduce large-scale flows? And if not, then can we place observational constraints to eliminate or differentiate one calculation from another?
In this effort, we have tried two other numerical setups (not reported here)—the first being a nonrotating high-resolution simulation (1536 × 3072 × 512), with the upper boundary at r/R⊙ = 0.99 and another with the same boundary and resolution, but with rotation included. Neither of these appeared satisfactory to us because the first had no rotation and was therefore isotropic across the 2-sphere; the second showed antisolar rotation (i.e., with poles spun up relative to the equator), and we therefore discarded both. The reason is that high-resolution calculations tend to have high Rossby numbers, lower viscosities and thermal diffusivities; the latter two effects introduce large (relatively speaking) convection velocities. Further, to obtain solar-like (fast equator) differential rotation, the Rossby number must be kept as low as is possible (28). Consequently, we settled on a low-resolution case with the upper boundary at r/R⊙ = 0.94 and solar-like rotation, i.e., equatorial spin up. To reiterate, this is just a small subset of the range of simulations that are possible with a variety of different codes.
The numerical simulation we use here (29) resolves rotating (Fig. 1, A and B) convection using 384 × 192 × 64 points, along the longitude, latitude, and radius, respectively, with r/R⊙ ≈ 0.96 serving as the upper boundary of the computational domain. The hydrodynamic equations in three-dimensional spherical geometry (with the entropy equation instead of the total energy equation) are solved using a fourth-order central derivative scheme and stepped forward in time with the Runge-Kutta method. We also apply a slope-limited artificial diffusion to all variables to stabilize the calculation. Explicit diffusion is applied for velocity and entropy, whose values at the top boundary are 1 × 1013 cm2 s−1 and 2 × 1013 cm2 s−1, respectively, and proportional to ρ−1/2 (where ρ is the density), to reproduce a rotation profile with a relatively fast rotating equator in comparison to the pole. In comparisons presented here, the static t = 0, σ = 0 component of the numerical data is removed.
RESULTS
Figure 3 (A and B) shows that, in contrast to simulations where power is largest at low spatial wave numbers and temporal frequencies, observationally inferred velocity monotonically decreases with decreasing wave number and frequency. This puzzling trend, in contrast to classical turbulence, where small wave numbers are energy containing, has been the focus, albeit inconclusive, of a variety of theoretical and computational efforts (4, 30–32). While the disparity in amplitude is not as significant as suggested by prior time-distance methods (5), the two behaviors are markedly different. To benchmark the seismic results against independent measurements of the surface toroidal-flow spectrum obtained from Doppler data, we compare the inferred power spectrum at r/R⊙ = 0.995 with surface-Doppler estimates (33) in Fig. 3C. Specifically, we use figure 12 of (33), multiplying it by
Shaded areas in (A to C) mark 1-σ errors in velocity. The legend in (A) applies to (B) as well. In contrast to the simulation, where velocity power peaks at s ≈ 15 and continuously decreases in either end, we see in (A to C) that power weakens with decreasing wave number (s), indicating that large-scale toroidal flows are suppressed. The comparison in (C) between toroidal flow amplitudes measured directly from surface Doppler measurements (33) and seismic inferences at r/R⊙ = 0.995 is favorable; the simulated power does not make an appearance since the upper boundary of the computational domain is r/R⊙ = 0.94. (D to E) Inferred RMS velocity compares reasonably with simulations [in contrast to prior estimates (5)]. Horizontal bars denote the width of the averaging kernel in (D to E). (F) RMS ratios of simulated-to-observed velocity as functions of harmonic degree s at different depths. Vertical error bars are computed using the variance in the spectra obtained from each of the eight 360-day periods analyzed here.
Spectra as a function of s − ∣t∣, which capture the geometry of convection, show a dominance of sectoral modes in the Sun, i.e., small s − ∣t∣ (Fig. 4, A and B), with prograde (t > 0) and retrograde (t < 0) components indistinguishable to within error bars (for the tracking frequency of 453 nHz). In contrast, simulations contain a much broader set of convective modes. This is more markedly seen in Fig. 4 (C and D), which highlight the differences between the latitudinal profiles of simulated and observed flow, where we plot power as a function latitude for simulations (Fig. 4C) and seismology (Fig. 4D). The error bars on the observed curves are invisible on this plot, whereas the simulated curve is very noisy: This occurs because the low spatial and temporal wave numbers are preferentially powered in the simulations, thereby requiring a much longer integration time to achieve convergence. It is interesting that the most vigorous toroidal flows are cospatial with the fastest rotating layers, i.e., differential rotation and toroidal power peak at the equator. Whether anisotropic toroidal flow power induces cospatial rotational shear, or the other way around, remains to be investigated. The similarity with recent measurements of Rossby waves in purely sectoral channels [e.g. (21, 22)] tempts one to speculate if there are more general underlying processes at play. The absence of nonsectoral modes in toroidal flow may be due to multiple reasons: low-amplitude excitation, exclusion of tesseral/zonal modes due to specific boundary conditions, or rapid damping (possibly by rotational shear). The most vigorous toroidal flows appear to lie within the spherical wedge delineated by the intersection of the cylinder tangent to the radiative interior with the sphere, especially so in the near-surface layers.
The legend in (A) applies to (B) as well. Low s − ∣t∣ correspond to sectoral harmonics, i.e., equatorially localized north-south–aligned features, whereas harmonics with large s − ∣t∣ are present at all latitudes. Having obtained flow coefficients
Figure 5, a reconstruction of flow images from observations, shows that equatorially localized north-south–aligned flows are preferentially powered—these are the high-amplitude sectoral modes, i.e., for s − ∣t∣ ≲ 12 in Fig. 4 (A and B). Simulated velocity images qualitatively differ from the inferences, with large-scale features distributed across all latitudes (Fig. 5, C and D). In contrast, the small scales appear to be preferred in the Sun (Fig. 5, A and B), due to the inferred power spectrum increasing continuously as a function of wave number. This spectral dependence results in the highest spatial frequencies having the most power, whereas the simulated spectra peak at low spatial frequencies and decrease with increasing wave number.
Instantaneous snapshot of toroidal flow, reconstructed from the inferred (A and B) and simulated (C and D)
In Fig. 6 (A and B), we show the dependence of toroidal flow power on temporal frequency—the observed velocity scale increases continuously with frequency, broadly in line with the idea that convection at higher frequencies is more vigorous than at lower frequencies, in contrast with computation (Fig. 6C). Figures 3 (A to C) and 6 (A and B) together indicate that solar flow power is continuously increasing as a function of wave number and frequency (over the ranges considered here), eventually overtaking and exceeding the simulations. Prior comparisons [e.g., (5)], which relied solely on amplitudes, thus miss the point that velocity-amplitude comparisons are sensitive to the spatial and temporal frequency bands under consideration and that, even if the comparisons were favorable, simulations and inferences could still be capturing intrinsically different physics (e.g., classical versus nonclassical turbulent systems, as discussed early on in this section). In the present work, we highlight the importance of making detailed comparisons with observations to improve faith in theoretical and computational models.
The spectrum
Figure 6D tracks “deep” and shallow convective velocities, i.e., the inferences at radii r/R⊙ = 0.90,0.99, as a function of the magnetic cycle. It is conceivable that these turbulent flows may be weakened during periods of high solar activity because magnetism suppresses free movement of plasma (34). However, Fig. 6D shows that there is no apparent correlation between surface and interior velocities with the sunspot cycle (35).
DISCUSSION
In summary, the state-of-the-art numerical simulation used here (29) reproduces certain features such as overall amplitudes—to within a factor of a few—of the toroidal flow but are not good models of solar turbulence, as evidenced by their inherent difference from seismic observations: (i) Large-scale flows are preferentially powered in simulations, whereas they are suppressed in the Sun, and (ii) simulated flows are weak at the equator and vigorous at high latitudes, in contrast to inferences, which show that they are confined to the equatorial regions and weak at high latitudes.
Differential rotation is thought to be sustained by a combination of a thermal wind and Reynolds stresses (3). A thermal wind could be induced by a latitudinal entropy gradient, in turn caused by differences in convection between the poles and equator (36), resulting in changes to the turbulence. Our inferences here of latitudinally nonuniform toroidal flow appear to match these symptoms, although they do not directly point to the cause. It would be very useful to theoretically study the impact of thermal winds on toroidal flows in future investigations. Inferences of Reynolds stresses would go a substantial distance toward discerning the origin of differential rotation, both in the Sun and stars (37), which state-of-the-art simulations are not successful at identifying. However, these are much more challenging to measure since they require the complete set of phases and amplitudes of toroidal and poloidal flows at all spatial and temporal scales (see the Supplementary Materials). Since inferences of toroidal flows here are obtained from studies of same-ℓ couplings, issues arising from leakage are significantly mitigated. Modes for ℓ < 150 or so are typically well separated in frequency (on the order of several to tens or more of microhertz), and attempts to infer perturbations varying on temporal frequencies much smaller than the nominal frequency separation between pairs of modes, as used in the analysis here, are hindered by the influence of spatial leakage and background power on the measurement. As a consequence, cross-ℓ couplings must be interpreted much more carefully and are deferred to a future study. The penalty for considering only same-ℓ couplings is that we are only able to infer odd harmonic degree toroidal flows and have practically no sensitivity to poloidal components.
Convection comprises spatially complex and temporally evolving structures; testing different theories of this process requires building a framework with which to make consistent comparisons. The theory and computation of convection are sensitive to a large variety of necessary inputs, e.g., parameters and boundary conditions, and eventual solutions may vary substantially. Determining whether a simulation is “correct” is for instance a difficult question—would the successful recovery of differential rotation and meridional circulation be sufficient cause to declare success or should we seek a greater degree of agreement, i.e., turbulence structure functions, Reynolds stresses, etc.? We have shown, for example, that convective amplitudes alone, while providing some broad indications of the success of a theory, are insufficient for making accurate and meaningful comparisons. To eliminate or validate a theory may thus demand agreement with both large-scale circulations (e.g., differential rotation) and other constraints, such as provided here. Thus, the present disagreement should be viewed as a significant step toward building a framework with which to make comparisons of theory and computation with observation and adding to the range of available measurements of convection. One lesson from the direct comparison with simulations here (29) is that high resolution alone is insufficient to achieve agreement and that a great distance still remains to be traversed between computation and observation.
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/30/eaba9639/DC1
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
REFERENCES AND NOTES
- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).