Using eigenmode-mixing to measure or constrain the Sun's interior B-field

TL;DR

This paper proposes a novel method using eigenmode-mixing of solar oscillations to measure or constrain the Sun's interior magnetic field, addressing a longstanding challenge in solar physics.

Contribution

It introduces a new detection statistic for non-axisymmetric mode-mixing in solar data and analyzes its effectiveness for measuring the Sun's interior magnetic field.

Findings

Detection statistic is robust against spatial aliasing.

Signal-to-noise ratio increases with the number of mode pairs.

B-induced mode-mixing should be detectable for various B-field strengths.

Abstract

Understanding the generation and distribution of the Sun's interior magnetic (B-) field is a longstanding challenge. Here we describe how measurements of the Sun's oscillation eigenfunctions might be used to measure the Sun's interior B-field. The B-field induces mode-mode couplings, causing the angular patterns of the eigenfunctions to differ from simple Y_{lm}'s We concentrate on the magnetic coupling between modes with the same (n,l) values and different but nearby $m$-values, since these non-axisymmetric couplings clearly cannot be due to the Sun's axisymmetric rotation and since for these cases, mode mixing is enhanced by the near-degeneracy of the mode frequencies. We analyze magnetically-induced mode mixing in two stages of increasing complexity: first neglecting mode damping, and then incorporating realistic damping rates. We introduce a novel detection statistic that tests for the presence of non-axisymmetric mode-mixing in Solar Doppler data. We show that our detection statistic is naturally robust against spatial aliasing. We estimate our statistic's signal-to-noise ratio (SNR) as a function of the mode-mixing amplitude. While B-induced mode-mixing is probably not detectable in a single mode pair, we argue that the phase of the B-induced mixing should be approximately the same across a wide range of modes. The total SNR then grows roughly as Np^{1/2}, where Np is the number of mode pairs. We conclude that B-induced mode-mixing should be detectable for a fairly wide range of B-field magnitudes and geometries.

Equations

$$\vec{v}(\theta ,\varphi ,t)=\sum _{nlm}i{\omega }_{nlm}{A}_{nlm}{\vec{\xi }}_{nlm}(R,\theta ,\varphi )\exp [i{\omega }_{nlm}t],$$

$$\displaystyle\left(\begin{array}[]{cc}1+\delta&\ \ 0\\ 0&\ \ 1\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{cc}1+\delta&\ \ \epsilon\\ \epsilon^{*}&\ \ 1\end{array}\right)$$

$$\rho {\ddot{\xi }}^i=(non\text{–}Bterms)+{\gamma }_j^i{\xi }^j+{\kappa }_{jk}^i{\nabla }^j{\xi }^k+{\eta }_{jkl}^i{\nabla }^j{\nabla }^k{\xi }^l$$

$$\ddot{a}+{\Gamma }_a\dot{a}+{\omega }_a^{2}a$$

$$\ddot{b}+{\Gamma }_b\dot{b}+{\omega }_b^{2}b$$

$$\tilde{a}(\omega )$$

$$\tilde{b}(\omega )$$

$$G_a(\omega )$$

$$G_b(\omega )$$

$$\tilde{f}(\omega )\equiv (2\pi {)}^{-1/2}{\int }_{-\infty }^{\infty }f(t)e^{-i\omega t}dt.$$

$$\int a^{\ast }(t)b(t)e^{-i\Omega t}dt,$$

$$\int {\tilde{a}}^{\ast }(\omega )\tilde{b}(\omega +\Omega )d\omega .$$

$$\int W(\omega ){\tilde{a}}^{\ast }(\omega )\tilde{b}(\omega +\Omega )d\omega .$$

$$\displaystyle\int G^{*}_{a}(\omega)G_{b}(\omega+\Omega)\bigg{[}\tilde{n}^{*}_{a}(\omega)\tilde{n}_{b}(\omega+\Omega)$$

$$\displaystyle+\lambda^{*}\bigg{(}|\tilde{n}_{b}(\omega+\Omega)|^{2}G^{*}_{b}(\omega+\Omega)+|\tilde{n}_{a}(\omega)|^{2}G_{a}(\omega)\bigg{)}\bigg{]}d\omega\,.$$

$$\hat{\lambda}^{*}\omega_{0}^{2}|\tilde{n}_{0}|^{2}\int{G^{*}_{a}(\omega)G_{b}(\omega+\Omega)\bigg{(}G_{a}(\omega)\,+\,G^{*}_{b}(\omega+\Omega)\bigg{)}\,d\omega}\,.$$

$$\frac{-4{\omega }_a\hat{\omega }+i\Omega \Gamma }{{\omega }_a^2(4{\hat{\omega }}^2+{\Gamma }^2)}$$

$$\bigg{[}\omega_{a}^{2}\,\bigg{(}\Omega^{2}+\Gamma^{2}-4\Omega\hat{\omega}-2i\,\Omega\,\Gamma\bigg{)}\bigg{]}^{-1}$$

$$\displaystyle S=\int{W(\omega)\frac{\big{(}-4\hat{\lambda}^{*}\omega_{0}^{2}\omega_{a}^{-3}|\tilde{n}_{0}|^{2}\big{)}\big{(}\hat{\omega}-\frac{i}{4}\Omega\Gamma\omega_{a}^{-1}\big{)}}{\big{(}4\hat{\omega}^{2}+\Gamma^{2}\big{)}\big{(}\Omega^{2}+\Gamma^{2}-4\Omega\hat{\omega}-2i\Omega\,\Gamma\big{)}}d\omega}\,.$$

$$\big{(}\Omega^{2}+\Gamma^{2}-4\Omega\hat{\omega}-2i\,\Omega\,\Gamma\big{)}\,,$$

