Frequency shifts of resonant modes of the Sun due to nearsurface convective scattering
Abstract
Measurements of oscillation frequencies of the Sun and stars can provide important independent constraints on their internal structure and dynamics. Seismic models of these oscillations are used to connect structure and rotation of the star to its resonant frequencies, which are then compared with observations, the goal being that of minimizing the difference between the two. Even in the case of the Sun, for which structure models are highly tuned, observed frequencies show systematic deviations from modeled frequencies, a phenomenon referred to as the “surface term”. The dominant source of this systematic effect is thought to be vigorous nearsurface convection, which is not well accounted for in both stellar modeling and modeoscillation physics. Here we bring to bear the method of homogenization, applicable in the asymptotic limit of large wavelengths (in comparison to the correlation scale of convection), to characterize the effect of smallscale surface convection on resonantmode frequencies in the Sun. We show that the full oscillation equations, in the presence of temporally stationary 3D flows, can be reduced to an effective “quietSun” wave equation with altered sound speed, Brünt–Väisäla frequency and Lamb frequency. We derive the modified equation and relations for the appropriate averaging of three dimensional flows and thermal quantities to obtain the properties of this effective medium. Using flows obtained from three dimensional numerical simulations of nearsurface convection, we quantify their effect on solar oscillation frequencies, and find that they are shifted systematically and substantially. We argue therefore that consistent interpretations of resonant frequencies must include modifications to the wave equation that effectively capture the impact of vigorous hydrodynamic convection.
1. Introduction
Measurements of oscillation frequencies of stars provide some of the strongest constraints on stellar structure models. Recent missions Kepler (Borucki et al., 2010) and COROT (Auvergne et al., 2009) have provided us with high quality asteroseismic data, which have revealed structural features of Sunlike stars in unprecedented detail. There have been several successful solar models which have reproduced the oscillation frequencies to within a few parts in thousand. Despite this degree of accuracy, the modeled frequencies display systematic deviations from the observed ones. Understanding the root of these deviations is important both from theoretical and observational point of view — not only is it a window into the missing physics, it is also a hindrance in fitting stellar models. With the ever increasing precision in measured frequencies, it is imperative to understand the physical reason behind the trend and to correct for it while extracting information about the Sun at the same time.
There are various reasons why these shifts are observed, which were categorized as model and modal effects by (Rosenthal et al., 1999) . Some of the possible sources for this bias are:

Errors in physics that go into the modeling, including turbulentpressure contributions (model effect).

Nonadiabaticity prevalent in the outer layers of the Sun (model+modal effect).

Imprecise modeling of propagation of waves through convective flows near the surface of the Sun (modal effect).