$$\displaystyle S=\int{W(\omega)\frac{\big{(}-4\hat{\lambda}^{*}\omega_{0}^{-1}|\tilde{n}_{0}|^{2}\big{)}\hat{\omega}}{\big{(}4\hat{\omega}^{2}+\Gamma^{2}\big{)}\big{(}\Omega^{2}+\Gamma^{2}\big{)}}d\omega}\,.$$

$$N=\sum _{i}W({\omega }_i)G_a^{\ast }({\omega }_i)G_b({\omega }_{i+h}){\tilde{n}}_a^{\ast }({\omega }_i){\tilde{n}}_b({\omega }_{i+h})\frac{2\pi }{T_{obs}}.$$

$$\big{<}|N|^{2}\big{>}=\sum_{i}|W(\omega_{i})|^{2}\,\frac{\omega_{0}^{-4}|\tilde{n}_{o}|^{4}\big{(}2\pi/T_{obs}\big{)}^{2}}{\big{(}\Omega^{2}+\Gamma^{2}\big{)}^{2}}\,$$

$$\big{<}|N|^{2}\big{>}=\int{W^{2}(\omega)\,\omega_{0}^{-4}\big{(}\Omega^{2}+\Gamma^{2}\big{)}^{-2}|\tilde{n}_{o}|^{4}\bigg{(}\frac{2\pi}{T_{obs}}\bigg{)}\,d\omega}\,$$

$$W(\omega )\propto \frac{\hat{\omega }}{4{\hat{\omega }}^2+{\Gamma }^2}$$

$${\int }_{{\omega }_a-5\Gamma }^{{\omega }_a+5\Gamma }{\tilde{a}}^{\ast }(\omega )\tilde{b}(\omega +\Omega )\frac{\omega -{\omega }_a}{4(\omega -{\omega }_a{)}^2+{\Gamma }^2}d\omega .$$

$$\int \frac{{\hat{\omega }}^2}{(4{\hat{\omega }}^2+{\Gamma }^2{)}^2}d\tilde{\omega }\approx \frac{\pi }{16\Gamma },$$

$$\big{<}|N|^{2}\big{>}^{1/2}=\frac{\pi}{\sqrt{8}}\omega_{0}^{-2}\Gamma^{-1/2}\big{(}\Omega^{2}+\Gamma^{2}\big{)}^{-1}|\tilde{n}_{o}|^{2}T_{obs}^{-1/2}\,.$$

$$\frac{|S|}{\big{<}|N|^{2}\big{>}^{1/2}}=\frac{1}{\sqrt{2}}\,|\hat{\lambda}|\,\bigg{(}\frac{\omega_{0}}{\Gamma}\bigg{)}^{1/2}\,\bigg{(}\omega_{0}T_{obs}\bigg{)}^{1/2}\,$$

$$\displaystyle SNR_{1-pair}\sim\frac{1}{\sqrt{2}}\,|\hat{\lambda}|\,Q^{1/2}\bigg{(}\omega_{0}T\bigg{)}^{1/2}$$

$$\displaystyle\sim 0.18\bigg{(}\frac{|\hat{\lambda}|}{10^{-5}}\bigg{)}\bigg{(}\frac{Q}{10^{3}}\bigg{)}^{1/2}\bigg{(}\frac{\nu_{0}}{3.3{\rm mHz}}\bigg{)}^{1/2}\bigg{(}\frac{T_{obs}}{1\,{\rm yr}}\bigg{)}^{1/2}$$

$$\sigma (\theta ,\varphi ,t)=\max \{0,\hat{n}\cdot \hat{r}\};$$

Full text

Sun: helioseismology – oscillations –magnetic fields

1 Introduction

Understanding the generation and distribution of the Sun’s interior magnetic field (B-field) is a major and longstanding goal of solar physics. In standard dynamo models, large-scale fluid flows –especially differential rotation in the shear layer near the bottom of the convection zone and the meridional flow from equator to poles –play central roles in generating and transporting magnetic field; see, e.g. Miesch (2005). Additionally, it has long been hypothesized that the Sun’s core contains a very strong toroidal B-field–possibly as large as $\sim 10^{6}G$–which is left over from the Sun’s formation, and which has persisted because of the extremely high conductivity in the Sun’s core; see, e.g., Gough & Taylor (1984).

Much of our knowledge of the Sun’s interior comes from helioseismology. Analysis methods for helioseismology data generally fall into one of two broad classes: local helioseismology and global helioseismology. The method proposed in this paper falls within the category of global helioseismology. Most (though certainly not all) work in global helioseismology has focused on the information encoded in the Sun’s oscillation frequencies. However, in this paper we focus on information about the Sun’s interior B-field that might be gleaned from measuring the mode shapes; i.e, the eigenfunctions as opposed to the eigenfrequencies.

In this paper we concentrate on couplings between modes that have the same $(n,l)$ values and $m-$values that are separated by a small (non-zero) integer. There are several motivations for this: a) sunspots tend to aggregate at preferred, evenly spaced longitudes (between one and four longitudes depending on the particular epoch and the northern versus southern hemispheres), which suggests the presence of substantial low-$m$ power in the B-field distribution, b) modes with the same $(n,l)$ values and neighboring $m$-values are nearly degenerate in frequency, and near-degeneracies amplify mode-mixing, and c) couplings between modes with different $m$-values are clearly not due to axisymmetric flows.

We analyze magnetically-induced mode mixing in two stages of increasing complexity. In the first stage, as a warm-up problem, we neglect the effects of mode damping. In the second stage, we account for realistic levels of mode damping. Our analysis motivates the introduction of a new, mode-mixing detection statistic that tests for the presence of non-axisymmetric mode-mixing in solar Doppler data. A nice feature of our detection statistic is that it is quite insensitive to the effects of spatial aliasing.

We estimate our statistic’s expected signal-to-noise ratio, as a function of the amplitude of the mode-mixing (which in turn depends on the amplitude and distribution of the Sun’s interior B-field). While we find that B-induced Solar mode mixing is likely too weak to be measurable with a single mode pair, we argue that combining results from a large number of mode pairs should dramatically increase the power of the method. The basic reason is that mode-mixing has a phase as well as an amplitude, so calculating our detection statistic produces a complex number. As we argue below, it seems reasonable to assume that over large regions in the space of mode-pairs, the ”signal contribution” to the statistic should produce nearly the same phase. So if one simply adds up these complex values for a large number of similar modes, the SNR for the mode-summed version of our detection statistic should grow roughly as $N_{pair}^{1/2}$, where $N_{pair}$ is the number of relevant mode pairs. With the mode-summed version of our detection statistic, we find that B-induced mode mixing should be measurable for a fairly wide range of plausible B-field magnitudes and geometries.

Of course, many of the ideas and methods in this paper already exist in the substantial literature on helioseismology. E.g., a formalism for using eigenmode-mixings to measure the Sun’s meridional flow, which emphasizes the importance of near-degeneracies, was developed and applied to SOHO data in a series of papers by Schad et al. (2011a, b, 2012, 2013). Also, inference from eigenfunction mixing is clearly closely related to inference from the analysis of cross-spectra, and the latter has been used to investigate meridional and zonal flows in a large number of papers by Woodard, Schou and others; see, e.g., Woodard et al. (2012, 2013). And of course there has been a good deal of other work exploring what helioseismology can say about the Sun’s interior B-field. E.g., Gough & Thompson (1990) considered the effects of an axisymmetric magnetic field that is misaligned from the Sun’s spin axis, while Antia et al. (2013) began the project of inferring the B-field in the convection zone from helioseismologically measured variations in the Sun’s local angular velocity. What is novel in our paper, we believe, is i) our emphasis on probing a truly non-axisymmetric B-field (not just a mis-aligned one) using measurements of eigenfunction-mixing between modes with different $m$ values, ii) our analysis of how results from nearly-degenerate perturbation theory get modified when realistic damping times for Solar modes are taken into account, and iii) our mode-mixing detection statistic, which we designed specifically to take realistic damping into account and to minimize the influence of spatial aliasing.

Finally, we mention that based on recent, asteroseismological measurements by the Kepler satellite Fuller et al. (2015) suggest that $\sim 20\%$ of red giants in the Kepler sample have very strong core B-fields: B\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10^{5}\,G and up to $\sim 10^{7}\,$G. However the relevance of these new discoveries to our own Sun is questionable, since at present the strong core B-fields are limited to red giants whose main sequence progenitors would have had masses in the range $1.6-2.0M_{\odot}$, and therefore would have had convective cores where the B-fields could have been generated by dynamo action; see Stello et al. (2016).

The rest of this paper is organized as follows. In Sec. 2 we summarize some necessary background material, partly to establish notation. In Sec. 3 we analyze the mixing of Solar oscillation modes with different $m$ values, resulting from small, non-axisymmetric perturbations. As a warm-up problem, in Sec. 3.1 we first examine this mixing in the limit of zero mode damping. In Sec. 3.2 we analyze Solar mode mixing for realistic mode damping rates and introduce our new mode-mixing detection statistic. In Sec. 3.2 we also estimate the signal-to-noise ratio (SNR) for this statistic, applied to a single mode pair, as a function of the strength of the B-field’s coupling between those modes. We find that mode-mixing is likely not observable for a single mode pair. However, we argue in Sec. 3.3 that combining the results from a large number of mode pairs could yield a detectable $SNR$ for an interesting range of B-field distributions. And while it might naively be supposed that spatial aliasing could significantly degrade our statistic’s detection power, in Sec. 3.4 we show why this is not the case, both analytically and via a (simplified) numerical simulation. In Sec. 4 we summarize our results, call attention to a couple caveats, and very briefly discuss next steps.

2 Background material

2.1 The signal and the noise

Our sign convention for modes is that a mode with frequency $\omega$ and azimuthal number $m$ has dependence $e^{i\omega t-im\phi}$ (with the Sun rotating in the positive-$\phi$ direction). Then to see the effect of rotation on mode frequencies, consider a mode that in the Sun’s co-rotating frame has the form $e^{i\omega_{co}t-im\tilde{\phi}}$, where $\tilde{\phi}\equiv\phi-\Omega t$ is the co-rotating azimuthal coordinate. The main effect of rotation on mode shapes and frequencies is a kinematic one: the mode gets dragged forward by the Sun’s rotation. Then in the inertial frame the mode has (approximately) the form $e^{i\omega_{co}t-im(\phi-\Omega t)}$, so the mode frequency in the inertial frame becomes $\omega\approx\omega_{co}+m\Omega$.

We will focus on the information encoded in the Doppler signal, $\delta\nu/\nu=(\vec{v}/c)\cdot\hat{n}$, where $\nu$ is the center of some absorption line, $\vec{v}$ is the velocity of the fluid at the Solar surface (more correctly, at the photosphere for the particular line of interest), and $\hat{n}$ is the line of sight to the detector. We denote this Doppler signal as $y=y(\theta,\phi,t)$, where $(\theta,\phi)$ are (inertial-frame) spherical coordinates on the Sun’s surface. The (time-dependent part of) the fluid velocity $\vec{v}$ at the Sun’s surface (radius $R$) can be decomposed as a sum of modes:

[TABLE]

where the $\vec{\xi}_{nlm}$ are the eigenmodes. (The $n=0$ modes are called f-modes, while the $n\geq 1$ acoustic modes are called p-modes. Solar g-modes have not been definitively detected yet, and we will not consider them in this paper.) For the most energetic modes, the surface velocity is dominated by its radial part, $v^{r}$, and for a spherically symmetric Sun, the radial components of the modes $\xi^{r}_{nlm}$ are $\propto Y_{lm}(\theta,\phi)$.