Nonlinearity in the wave equation (modal effect).
It is unlikely that any one of these mechanisms would explain the observed deviations completely, rather a combination of these corrections would lead to a better estimate of oscillation frequencies, and possibly eliminate the observed bias. We would like to differentiate between the first and third points — the first assumes a hydrodynamic equilibrium model, but treats oscillations as if they were about an equivalent hydrostatic background, while the third treats the oscillations about a dynamic background. In this article we focus on the third aspect, and propose a method to correct the modal effect. Strong convective flows with speeds close to the sound speed exist near the surface of the Sun. We show that advection of traveling waves by these flows can have a significant impact on the resonant mode frequencies, and it is therefore important to account for this effect in order to extract the correct physics from measurements.
Acoustic modes (p modes) are standing waves bounded from below by a lower turning point set by the Lamb frequency, and from above by an upper turning point set by the acousticcutoff frequency. The upper turning point for modes with frequencies below is significantly below the surface of the Sun, but modes with frequencies lying in the range to encounter an upper turning point close to the surface. These are the modes for which the measured frequencies display a systematic deviation from model, hence it seems likely that the missing physics lies close to the surface, and the deviation is often termed ‘surface effect’. The deviation is predominantly a function of the frequency of a mode, hence corrective terms should also be functions of the mode frequency. An early attempt at mitigating the frequency differences was by Gough (1990), whose prescription was in the form of a power law correction based on the physical process involved: fibrillated magnetic fields lead to a change in sound speed, and hence an observed frequency would be shifted by , where is the mode inertia; while alternate parametrizations of convection lead to a shift . Kjeldsen et al. (2008) had studied the shifts in the Sun, and with an empirical power law model , were able to correct for the shifts in , and . Recently Ball & Gizon (2014) have shown that a combination of the two effects studied by Gough — the cubic shift and the inverse shift — is able to fit the deviations better than Kjeldsen’s power law. These corrections are satisfactory from a modelfitting viewpoint, but by eliminating the frequency differences altogether we lose information about their origin.
Despite the high flow speeds, there still exists a small parameter in the wave equation — the granulation length scales are usually much smaller than horizontal wavelengths of low modes, which are of the order of the solar radius. Hanasoge et al. (2013) had designed an alternate approach to deal with smallscale sound speed perturbation, by expanding the wave equation in terms of the ratio of the two length scales, followed by averaging over small scales. They had shown that if the horizontal wavelength of a propagating wave is much larger than the length scale over which the background medium fluctuates, it is possible to replace the background by a homogeneous, coarsegrained one. In this paper, we use a similar analysis in presence of smallscale flows, and derive an expression for an effective speed. A way to think of homogenization is as follows: nearsurface convection is a stochastic process, evolving both temporally and spatially. Here we assume temporally stationarity, with waves propagating through a frozen configuration of convective flows. Precise measurements of resonant frequencies require taking spectra of long sequences of raw observations of the wavefield — at a minimum, several oscillation periods but typically much longer than that, on the order of a million periods, so as to resolve the mode linewidth and beat down the noise sufficiently. Over this time, waves will have propagated through a large number of realizations of convective flows. Thus one may think of this as an ensemble of wave equations, each ascribed a set of coefficients corresponding to different realizations of the process of convection, drawn from the same probability density function. Each realization of convective flows corresponds to a wave solution, and the question then becomes: how to obtain an ensemble average of these wave solutions? Homogenization provides an answer to this question, specifically in the limit where the horizontal wavelength is much larger than the length scale corresponding to the flows. The principle behind homogenization is that in cases where waves and flows are scaleseparated, measurables such as frequencies can be reproduced by studying wavepropagation through a quiet, homogeneous medium obtained by appropriately averaging over smallscale flows. The method has an advantage that the effective speed which governs wave behaviour on a largescale emerges naturally out of the averaging process. In the asymptotic limit, this is the speed which shall appear in the dispersion relation, and hence shall determine the observed frequency differences. The method of homogenization allows an extension to random media where the distribution function is translationally invariant and bounded (Papanicolaou & Varadhan, 1982), and thus holds promise in progressing beyond our simplified analysis to real solar models.
Several authors such as Murawski & Roberts (1993); Duvall et al. (1998) had studied the effect of random scatterings by convective flows on surface gravity waves (f modes), arriving at corrections to the dispersion relation in order to match the observed frequencies and linewidths. A different treatment to capture the effect of flows on p modes had been carried out by Brown (1984), who had proposed a correction of the form by carrying out a horizontal average over the wave equation. This form of the correction comes closest to our approach, which generalizes the analysis to arbitrary smallscale vector flow velocities. An asymptotic analysis of pmode frequency shifts for small Mach numbers had been carried out by Stix & Zhugzhda (2004), but whether the analysis extends to transonic flows near the solar surface remains unclear.
2. Governing Dynamics
Waves on the Sun are modeled as smallamplitude perturbations about a background, which we assume to be in a steady state. The background medium is characterized by its pressure , density , flow velocity , and gravitational potential . Waves shall induce oscillations about the temporarily stationary background, and the velocity of a moving blob would be given by a sum of the background flow velocity and the mode velocity. We work with the vector flow velocity without splitting it up into horizontal and vertical components, since both the components repeat over small scales horizontally and hence can be treated analogously. In further analysis, we shall suppress explicit spatial and temporal dependence for notational convenience. The full magnetohydrodynamic equation for waves traveling through nonuniform flows was analyzed by Webb et al. (2005), we replicate the equation of motion here. The background flow field satisfies the equation
(1) 
where denotes the secondrank identity tensor defined as . In presence of a flow field , the material derivative of a smallamplitude displacement is defined as
(2) 
We assume that there exists a flow field which satisfies the background equations, and concentrate on the wave equation. The acoustic wave equation, to linear order in the perturbations, is given by
(3) 
where in the exponent denotes the transpose operator (Webb et al., 2005). We shall stick to the Cowling approximation as our focus is on correcting for the wave speed; the analysis can be extended to a general case in an analogous manner.
A propagating wave is scattered out in all direction by the background flow field, which results in a change in the phase of the forwardtraveling wave. Multiple scattering would lead to repeated phase shifts, leading to an altered phase velocity of the wave. This modified velocity would govern the dispersion relation, which would in turn determine the oscillation frequencies. In order to model this scattering perturbatively, we resort to the multiple scale approach of spatial homogenization.
3. Spatial Homogenization
Lowdegree resonant modes on the sun have horizontal wavelengths () comparable to the solar radius, which is much larger than the typical length scales of granulation on the Sun (). In such a case, the dynamics can be parametrized in terms of the scale ratio . In the limit , we can replace the rapidly fluctuating medium by a homogeneous, albeit anisotropic one. To illustrate this approach, we study a special case of a periodic flow in a Cartesian box, where we denote the vertical depth by the coordinate , and label a point on a horizontal plane by the coordinates . The background medium repeats after a length along both the horizontal directions and , and is stratified in the vertical direction. On a particular horizontal plane, the flow takes the form of tessellated square cells. To study the dynamics in the asymptotic limit , we introduce two spatial coordinates which are scaleseparated: a slow coordinate and a horizontal fast coordinate . Note that the vertical component of mode wavelength can be comparable to the scale of flows, so the advection terms might not be scaleseparable vertically. Bearing this in mind, scale separation is being carried out strictly in the horizontal directions, labeled by or equivalently by . The depth coordinate is common to the flows and the wave, and it characterizes the dynamics completely in the vertical direction. In terms of these new coordinates, the spatial gradient can be rewritten as
(4) 
where acts over the coordinates denoted by , and acts over the fast horizontal coordinates denoted by .
We expand the wave field as a perturbative series in the small parameter as
(5)  
The term at order represents the largescale characteristics of the wave field, while the higherorder terms arise as the incoming wave is scattered by the flows. This approach differs from standard perturbative expansions in two ways — firstly the expansion parameter is not in terms of the flow speed but the length scale ; secondly the zerothorder field is not a solution to an equation without flows, but senses an average effect of the flows. Our approach would be to replace the spatially varying medium by an equivalent quiet background, which is different from what would have been there had there been no flows. We substitute Equation (5) in Equation (3) and use Equation (2) to write down the equations order by order in .
At order we get
(6) 
where
(7) 
We introduce an alternate definition of the operator by defining the transpose map explicitly. We introduce the following tensors and tensorial relations:

The doubledot product of two dyadic tensors and defined as , where repeated indices are summed over. This definition can be extended to higher order tensors and , where the last two indices of and the first two of shall be contracted.

A fourth order tensor defined as . The tensor has the property for a dyadic tensor . The definition can be extended to higherorder tensors, where the transpose operator flips the first two indices.

A fourthrank tensor , defined as , that satisfies for any dyadic tensor .
We define the tensor
(8) 
in terms of which we can rewrite the operator as
(9) 
We shall assume henceforth that the gravitational potential is independent of the fast coordinate , which lets us drop any term like or . This assumption is reasonably justified, since density in the nearsurface layers is much smaller than the interior and hence the gravitational potential would depend only weakly on the flows near the surface. Under this assumption, Eq. (6) becomes
(10) 
A particular solution to Equation (10) is a zerothorder field that is independent of the fast coordinate . This is a sufficient, but not a necessary condition. A general solution to Equation (10) can be a function of , which means the operator might be singular. We shall need to invert this operator at order , so it might be necessary to define the inverse precisely by imposing boundary conditions on the fields, or as a generalized inverse, which is beyond the scope of this paper. We shall proceed assuming that is independent of , which fits into our picture of perturbations about a homogeneous background.
Collecting terms of the order , we obtain
(11) 
where
(12) 
and is as defined in Equation (9). We use the fact that is independent of to rewrite Equation (11) as
(13) 
We note two facts  firstly, Equation (13) is linear in , secondly, the right hand side is a product of two terms, one of which varies horizontally over the slow scale while the other varies over the fast scale . This motivates us to seek solutions for by separating the slowly varying horizontal components of from the rapidly varying components of . We substitute
(14) 
where is a thirdrank tensor. Substituting Equation (14) in Equation (11), and using the fact that is independent of , we obtain
(15) 
We seek a solution which satisfies
(16) 
This is conventionally known as a cell problem (see e.g. Bensoussan et al., 1978).
At order we get the equation of motion for the zerothorder field
(17) 
where
(18)  
(19) 
We know that the first order field depends on the zerothorder field through Equation (14); we substitute this in Equation (17) and average over the fast coordinate . The averaging process restricts us to the zerospatialfrequency mode, however we do retain information about smaller length scales through . Given a horizontal cell , we define the cell average of a function as
(20) 
The terms and in Equation (17) are total divergences in , so they drop off on averaging as we impose periodic boundary conditions. We obtain
(21) 
where
(22) 
resembles a wave speed squared term that governs wave propagation on large scales. This is similar in form to the result obtained by Brown (1984), except the term which arises from the separation of scales. We can substitute for from Equation (8) to recast Equation (22) in a more tangible form as
(23)  
Convective cells on the surface of the Sun have a structural form where hot gas rises up near the center and cooler gas sinks down around the edges of the cell. We assume that the flow in a particular cell is antisymmetric about the cell center, hence . This lets us drop the advective term in Equation (21) and obtain
(24) 
Equation (24) is the homogenized wave equation, wherein we have replaced the spatially fluctuating background medium by an effective homogeneous, anisotropic one. The effective speed depends on both the horizontal and vertical flow profiles, the sound speed profile and the density. If density variations in a cell are spatially correlated with the flow speeds, a description in terms of the average speeds and average densities which ignores correlations might not be sufficient.
Periodicity of flows along specific directions leads to a tensorial wave speed, hence we expect anisotropy in wave propagation. Equation (24) is a standard eigenvalue problem which we can solve to find the dispersion relation for vertical order and horizontal wave vector . This dispersion relation captures the extent to which the frequency of a monochromatic wave would appear to differ from the canonical dispersion relation , which is governed only by the adiabatic sound speed. As an example, we consider waves propagating in a 2dimensional uniformly dense fluid. We assume that the fluid film is tiled with smallscale, periodic flows . A plane wave propagating along the axis would have the form
(25) 
where is the amplitude, is the component of the wavevector and is the frequency of the wave. Neglecting changes in gravitational field in the plane, we can rewrite Equation (24) in Fourier space as
(26) 
where the terms referred to as to arise from the solution to the cell problem. The corresponding dispersion relation can be obtained by diagonalizing the matrix in the second term of Equation (26). If the terms to are negligible, the dispersion relation would have the form , but this might not be a good approximation in general.
We note that Equation (24) is not damped at the leading order, so attenuation, if any, would be captured at a higher order. This is consistent with observations, since the observed pmode linewidths are much smaller than the corresponding frequencies.
4. Numerical test of flow homogenization
In this section we take a simple example of waves propagating through flows in a 2D plane, and show that we can effectively reproduce the wave speed by homogenizing over small scale flows. A thorough analysis of the approach would include realistic pressure, density and soundspeed profiles, possibly starting from simulations of nearsurface structure of stars. In this work, we focus only on soundspeed variations induced by flows, and assume that the background pressure and gravitational field are unchanged in presence of the flow. We also assume that the background medium is uniformly dense throughout. Under these assumptions, Equation (3) simplifies considerably as we can drop the terms and . We combine Equations (2) and (3) to obtain
(27) 
We assume that the waves are propagating in a medium that is periodic along both the directions; in other words, the unit box is tessellated to fill up the entire twodimensional plane. We choose to work in Cartesian coordinates , with the origin at the center of a cell. As an example of a flowfield in a cell, we choose a vortex centered about the cellcenter. The functional form of the flow is given by
(28) 
where is a velocity scale, is a length scale of the flow related to the cellsize . The sound speed would depend on several factors such as the pressure in a cell and the equation of state. For our illustrative purpose, we choose a simple model given by
(29) 
where is the sound speed at the center of the cell and is the sound speed far away from the center of the cell. We choose a cell lengthscale , and values of sound speed and flow parameters as listed in Table 1, where refers to the magnitude of the flow velocity .
Run  
Run  
Run  
Run 
We recast Equation (27) in the form of an eigenvalue problem by using Equation (2) along with Equation (27). Assuming harmonic time dependence, we rewrite the system of equations as
(30) 
where the operator encapsulates the restoring force on the wave field and is given by . We note that both and are periodic with the same length scale as the flow, hence we use Bloch’s theorem to reduce the problem to one cell. According to Bloch’s theorem, the wavefield can be written as , where is a horizontal, as yet arbitrary, wavevector, and the function has the periodicity of the flow. We study waves propagating along the xaxis which have . We impose a largescale periodicity condition , which restricts us to discrete values of the wavevector for integral . Corresponding to each value of we obtain an eigenvalue problem from Equation (30), which we solve numerically for the frequency . We refer to this relation as ‘true dispersion relation’. In the limit of infinite wavelength, we can replace the fluctuating medium by a homogeneous one, given by Equation (24). We make one further approximation — because of potential nonuniqueness in computing , we restrict ourselves to the zerothorder homogenized speed given by
(31) 
The homogenized equivalent to Equation (24) in this case is
(32) 
where is as defined in Equation (31). We compute the dispersion relation corresponding to Equation (32), and refer to it as ‘homogenized relation’. It is worth pointing out that the approximation made in Equation (31) is crude; one can not neglect if Mach numbers are close to , and it might be necessary to compute by solving Equation (16) to achieve a higher level of accuracy.
We compare the zerothorder homogenized dispersion relation with the true dispersion relation in Figure 1. We also compare it with dispersion relations of the form and , which do not explicitly account for flows. We have listed frequency differences for several different parameter values in Table 2. We find that the homogenized relation replicates the true dispersion relation better than the alternate relations. The relative frequency differences corresponding to Run 2 have been plotted in Figure 2. Comparing with Run 1, we note that the frequencies predicted by the zerothorder approximation to the homogenized dispersion reduce in accuracy as the Mach number increases.
(Run )  (Run )  (Run )  (Run )  