We will find it useful to separate the Sun’s modes into two classes: $\sim 10^{7}$ resolvable modes, and then all the unresolvable ones. For our purposes, the unresolvable modes form a ”confusion noise” background, in which the resolvable foreground modes are imbedded. Besides unresolved oscillation modes, there are several noise sources that one might call ”seeing noise”: finite photon counts, Earth atmospheric effects, pixelization, etc. (And there are other solar phenomena that are often not included in one’s helioseismic model, such as sunspots.) For the current study, it seems reasonable to assume that surface motions from unresolved modes dominate the ”noise”. And because the number of unresolved modes is vast, by the Central Limit Theorem that noise should be approximately Gaussian.

2.2 Coupling of nearly degenerate modes: a simple, two-mode illustration

We emphasize that the main physical reason that the Sun’s B-field can induce significant coupling is because modes that have the same $(n,l)$ values and nearby $m$ values are nearly degenerate in frequency. Therefore in this section we review some basic results from nearly-degenerate perturbation theory. The effect is simplest to understand for non-dissipative systems. For zero damping the perturbed fluid equations can be written in Hamiltonian form, so consider the following 2-D Hamiltionian:

[TABLE]

(with $\delta$ real), whose eigenvectors are obviously $(1,0)$ and $(0,1)$. Now we add a small, off-diagonal perturbation:

[TABLE]

For simplicity, we will consider just the two limiting cases where either $|\epsilon|<<|\delta|$ or $|\epsilon|>>|\delta|$. Case (i): $|\epsilon|<<|\delta|$. Then the two new eigenvalues are approximately $1+\delta+|\epsilon|^{2}/\delta$ and $1-|\epsilon|^{2}/\delta$ and the corresponding (unnormalized) eigenvectors are approximately $(1,\epsilon^{*}/\delta)$ and $(-\epsilon^{*}/\delta,1)$, respectively. E.g., if $\epsilon\sim 10^{-4}$ and $\delta\sim 10^{-3}$, then while the eigenvalues are perturbed by only $\sim|\epsilon|^{2}/\delta\sim 10^{-5}$, the components of the eigenvectors are shifted by $\sim\epsilon/\delta\sim 10^{-1}$. Case (ii): $|\epsilon|>>|\delta|$. The two new eigenvalues are approximately $1+\delta/2+|\epsilon|$ and $1+\delta/2-|\epsilon|$ and the corresponding (unnormalized) eigenvectors are approximately $(\epsilon/|\epsilon|,1-0.5\delta/|\epsilon|)$ and $(1-0.5\delta/|\epsilon|,-\epsilon/|\epsilon|)$, resp. That is, while the shift in the eigenvalues is only $\emph{O}(\epsilon)$, the components of the eigenfunctions are shifted by $\emph{O}(1)$.

2.3 B-field terms in the equations of motion, and their expansion into spherical harmonics

Let $v^{i}$ be the perturbed (time-dependent) velocity field. Very schematically, the interior B-field adds an extra ”forcing” term to the perturbed fluid equations of motion Schnack (2009):

[TABLE]

where $\gamma^{i}_{j}$ is a sum of terms of the general form $\nabla B\,\nabla B$ (with two of the indices contracted; e.g., $\nabla_{k}B^{i}\,\nabla_{k}B_{j}$) and $B\nabla\nabla B$. Similarly $\kappa^{i}_{jk}$ is a sum of terms of the form $B\,\nabla B$ (e.g., $B^{i}\,\nabla_{k}B_{j}$) and $\eta^{i}_{jkl}$ is a sum of terms of the general form $eBB$, where $e_{ij}$ is the spatial 3-metric. The tensors $\gamma^{i}_{j}$, $\kappa^{i}_{jk}$, and $\eta^{i}_{jkl}$ can all be expanded in tensor spherical harmonics: basically tensorial versions of the $Y_{lm}$’s. Based on the preferred longitudes for sunspots, it seems a reasonable guess that the largest (non-axisymmetric) couplings will come from low-$m$ terms: $m=\pm 1,\pm 2,\pm 3,{\rm or}\pm 4$. Recall, from the standard rules of composition of angular momentum, that a perturbation $\propto Y_{lm}$ (with $m$ non-negative, say) can couple two modes $\propto Y_{l^{\prime}m^{\prime}}$ and $\propto Y_{l^{\prime\prime}m^{\prime\prime}}$, respectively, only if (i) $|m^{\prime\prime}-m^{\prime}|=m$, (ii)$|l^{\prime}-l^{\prime\prime}|<l<l^{\prime}+l^{\prime\prime}$, and (iii)$l+l^{\prime}+l^{\prime\prime}$ is even ( with the last condition essentially representing parity conservation). The lowest order spherical harmonic connecting an $(n,l,m)$ and $(n,l,m+1)$ mode-pair is a $Y_{21}$ term (since the coupling from a $Y_{11}$ term must vanish by condition (iii) ), so it seems likely that this would usually be the dominant magnetic-perturbation term Similarly, it seems likely that the dominant magnetic perturbation that connects modes $(n,l,m)$ and $(n,l,m+2)$ modes is a $Y_{22}$ term. Of course, the amplitude of B-field forcing terms will have some radial dependence, which should be reflected in the relative strengths of the mode couplings for different $n$ values.

3 The magnitude of eigenmode mixing

In this section we estimate the magnitude of B-field-induced eigenmode mixing in the Sun. We do this in two steps: first, as a warm-up problem, in the limit of zero mode damping, and afterwards accounting for realistic damping rates. Throughout this section we will analyze the coupling between mode pairs ”as if no other modes existed”, but, after some development, it will be easy to see that including all the modes in our analysis would have a negligible effect on any particular two-coupling. The basic reason is that dimensionless mode-coupling parameters will typically have magnitude $\hat{\lambda}\sim 10^{-6}$ (or less), and the influence of any other modes on the coupling of some given mode-pair will be quadratic in these coupling parameters.

3.1 Mode-mixing estimate for zero damping

Observed p-modes typically have quality factors $Q\sim 10^{3}$, and therefore have damping times of order days. This observation motivates making some initial, ballpark estimates in the limit of zero damping. (We will find that our estimates based on this limit are not at all accurate, but they will nevertheless provide an interesting basis for comparison.) Here we will also make several other simplifying approximations. First, we will approximate the Sun as uniformly rotating at frequency $\Omega$, so the Sun’s ”frozen in” $B$-field is time-independent in the Sun’s rotating frame. Therefore in this subsection we find it simplest to work in this rotating frame, since it allows us to use time-independent perturbation theory to estimate the amplitude of mode-mixing.

In an inertial rest frame, the frequency splitting of Solar modes with the same $(n,l)$ values, but different $m$ values, is dominated by the kinematic term: $\Delta\nu\approx(\Omega/2\pi)\,\Delta m$. However, in the rotating frame of the Sun (and ignoring for the moment the effect of the B-field), the mode-splitting is primarily due to the Coriolis effect: $\Delta\nu\approx c_{nlm}(\Omega/2\pi)\,\Delta m$, where the size of the coefficient $c_{nlm}$ is typically $\sim 1\%$ (Kosovichev (1996)). Thus the $\delta$ term in our simple example, Eqs. (4)–(7), would in this case be $\delta\sim|\Delta\nu/\nu|\sim 10^{-6}\,$ (for $|\Delta m|=1$). Thus, for the idealized case of zero damping, magnetic cross-coupling shift frequencies by only $\sim 10^{-6}$ (fractionally) could still lead to ${\cal O}(1)$ mixing between the eigenfunctions.

3.2 Mode-mixing estimate for realistic damping

We will now see that the mode-mixing estimates we made in the previous subsection are substantially altered when we incorporate realistic mode damping rates. We will continue to approximate the Sun’s rotation as uniform, but while in Sec. 3.1 we worked in the Sun’s co-rotating frame, in this section we find it more convenient to work in the Sun’s inertial frame. Also in this subsection we will neglect spatial aliasing–effectively assuming that we are observing all $4\pi$ steradians of the Solar surface; we will analyze the effects of spatial aliasing in Sec. 3.4.

We return to our two-mode system, but this time described by the Langevin equation, which includes both damping and driving terms in addition to a mode-mode coupling term. For simplicity we will consider a pair of modes with quantum numbers $(n,l,m)$ and $(n,l,m+1)$. The extension to couplings between modes with $|\Delta m|=2,3,{\rm or}\,4$ is trivial, and at the end of this subsection we will describe how to modify our equations to account for different values of $\Delta m$.

Call the complex amplitudes of the two modes $a(t)$ and $b(t)$, respectively. The Langevin equations for the coupled modes are then

[TABLE]

Here $\omega_{a,b}$ are the two modes’ oscillation frequencies, $\Gamma_{a,b}$ are their damping rates, and the $n_{a,b}(t)$ are driving terms that have the statistical characteristics of noise. Since $\lambda$ is always small compared to $\omega_{a,b}^{2}$, we will treat it as a small perturbation. Fourier transforming the solutions to Eqs. (9)-(10) and expanding them through first order in $\lambda$, we find:

[TABLE]

where the Green’s functions $G_{a,b}(\omega)$ are explicitly given by

[TABLE]

and where our convention for Fourier transforms is

[TABLE]

Now, as a starting point for introducing our mode-mixing detection statistic, consider the following integral:

[TABLE]

which by the convolution theorem equals

[TABLE]

To maximize the signal-to-noise, it will eventually prove useful to introduce an additional weighting factor $W(\omega)$ in the integrand of (17), so the modified version will be

[TABLE]

But for now we will stick with the version in Eq. (17). We will also eventually want to restrict the limits of integration in Eq. (18), but for now we just leave them unspecified. An expansion of Eq. (17) through first order in $\lambda$ yields:

[TABLE]

Next we will compare the relative sizes of the terms that are linear in $\lambda^{*}$ (i.e., the B-dependent, ”signal” terms”) with the sizes of the $\lambda$-independent terms (i.e., the B-independent, ”background noise” terms). It seems safe to approximate the noise as stationary, and to approximate its spectrum as flat over the very narrow region of interest. Stationarity implies that noise amplitudes at different frequencies are uncorrelated. Likewise $|\tilde{n}_{a}(\omega)|^{2}$ and $|\tilde{n}_{b}(\omega+\Omega)|^{2}$ should be approximately equal, so we will refer to both as $|\tilde{n}_{0}|^{2}$. Similarly, $\Gamma_{a}$ and $\Gamma_{b}$ should be nearly the same, so for simplicity, in our estimates we will also take them to be equal: $\Gamma_{a}=\Gamma_{b}\equiv\Gamma$. Using the fact that $\tilde{n}_{a}(\omega)$ and $\tilde{n}_{b}(\omega+\Omega)$ are uncorrelated, the expectation value of (19) becomes

[TABLE]