5. Impact on solar frequencies
Correcting the solar oscillation equations involves correcting for both modal and model effects, as described by Rosenthal et al. (1999). In a recent paper, Piau et al. (2014) have explored the impact of various parametrization of stellar convection — from mixinglength theory (Spruit, 1977) to patched hydrodynamic models — on oscillation frequencies. They have used horizontallyaveraged sound speeds from these models in the wave equation. We have shown that the presence of flows introduces further advective corrections in the wave speed, and the horizontallyaveraged sound speeds do not represent the true wave speed. The actual speed of the waves is obtained from Equation (24), and the difference between this speed and the horizontallyaveraged sound speed can be substantial near the solar surface.
In this section, we start from Model S (ChristensenDalsgaard et al., 1996) which uses the mixinglength formulation of convective flux. We correct wave speeds using flow profiles obtained from simulations by Schüssler (personal communication) using the MURaM code (Vögler et al., 2005). The flowspeed profiles have been plotted in Figure 3. We start from equation (24) and work with lowdegree modes , which are predominantly radial. We assume that the Mach number is low enough for the zerothorder homogenization approximation (Equation 31) to hold, and retain only the dominant component of which goes roughly as . The exact form of which goes into the correction depends on the direction of advection at every depth; for our purpose we shall restrict ourselves to the special cases of radial flow and tangential flow . We refer to this modified model as Model Hom, and compute the eigenfrequencies corresponding to this model assuming the same boundary conditions as Model S. We compare these frequencies with observed solar frequencies obtained by the Birmingham Solar Oscillation Network (BiSON) (Chaplin et al., 2007). We find that the presence of flow advection in the wave equation leads to a reduction in mode frequencies, which reduces the observed difference between observed and Model S frequencies.
To study the impact of modal versus model corrections on frequencies, we compare the Model Hom frequencies with another hydrostatic model (Basu & Antia, 1994, henceforth called Model BA), which uses the CanutoMazzitelli formulation of convective flux (Canuto & Mazzitelli, 1991). The frequency differences between these models and the observations by BiSON have been plotted in Figure 4. We find that the modal frequency shifts due to the presence of advective flows is comparable to those obtained from modeling differences. Both the effects need to be looked into while fitting solar frequencies.
Some of the assumptions made in this section, such as diagonallydominated and isotropic , are expected to be violated considerably near the surface, but the present example demonstrates that the effect of flow correction can be substantial. A detailed computation, starting either from a fully hydrodynamic background or a patched hydrodynamic background, is necessary to make quantitative statements about frequency shifts due to convective flows.
6. Conclusion
Observed solar oscillation frequencies can differ from modeled frequencies because of an inaccurate background model, as well as incomplete description of mode physics. The presence of nearsurface convective flows on the Sun results in mode frequencies which are different from ones computed starting from a quiet background. Threedimensional convection simulations generally result in an increased size of the acoustic cavity, thereby reducing the frequencies (Rosenthal et al., 1999; Piau et al., 2014). In this paper, we have shown that the impact of advection on waves might result in frequency shifts of similar magnitudes. In presence of flows, the appropriate wave speed is not the horizontallyaveraged sound speed, rather it is a combination of sound speed and an appropriate projection of the flow speed in the direction of propagation.
Granulation length scales on the solar surface are much smaller than horizontal wavelengths corresponding to low modes, which lets us replace the irregular background medium by an effective, homogeneous one. Acoustic waves propagating through steady background flows experience a change in phase speed on being scattered by the flows. The effective wave equation in the largewavelength approximation has the form
where the effective speed is a tensor, given by Equation (22).
The impact of time evolution of flows, missing in the present analysis, needs to be studied in detail. The spatial scales of waves and background flows are significantly different, thereby allowing a two scale analysis. The timescales of granule evolution are, however, not very different from time periods of the wave. Overlap of scales generally results in strong scattering and forcing, thus a picture of flows evolving in time might be necessary to capture frequency shifts and damping. Added to this are the impact of superadiabaticity and magnetic fields, which have been ignored in the present analysis, but might be important in governing the spectrum. Nevertheless the present analysis provides a useful starting point to model wave propagation through convective flows and resultant frequency shifts.
7. Acknowledgment
We thank the anonymous referee for his useful suggestions which improved the quality of this work significantly. JB would like to acknowledge the financial support provided by the Department of Atomic Energy, India.
References
 Auvergne et al. (2009) Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411
 Ball & Gizon (2014) Ball, W. H., & Gizon, L. 2014, ArXiv eprints, arXiv:1408.0986
 Basu & Antia (1994) Basu, S., & Antia, H. M. 1994, Journal of Astrophysics and Astronomy, 15, 143
 Bensoussan et al. (1978) Bensoussan, A., Lions, J., & Papanicolaou, G. 1978, Asymptotic analysis for periodic structures (NorthHolland Publishing Company)
 Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
 Brown (1984) Brown, T. M. 1984, Science, 226, 687
 Canuto & Mazzitelli (1991) Canuto, V. M., & Mazzitelli, I. 1991, ApJ, 370, 295
 Chaplin et al. (2007) Chaplin, W. J., Elsworth, Y., Miller, B. A., Verner, G. A., & New, R. 2007, ApJ, 659, 1749
 ChristensenDalsgaard et al. (1996) ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286
 Duvall et al. (1998) Duvall, Jr., T. L., Kosovichev, A. G., & Murawski, K. 1998, ApJ, 505, L55
 Gough (1990) Gough, D. O. 1990, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 367, Progress of Seismology of the Sun and Stars, ed. Y. Osaki & H. Shibahashi, 283
 Hanasoge et al. (2013) Hanasoge, S. M., Gizon, L., & Bal, G. 2013, ApJ, 773, 101
 Kjeldsen et al. (2008) Kjeldsen, H., Bedding, T. R., & ChristensenDalsgaard, J. 2008, ApJ, 683, L175
 Murawski & Roberts (1993) Murawski, K., & Roberts, B. 1993, A&A, 272, 595
 Papanicolaou & Varadhan (1982) Papanicolaou, G. C., & Varadhan, S. R. 1982, Statistics and probability: essays in honor of C. R. Rao, 547
 Piau et al. (2014) Piau, L., Collet, R., Stein, R. F., et al. 2014, MNRAS, 437, 164
 Rosenthal et al. (1999) Rosenthal, C. S., ChristensenDalsgaard, J., Nordlund, Å., Stein, R. F., & Trampedach, R. 1999, A&A, 351, 689
 Spruit (1977) Spruit, H. C. 1977, PhD thesis, Thesis University of Utrecht, The Netherlands.
 Stix & Zhugzhda (2004) Stix, M., & Zhugzhda, Y. D. 2004, A&A, 418, 305
 Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
 Webb et al. (2005) Webb, G. M., Zank, G. P., Kagashvili, E. K., & Ratkiewicz, R. E. 2005, Journal of Plasma Physics, 71, 785