where $\omega_{0}\equiv(\omega_{a}+\omega_{b})/2$, $\hat{\omega}\equiv\omega-\omega_{a}$ and $\Delta\equiv(\omega_{b}-\omega_{b})/2=(\Omega+\epsilon)/2$ (so $\Delta$ is positive), and where we have defined $\hat{\lambda}\equiv\lambda/\omega_{0}^{2}$ (so $\hat{\lambda}$ is dimensionless). Again, $\epsilon$ is typically $\sim 1\%$ the size of $\Omega$, so $\Delta$ is nearly equal to $\Omega$. Then, neglecting terms that are cubic (or higher) in the small (compared to $\omega_{a}$) quantities $\hat{\omega}$, $\Omega$ and/or $\Gamma$, one easily shows that the sum \big{(}G_{a}(\omega)\,+G^{*}_{b}(\omega+\Omega)\big{)} in Eq. (20) becomes

[TABLE]

Similarly the term \big{[}G^{*}_{a}(\omega)G_{b}(\omega+\Omega)\big{]} can be approximated as

[TABLE]

and so the expectation value of the ”signal part” of our modified statistic (18) becomes

[TABLE]

We asserted at the beginning of this section that, while we have restricted attention to only two modes, including the effects of couplings to other modes would affect any 2-mode SNR only by higher order terms in the magnetic coupling parameters. At this point, an easy way to see that is just to add a third mode, with amplitude $c(t)$, to our dynamical system system Eqs. (9)–(10), and repeat our calculations down through Eq. (19). We leave that as an exercise for the reader, after which the generalization to an arbitrary number of modes should be obvious.

Again, we want to choose $W(\omega)$ to maximize the signal-to-noise, but before doing so, to simplify the analysis, we will make a couple more approximations. First, because of the \big{(}4\hat{\omega}^{2}+\Gamma^{2}\big{)} factor in the denominator of the integrand in Eq. (23), the integral will be dominated by the region \hat{\omega}\equiv(\omega-\omega_{a})\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}{\rm few}\times\Gamma. In this region, the rms value of $\hat{\omega}$ is $\sim\Gamma$ (and below we will see that it is this rms value that matters for the SNR) so in the factor \big{(}\hat{\omega}-\frac{i}{4}\Omega\Gamma\omega_{a}^{-1}\big{)} the ratio of the second term to the first is typically $\sim\Omega/\omega_{a}\sim 10^{-4}$, so we will neglect that second term. For the same reason, we can approximate $\omega_{0}^{2}\omega_{a}^{-3}$ by $\omega_{0}^{-1}$.

To further simplify the calculation, we will assume for the moment that one of the terms $\Gamma$ or $\Omega$ is significantly larger than the other; i.e., either $\Gamma>>\Omega$ or $\Omega>>\Gamma$. In either case, in the term

[TABLE]

the pieces \big{(}-4\Omega\hat{\omega}-2i\,\Omega\,\Gamma\big{)} can be neglected with respect to \big{(}\Omega^{2}+\Gamma^{2}\big{)}. So Eq. (23) has now been approximated as

[TABLE]

To choose the optimal $W(\omega)$, we next need an expression for the ”background noise” piece $N$, which arises from random correlations between the two modes.

We find $N$ easiest to estimate if we approximate the $\lambda$-independent part of the continuous integral (18) by the corresponding discrete sum over frequency bins, with bin width $\Delta\omega=2\pi/T_{obs}$:

[TABLE]

where $h$ is the integer nearest to $\Omega/\Delta\omega$. The terms $\tilde{n}^{*}_{a}(\omega_{i})$ and $\tilde{n}_{b}(\omega_{i+h})$ are statistically independent, so the sum accumulates like a random walk. Hence, using the same approximations as above, we have

[TABLE]

Converting this discrete sum back to an integral, we obtain

[TABLE]

Optimizing $W(\omega)$ means maximizing $|S|$ for fixed \big{<}N^{2}\big{>}^{1/2}. Using the method of Lagrange multipliers, one finds that the optimum choice is

[TABLE]

Thus the final form of our mode-mixing detection statistic (for $|\Delta m|=1$) is

[TABLE]

The limits of integration here are somewhat arbitrary, but since the integrand falls off rapidly for $|\omega-\omega_{a}|>>\Gamma$, the SNR should depend only weakly on the choice. We should explain that the reason we restrict the integration range at all is just our intuition that restricting the region of integration limits the possibility of contamination from any artifacts in the spectra (e.g., from data gaps, any instrumental lines, etc.) .

Plugging Eq. (29) into Eq. (25) and using the approximation

[TABLE]

(the rhs of Eq. (31) is actually the exact value of the integral when the limits of integration are taken to $\pm\infty$), we find that |S|=(\pi/4)|\hat{\lambda}|\omega_{0}^{-1}\Gamma^{-1}\big{(}\Omega^{2}+\Gamma^{2}\big{)}^{-1}|\tilde{n}_{o}|^{2} and that

[TABLE]

Hence we arrive at

[TABLE]

so the SNR for a single mode-pair is

[TABLE]

where $Q\equiv\omega_{0}/\Gamma$ and $\nu_{0}\equiv\omega_{0}/(2\pi)$ Note that while we started with the simplifying assumption that either $\Gamma>>\Omega$ or $\Omega>>\Gamma$, or final result Eq. (35) is independent of which of these limits we are in. It therefore seems reasonable to assume that Eq. (35) is also a fairly accurate estimate of the SNR for the intermediate case $\Gamma\approx\Omega$. Finally, while for simplicity of exposition, we so far have restricted to the $\Delta m=1$ case, the generalization to $\Delta m=2,3,4$ is trivial: in every numbered equation in this subsection, just replace $\Omega$ by $\Delta m\,\Omega$.

3.3 Summing the signal over mode pairs

At first glance, the estimate (35) does not seem very promising, unless the B-field is near the upper range of expectations. But Eq. (35) represents the SNR for just a single pair of modes, and it seems likely that the situation improves quite dramatically when one combines results from many mode pairs with the same $\Delta m$. The reason is that the B-field coupling should be nearly phase-coherent over large numbers of mode pairs. To be concrete, assume, e.g., that the dominant coupling for mode pairs with $\Delta m=1$ comes from a $Y_{2}^{1}$ tensor perturbation. This is a large-angular-scale perturbation, and so it seems likely that the phase of the perturbation is also coherent over a large range of radii. (Of course, as the Sun rotates and drags the B-field with it, the complex amplitude of the perturbation rotates in the complex plane at the rate $\Omega$.) But then the signal part of our complex detection statistic should have (approximately) the same phase for some large number of mode pairs. Therefore if we add up the complex amplitudes of our detection statistic over mode pairs, the signal parts add coherently, and the total SNR (over all pairs) scales like $N_{pair}^{1/2}$, where $N_{pair}$ is the number of (phase-coherent) pairs. There are of order $10^{7}$ measured modes, and so about $10^{7}$ neighboring mode pairs, $(n,l,m)$ and $(n,l,m+1)$, and so the SNR enhancement factor could be as high as $\sim 3000$. That’s a rough upper limit on the enhancement factor. And clearly, the extent to which this possible enhancement is realized depends in part on the radial scale over which the angular pattern of the B-field changes substantially. But even plugging in $N=10^{4}$, with $Q=3\times 10^{3}$ and $T_{obs}=4\,$yr, one finds B-field couplings with $|\hat{\lambda}|$ as small as $2\times 10^{-6}$ could be measured with $SNR\approx 12$- which is much more promising! Based on the fact that sunspots cluster at up to $4$ longitudes (at any time), it seems quite plausible that our method could yield information on spherical harmonics of the forcing tensor, Eq. (8), up to $|\Delta m|\approx 4$.

3.4 Robustness of our mode-mixing detection statistic against spatial aliasing

In our formula for the SNR of our 2-mode detection statistic, Eq. (35), the amplitudes $a(t)$ and $b(t)$ are the true amplitudes of the modes corresponding to spherical harmonics $(l,m)$ and $(l,m^{\prime})$, respectively. But to date all helioseismology observations have been made either on the Earth or from satellites whose distance from the Earth is a very small fraction of $1\,$AU. So with current telescopes we only have access to the half of the Sun facing us, and, when one accounts for the fact that we only measure the component of the Sun’s surface velocity that is along the line of sight, one finds that our effective viewing area is closer to one-third of the Sun’s surface. This is the origin of spatial aliasing – a mixing of the measured spherical harmonics. Spatial aliasing is a rather large effect at any instant, so one might worry that it will swamp the mixing due to the Sun’s interior B-field. The aim of this subsection is to show that, in fact, spatial aliasing does not substantially degrade the power of our detection statistic. First we will show analytically why this is the case. We have also performed some simple simulations of the effects of aliasing, and we will show that our simulation results are consistent with our analytic estimates. Let $(x,y,z)$ be inertial coordinates, with origin at the center of the Sun and $z$ along the Sun’s spin axis, and then define $(\theta,\phi)$ on the surface of the Sun in the usual way: $cos\theta=z/r$, $sin\theta\,cos\phi=x/r$ , and $sin\theta\,sin\phi=y/r$, with $r=\sqrt{x^{2}+y^{2}+z^{2}}$. (Here and below we shall neglect the Sun’s slight oblateness; i.e., we model its surface as a sphere.) Let $\hat{n}(t)$ be the unit vector from the center of the Sun to the observer (here assumed to be on or near the Earth), and let $\hat{r}$ be the unit radial vector from the Sun’s origin to the location $(\theta,\phi)$. On the Sun’s surface, define $\sigma(t)$ by

[TABLE]

i.e., $\sigma(\theta,\phi,t)$ is $\hat{n}\cdot\hat{r}$ on the ”front side” of the Sun (facing the observer) and zero on the ”back side”. Define the time-dependent inner product between any two complex functions $f(\theta,\phi,t)$ and $g(\theta,\phi,t)$ on the Sun’s surface by:

[TABLE]

Next we define the overlap function $\gamma_{l,m}^{m^{\prime}}(t)$ by

[TABLE]

where \big{(}N_{lm}N_{lm^{\prime}}\big{)}^{-1/2} is just a time-averaged normalization factor, specifically:

[TABLE]

The overlap function $\gamma_{lm}^{m^{\prime}}(t)$ varies on the timescale of a year because that is the timescale on which the ”visible half” of the Sun (from the Earth) varies. Also, clearly, $\gamma_{lm}^{m^{\prime}}(t)=\gamma_{lm^{\prime}}^{m\,*}(t)$.

Now, continuing with our 2-mode example, imagine for simplicity that the Sun’s perturbed radial velocity $\delta v^{r}$ at the surface $R$ is the sum of only two modes:

[TABLE]

where the amplitudes (but not the phases) of $a(t)$ and $b(t)$ are slowly varying (due to damping and excitation) and where, again, $\omega_{b}-\omega_{a}\approx(m^{\prime}-m)\,\Omega$. Now let us define the observed values of these amplitudes, $a_{ob}(t)$ and $b_{ob}(t)$ by

[TABLE]

so that

[TABLE]

Transforming to the Fourier domain, we then have

[TABLE]

Now let us see what happens if in Eq.(18) we replace $\tilde{a}(\omega)$ and $\tilde{b}(\omega)$ by the corresponding ’observed’ amplitudes, given by Eqs. (44)-(45). The ”background noise” piece

[TABLE]

gets augmented by one cross-term that is quadratic in $\tilde{n}_{a}$:

[TABLE]

as well as a similar term that is quadratic in $\tilde{n}_{b}$. The convolution factors $\tilde{\gamma}_{lm}^{m\,*}(\omega-\omega^{\prime})$ and $\tilde{\gamma}_{lm^{\prime}}^{m}(\omega-\omega^{\prime\prime})$ effectively ”smear” $\tilde{n}^{*}_{a}(\omega)$ and $\tilde{n}^{*}_{a}(\omega+\Delta m\,\Omega)$ over frequency bands of width $\sim 2\pi/{\rm yr}$. However the key point is that $\Omega>>2\pi/{\rm yr}$, so even these smeared bands are non-overlapping, and hence their values are not correlated – unlike the case the for signal terms. It is easy to see that the same is true for the aliasing terms quadratic in $\tilde{n}_{b}$. So aliasing modifies the background noise piece by only a small fraction of its value. To illustrate how this works, Fig. 1 simulates the build-up over frequency of three pieces of the our detection statistic: the $\lambda$-independent noise piece given by (the continuous version of) Eq.(26), the $\lambda$-dependent signal piece given by Eq.(25), and the aliasing contribution given by the sum of Eq.(47) and the corresponding term that is quadratic in $\tilde{n}_{b}$. The specific parameters chosen for his particular simulation were $l=10,m=4,m^{\prime}=5$, $T_{obs}=1/$yr, $\omega_{0}=0.02s^{-1}$, $\Omega=3\times 10^{-6}{s}^{-1}$, $\Gamma=10^{-3}\omega_{0}$, and $\omega_{a,b}=\omega_{0}\pm\Omega/2$. We also take $n_{a}$ and $n_{b}$ to be uncorrelated white noises with the same amplitude, which should be a good approximation over the narrow frequency band of interest. For the sake of visual clarity, we took $\hat{\lambda}=3\times 10^{-4}$, with is probably unphysically large, but does not affect the relative sizes of the ”unaliased” background noise piece and the aliasing piece. The take-home point of Fig. 1 is that the aliasing contribution is only a modest fraction of the full noise, and so has little effect on our SNR estimates. Fig. 1 displays just one realization of $n_{a}$ and $n_{b}$, but is a typical result.

To summarize, if one Fourier transforms the observed mode amplitudes, aliasing adds sidebands to the mode ”lines”. But for integration times of order a year or more, the separation of these sidebands from the carrier is much less than the separation between the mode lines, and so aliasing does not induce the mode correlations that our detection statistic searches for. We suggest that this method of mitigating the effects of spatial aliasing might prove useful in other sorts of helioseismological studies as well.

4 Summary, caveats and future work

In this paper we investigated the possibility of using mode-mixing to probe the Sun’s interior B-field. We concentrated on modes with the same $(n,l)$ values and nearby $m$ values, since such modes are nearly degenerate in frequency, which enhances mode-mixing. We constructed a novel mode-mixing detection statistic for this effect, Eq. (30), and we showed that for long observation times (of order a year or more), our statistic is quite robust against the effects of spatial aliasing. We estimated the SNR for our detection statistic, for realistic mode damping. The detectability of mode-mixing is enhanced by a couple of factors, in addition to the near-degeneracies. First, the SNR grows like $N_{cyc}^{1/2}$, where $N_{c}=\nu_{0}T_{obs}$ is the number of observed oscillation cycles. For a year’s worth of observation of five-minute oscillations, $N_{cyc}^{1/2}\sim 300$. Second, we argued that the phase of the mode-mixing is likely approximately constant over a large range of $(n,l)$ values, for fixed $\Delta m$. So by adding up the complex values of our detection statistic Eq. (30) over $N_{pair}$ mode pairs, the total SNR should grow roughly as $N_{pair}^{1/2}$, where $N_{pair}$ is the number of mode pairs with similar mixing angle. $N_{pair}$ could be as large as $\sim 10^{7}$ , but even assuming $N_{pair}\sim 10^{4}-10^{5}$, Eq. (35) suggests that couplings as small as |\hat{\lambda}|\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10^{-6} should be detectable.

As caveats, we here remind the reader of some effects that we have not yet taken into account in our analyses. First, for simplicity, all our analyses took the Sun to be uniformly rotating. It would be more accurate to describe the Sun’s angular velocity as a constant $\Omega_{0}$ (some weighted average of the angular velocity field) plus an axisymmetric perturbation $\delta\Omega(r,\theta)$. However we would argue that this improvement would not substantially affect our estimate of the mode-mixing SNR. The reason is that, as we have seen, the physical mechanism that is most important for causing mode-mixing to saturate is mode damping, and typically $\Gamma>\delta\Omega$. I.e., on the damping timescale $\Gamma^{-1}$, the non-uniform part of the angular velocity, $\delta\Omega$, is too small to cause much ”re-arrangement” of fluid and magnetic field inside the Sun. But we have not actually demonstrated this, so include that as a caveat.

Secondly, and perhaps most importantly, we have not yet tried to assess the likely impact of systematic errors on measurements of mode mixing. There are quite a few known instrumental effects, such as pixelization, whose impact we could reasonably try to assess. However, inversions of helioseismology data today also reveal effects that are clearly spurious but of unknown origin, such as the infamous ”center-to-limb” effect, see  Duvall & Hanasoge (2009); Baldner & Schou (2012); Zhao et al. (2013). It is hard to assess the impact of systematics that are not understood, which is a problem that our proposed method shares with much of the rest of helioseismology.

Regarding future work, Cutler & Woodard have recently begun to calculate the summed version of our mode-mixing detection statistic using SDO/HMI data.

Acknowledgments

This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract to the National Aeronautics and Space Administration. Special thanks go to Martin Woodard for a great many helpful and informative discussions, general encouragement, and for carefully reading and improving a draft of this manuscript. I also owe special thanks to Marco Velli and Neil Murphy for getting me involved in this subject and for many useful discussions and overall encouragement. Also, Stuart Jeffries, Douglas Gough, Charles Baldner and Tim Larson were all very generous with their time in helping educate me about this subject.