Analyzing supergranular power spectra using helioseismic normal-mode coupling
DDraft version February 18, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Analyzing supergranular power spectra using helioseismic normal-mode coupling
Chris S. Hanson, Shravan Hanasoge,
2, 1 and Katepalli R. Sreenivasan
3, 1 Center for Space Science, NYUAD InstituteNew York University Abu DhabiAbu Dhabi, UAE Department of Astronomy and AstrophysicsTata Institute of Fundamental ResearchMumbai, India Department of Physics, Courant Institute of Mathematical Sciences, Tandon School of EngineeringNew York University, New York, 10003 (Received February 18, 2021; Revised; Accepted)
Submitted to APJABSTRACTNormal-mode coupling is a technique applied to probe the solar interior using surface observationsof oscillations. The technique, which is straightforward to implement, makes more use of the seismicinformation in the wavefield than other comparable local imaging techniques and therefore has the po-tential to significantly improve current capabilities. Here, we examine supergranulation power spectrausing mode-coupling analyses of intermediate-to-high-degree modes by invoking a Cartesian-geometricdescription of wave propagation under the assumption that the localized patches are much smaller insize than the solar radius. We extract the supergranular power spectrum and compare the results withprior helioseismic studies. Measurements of the dispersion relation and life times of supergranulation,obtained using near surface modes (f and p ), are in accord with the literature. We show that thecross-coupling between the p and p acoustic modes, which are capable of probing greater depths, arealso sensitive to supergranulation. Keywords:
Methods: data analysis – Sun: helioseismology – convection – waves INTRODUCTIONHelioseismology is applied to image the solar inte-rior through the analysis of the Sun’s observed sur-face wave field (e.g. Christensen-Dalsgaard 2002). Theoscillations, which comprise acoustic p modes andsurface-gravity f modes, are thought to be generatedby near-surface turbulent stresses (Goldreich & Keeley1977). Spherically symmetric models (e.g. Christensen-Dalsgaard et al. 1996) are often used to compute modeeigenfunctions. Inhomogeneities and 3D structure in theSun, such as flows and magnetic fields, cause the solareigenfunctions to be different from those of the referenceset. However, the reference eigenbasis has properties of
Corresponding author: Chris S. [email protected] completeness and orthogonality, allowing us to write so-lar eigenfunctions as weighted linear combinations of thereference set. Thus, the normal modes of the Sun aresaid to be “coupled” with respect to the reference, andmay be used to probe non-axisymmetric properties ofthe solar interior (Woodard 1989).Over the past fifty years, a number of helioseismictechniques have been developed to image the subsurface.Some of the most successful results include inferencesof the Sun’s radial and latitudinal differential rotation(Schou et al. 1998) using global-mode analysis, sectoralRossby waves (L¨optien et al. 2018) through ring-diagramanalysis (Hill et al. 1996), and meridional circulation(Gizon et al. 2020) using time-distance helioseismology(Duvall et al. 1993b). Each of these techniques tends toutilize one component of the wavefield, e.g. frequencyshifts or travel times, whose estimation often requiressignificant data processing. The technique of normal- a r X i v : . [ a s t r o - ph . S R ] F e b Hanson, Hanasoge and Sreenivasan mode coupling (Woodard 1989), on the other hand, con-tains straightforward prescriptions for retrieving seismicmeasurements from observations, i.e., weighted linearsums over Fourier-domain wavefield correlations, whichare then ready for making inferences of subsurface per-turbations. The technique has proven successful in de-tecting and inferring global-scale perturbations such asRossby waves and toroidal flows (e.g. Hanasoge & Man-dal 2019; Hanasoge et al. 2020). While mode-couplingtheory for local helioseismic analysis has been developed(Woodard 2006), with some inferences (Woodard 2007),validation and detailed comparisons to other techniqueshas yet to be performed.The local helioseismic application of mode coupling as-sumes the observed wave field, over a small patch of thesolar surface, is described by an orthonormal and com-plete basis of eigenfunctions ξ α ( x ), where x = ( x, y, z )describes the Cartesian space coordinate . The descrip-tion in Cartesian geometry is valid in the limit where theobserved patch size used in the analysis is much smallerthan the solar radius. Here we gather the mode charac-teristics in the subscript α ≡ { n, k } which denotes theradial order n and horizontal wavenumbers k = ( k x , k y )of the mode. As mentioned earlier, subsurface pertur-bations in the Sun render the solar eigenfunctions ξ (cid:12) α different from those of the reference model ξ α . How-ever, the orthonormality of the reference eigenfunctionsallow us to express the former as a weighted linear sumof the latter, ξ (cid:12) α = (cid:80) α (cid:48) c α (cid:48) α ξ α . The coupling coefficients c α (cid:48) α are dependent on the unknown subsurface structuresin the sun and can be estimated from observations.Under the first-Born approximation, which assumessmall perturbations, previous studies have shown thatcoupling coefficients c α (cid:48) α are proportional to the cross-correlation of the complex Fourier components of theobserved surface wave field φ ω k = φ ( k , ω ) (e.g., Woodard2006, 2014). In Cartesian geometry, this relationship isexpressed as (cid:104) φ ω ∗ k φ ω + σ k + q (cid:105) = H ω | k || k + q | σ c α (cid:48) α (1)where the angular brackets define the expected value, φ ∗ is the complex conjugate of φ , and ( ω, k ) and ( σ, q )refer to the angular frequency and vector wavenumbersof the modes and perturbations, respectively. The co-efficient H contains information on the power-spectralmodel (see Eq. 4 below) and incorporates the scatteringphysics described by the first-Born approximation (e.g.Hanasoge et al. 2017). Although H does not explicitly x is oriented along the direction of rotation, y along the towardthe solar north pole and z pointed outwards along the local radialcoordinate contain the radial order n as a sub- or superscript, thedependence on n appears through the selection of thefrequency window used in computing the H .The observed wave field correlations, and in turn thecoupling coefficients, are dependent on the spatial struc-ture and temporal frequency of subsurface perturbation U . The relationship between the correlations and U isdetermined through the sensitivity kernels K , which arederived using the reference solar model, φ ω ∗ k φ ω + σ k + q = H ω | k || k + q | σ (cid:90) R (cid:12) K k , q ( z ) · U σ q ( z )d z + η, (2)where R (cid:12) is the solar radius and η is a noise real-ization. The left-hand side of Eq. 2 is trivially com-puted from the Fourier components of a small patch oftracked Doppler images, the standard data used in mostlocal helioseismic studies. Unlike ring-diagram analy-sis, which is limited to computing the mean flow field( σ, q ) = (0 ,
0) within the patch, mode coupling can beused to retrieve different spatial and temporal scales.Furthermore, the method does not require the compli-cated geometric point-to-point cross-correlation averag-ing schemes of time-distance helioseismology (e.g. Du-vall et al. 1993b; Gizon et al. 2010).In this Paper, we estimate the supergranulation powerspectrum using mode coupling (i.e. computing the left-hand side of Eq. 2) and compare the results from othertechniques. Supergranulation is a distinct feature seenin the Sun’s near-surface flow field and is regarded as anintermediate scale of thermal convection. Supergranulesare approximately 35 Mm ( | q | R (cid:12) ∼ nalyzing supergranular power spectra using helioseismic normal-mode coupling DATA ANALYSISWe use seven Carrington rotations (2197-2203) ofDopplergram images taken by the Helioseismic andMagnetic Imager (HMI, Schou et al. 2012). We trackedregions of size ∼ along the equator for 11days, from 70 ◦ East to 70 ◦ West. For each Carringtonrotation, we generate two data cubes where the centralmeridian crossing time for each cube coincides with theCarrington longitude of disk center located at 90 ◦ and270 ◦ . In total, we have 14 data cubes. The spatial reso-lution is downsampled by a factor of 2 and the images arePostel projected. The Dopplergrams were tracked at theSnodgrass (1984) rotation rate, which is 0.02893 µ rad/sslower than Carrington rotation rate. In practice, thetracking and generation of data cubes is performed us-ing mtrack . To avoid spatial aliasing, we apply a 2Dspatial-apodization function to the cubes, which is equalto 1 for pixels within 90 Mm of the patch center and ta-pers to zero over a distance of 7 Mm. The data cubesare then Fourier transformed in each dimension using aDiscrete Fourier Transform to generate φ ω k .Here, we compute B σ kq coefficients which are linear-least-squares fits to raw wavefield correlations (Woodard2016), B σ kq = (cid:80) ω H ω ∗ kk (cid:48) σ φ ω ∗ k φ ω + σ k + q (cid:80) ω |H ωkk (cid:48) σ | , (3)where k = | k | and k (cid:48) = | k (cid:48) | = | k + q | . The B σ kq co-efficients are equivalent to the integral of Eq. 2. Inthe above equation (and Eq. 1 and 2), H is the power-spectral model H ωkk (cid:48) σ = − ω ( N k | R ωα | R ω + σα (cid:48) + N k (cid:48) | R ω + σα (cid:48) | R ω ∗ α ) , (4)where N k is the mode amplitude and R ωα is theLorentzian profile of the mode resonance (Andersonet al. 1990; Duvall et al. 1993a), R ωα = 1( ω α − i Γ α / − ω . (5)Here, Γ α is the full width at half maximum of the modeand ω α is the resonant frequency. The form of Eq. 4 as-sumes the wave operator is Hermitian (Woodard 2014),and for simplicity, we have neglected leakage terms. Thereference mode width and frequency are taken from themode-fit parameter outputs of the ring-diagram module http://hmi.stanford.edu/rings/modules/mtrack.html rdfitc (Basu et al. 1999). We take fits from all tiles lo-cated at disk center from 2010 to 2018, and compute theaverage. The mode amplitude N k is computed through N k = (cid:80) ω | φ ω k | (cid:80) ω | R ωα | , (6)where the sum in ω is over the frequency bins withinfive line widths of ω nk .The computation of all possible B σ kq would populatea five dimensional array [ k x , k y , k x + q x , k y + q y , σ ] foreach possible n, n (cid:48) combination. This requires signif-icant computational and storage capacity. However,most of these couplings would be uninformative due ei-ther to the prominence of the background or becausethey are insensitive to the spatial scales of interest. Fur-thermore, due to Eq. 5, the coupling between modes de-creases as the frequency spacing of the modes increasesand the background begins to compromise the measure-ment. To ensure that we capture significant couplings,we compute B σ kq coefficients for modes sufficiently closeto each other in frequency and perform the sums in Eq. 3over a number linewidths. Specifically, | ω α − ω α (cid:48) − σ | ≤ (cid:15) Γ α , (7)and ω ∈{ ω α − (cid:15) Γ α / , ω α + (cid:15) Γ α / }∪ (8) { ω α (cid:48) − (cid:15) Γ α (cid:48) / , ω α (cid:48) + (cid:15) Γ α (cid:48) / } , where (cid:15) depends upon the n and n (cid:48) used. Here, (cid:15) isdetermined empirically to obtain reasonable signal-to-noise ratio. As a future study, it is useful to performparameter searches to determine the ideal (cid:15) . Finally, welimit ourselves to intermediate frequencies 2 ≤ ω/ π ≤ ∼ . ≤ ω/ π = 3 mHz, as well as the squared modulus of the B coefficients computed for an f mode of wave number | k | R (cid:12) = 900. Given the frequency limits of Eq. 7 and 8,the computed couplings follow the f mode ring. In thisstudy, we estimate couplings in the range of | q | ≤ | q | R (cid:12) ≈ NOISE MODEL http://hmi.stanford.edu/teams/rings/modules/rdfitc/v13.html Hanson, Hanasoge and Sreenivasan −1000 −500 0 500 1000k x R ⊙ −1000−50005001000 k y R ⊙ q x R ⊙ q y R ⊙ Figure 1.
Slice of Doppler power spectrum at fixed fre-quency ω/ π ≈ | B σ kq | computed using the couplings of an f mode [ k x R (cid:12) , k y R (cid:12) ] =[634 , − (cid:15) = 1. The measurement B σ kq contains a non-trivial sys-tematic correlated-noise component (e.g., Gizon 2004;Woodard 2007), which may be estimated through (cid:104)| B σ kq | (cid:105) = (cid:80) ω,ω (cid:48) H ω ∗ kk (cid:48) σ H ω (cid:48) kk (cid:48) σ (cid:104) φ ω ∗ k φ ω + σ k (cid:48) φ ω (cid:48) k φ ω (cid:48) + σ ∗ k (cid:48) (cid:105) ( (cid:80) ω |H ωkk (cid:48) σ | ) . (9)Invoking Isserlis’s theorem (Isserlis 1918), the expectedvalue , denoted by angle brackets in Eq. 9, is given by (cid:104) φ ω ∗ k φ ω + σ k (cid:48) φ ω (cid:48) k φ ω (cid:48) + σ ∗ k (cid:48) (cid:105) = (cid:104) φ ω ∗ k φ ω (cid:48) k (cid:105)(cid:104) φ ω + σ k (cid:48) φ ω (cid:48) + σ ∗ k (cid:48) (cid:105) + (cid:104) φ ω ∗ k φ ω + σ k (cid:48) (cid:105)(cid:104) φ ω (cid:48) k φ ω (cid:48) + σ ∗ k (cid:48) (cid:105) + (cid:104) φ ω ∗ k φ ω (cid:48) + σ ∗ k (cid:48) (cid:105)(cid:104) φ ω (cid:48) k φ ω + σ k (cid:48) (cid:105) . (10)We assume multivariate Gaussian noise and hence themodes are uncorrelated across different frequencies ω and wavenumbers k , i.e., (cid:104) φ ω ∗ k φ ω (cid:48) k (cid:48) (cid:105) = (cid:104) φ ω ∗ k φ ω k (cid:105) δ ωω (cid:48) δ kk (cid:48) (Gizon & Birch 2004). Furthermore, we neglect termsin Eq. 10 that contribute when σ = 0, since we are inter-ested in finite σ . Therefore, the noise model is expressedas (cid:104)| B σ kq | (cid:105) = (cid:80) ω |H ωkk (cid:48) σ | (cid:104)| φ ω k | (cid:105)(cid:104)| φ ω + σ k (cid:48) | (cid:105) ( (cid:80) ω |H ωkk (cid:48) σ | ) , (11)where the mode power spectrum is given by (cid:104)| φ ω + σ k (cid:48) | (cid:105) = N k (cid:48) | R ω + σk (cid:48) | and (cid:104)| φ ω ∗ k | (cid:105) = N k | R ωk | . ANALYZING SUPERGRANULATION POWERSPECTRA
Table 1.
Mode coupling parameters n - n (cid:48) coupling | k | R (cid:12) range ω/ π range [mHz] (cid:15) f-f 400-1000 2.0-3.1 2p -p -p We compute power spectra of the B σ kq coeffcients us-ing n ≤ -p and p -p couplings. The firsttwo were chosen in order to compare to previous studies,while the cross n pair will demonstrate that supergran-ulation is also encoded in correlations of these modes.Table 1 shows the range of | k | and ω and the widths ofthe frequency windows used to compute the B σ kq coeff-cients. The power spectrum P ( σ, q ) for each data cubeis computed from Eq. 3 through P ( σ, q ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) k B σ kq (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (12)In this equation, we attribute equal weight to all k whencomputing P . Because all mode eigenfunctions are max-imally sensitive to the surface, this simple averagingscheme will result in the observational sensitivity peak-ing near the surface.Figure 2 shows slices of P ( σ, q ), averaged over the14 data cubes, computed from three different mode-coupling pairs. Cuts of the spectra at | q | R (cid:12) ≈ σ and angular distance ψ from the positive q x axis, show the so-called super-granular waves. These dual peaks in frequency can-not coalesce under a Galilean transform and have thusbeen interpreted as traveling waves (Gizon et al. 2003).Each of the mode-coupling spectra are consistent witheach other and follow the measured dispersion relationof Langfellner et al. (2018). At low temporal and spatialfrequencies, especially σ = 0, which we neglect, there isa systematic correlation noise that is equal to or greaterthan the signal.Slices of the spectra at σ/ π = 3 . µ Hz (middlepanels of Fig. 2) indicate that the f-f spectrum is con-sistent with the time-distance and LCT measurementsof Langfellner et al. (2018) that show enhanced powerin the prograde direction and peak power at | q | R (cid:12) ≈ -p spectrum is similar to the f-f spec-trum, though with a broader (approximately 20%) ringof power. The p -p spectrum differs considerably fromthe f-f and p -p spectra, with a significant depression nalyzing supergranular power spectra using helioseismic normal-mode coupling
180 90 0 90 180 [deg] / [ H z ] f-f coupling
180 90 0 90 180 [deg] p -p coupling
180 90 0 90 180 [deg] p -p coupling
200 0 200 q x R q y R
200 0 200 q x R q x R | q | R p o w e r x | q | R | q | R Figure 2.
From left to right are the results from f-f, p -p and p -p mode coupling, respectively. Top row: Slice of powerspectrum at | q | R (cid:12) ≈
172 as a function of frequency and azimuthal distance ψ from the positive q x axis. Overlaid is the dispersionrelation measured by Langfellner et al. (2018) (red line). Middle row: Slice of supergranular power spectrum at σ/ π = 3 . µ Hz.Bottom row: average power of supergranular spectrum (solid lines) and noise model (dashed lines). Three-sigma errors on themean are indicated by the shaded regions. in power for q R (cid:12) ≤ n - n (cid:48) coupling, where the smallest q measurable(largest flow scale) is limited by the difference in k be-tween p and p at similar frequencies. For example, at ω/ π = 3 mHz, wavenumbers associated with p and p differ by k R (cid:12) = 126, and hence those couplings are lesssensitive to larger scale flows ( | q | R (cid:12) ≤ n couplingsare inappropriate for supergranulation studies (e.g. f-p and p -p which are insensitive at 3 mHz to | q | R (cid:12) ≤ | q | R (cid:12) ≤ µ Hz ≤ σ/ π ≤ µ Hz) and the noise models are also shown inFig. 2 (bottom panels). These averages exclude | σ | ≤ . µ Hz which are dominated by systematic correlatednoise. The results show that the f-f couplings peak at q R (cid:12) = 120 and become consistent with the noise modelat q R (cid:12) ≈ -p couplings show a 60% greateramplitude than the f-f couplings. The peak in the p -p Hanson, Hanasoge and Sreenivasan σ / π [ μ H z ] |q|R ⊙ H W H M [ μ H z ] Figure 3.
Dispersion relation (top) and half-width at halfmaximum (HWHM, bottom) obtained from Lorentzian fitsof the mode-coupling supergranulation spectra. Fit parame-ters for the f-f spectra (black), p -p (red) and p -p (green)spectra are shown. For reference, we show the dispersionrelation of Gizon et al. (2003) (yellow curve) and the fits ofLangfellner et al. (2018) (blue diamonds), both computedfrom time-distance f-mode measurements. One-sigma errorsare shown for the coupling-spectra fits. spectrum occurs at q R (cid:12) = 145 and is consistent with thenoise model at q R (cid:12) = 300. Finally, the p -p couplingspeak at q R (cid:12) = 190, are consistent with noise at q R (cid:12) =300, and have amplitudes (and noise) ten times smallerthan the f-f spectrum.To further analyze these results we fit the spectra,for different | q | , using the parametric wave model ofLangfellner et al. (2018) (see their Eqs. 1 and 2). Inthe coupling spectra, there is correlated noise at small σ and | q | (e.g. mean flow fields), which introduces a biasin the fits. To reduce this bias and improve the stabil-ity of the fitting, we increase the frequency resolution ofthe Doppler cubes by a factor of ∼ .
5. This is achievedby taking each pair of data cubes in a Carrington rota-tion and stacking them in time, with seven days of zeropadding between the cubes. This reduces the numberof cubes in the averaged power spectra to seven, butconfines the noise to smaller bins near zero frequency,thus improving the fit. Figure 3 shows the dispersionrelation of the so-called supergranular waves, and thecorresponding half-width at half maximum (HWHM,1/lifetime). The fits from all three coupling spectra areconsistent in the range of 100 ≤ | q | R (cid:12) ≤ | q | R (cid:12) ≤
100 the fits begin to be biased by the lowfrequency noise. Below | q | R (cid:12) = 100 the p -p fits be-come unstable, while the p -p fits are unstable below | q | R (cid:12) = 75. Above | q | R (cid:12) = 225, the Lorentzian profilesbecome indistinguishable from diffuse power, and hence we do not fit beyond this scale. The measured HWHM ofall three spectra remain consistent with each other for all | q | R (cid:12) . Our measured dispersion relation and HWHMshow general agreement with the dispersion relation ofGizon et al. (2003) ( σ/ π = 1 . | q | R (cid:12) / . µ Hz)and the fits of Langfellner et al. (2018), both of whichare computed from time-distance f mode measurements. CONCLUSIONSWe have demonstrated the straightforward nature ofnormal-mode coupling for local helioseismic studies, andused it for detecting supergranules. The input datafor these calculations are the 3D Fourier-transformedtracked Doppler images, which also form the stan-dard input for other local helioseismic techniques.We have shown that supergranulation power spectracan be obtained by computing the cross-correlations φ ∗ ( ω, k ) φ ( ω (cid:48) , k (cid:48) ) of the wave field. We examined fand p self-couplings, finding that the resulting spectraare consistent with previous results (Gizon et al. 2003;Langfellner 2015). We also show that supergranulationmay be measured in cross couplings between p and p ,which are capable of probing greater depths than the for p modes.In all three coupling pairs, we identified the super-granular waves of Gizon et al. (2003), measuring thedistinct dispersion relation that is in general agreementwith that of previous studies. Our measurements of thelifetimes of these waves are also in agreement with pre-vious findings. The nature and origin of these super-granular waves is yet to be understood, but the resultsof this Paper show that there is information in vari-ous p mode coupling that could be used to shed furtherlight on this phenomenon. The agreement between ourmeasurements of the complex spatio-temporal dynamicsof supergranules and those of prior studies implies thatmode coupling will be an effective technique with whichto probe new and exciting questions.While we show the spectra for only a few n - n (cid:48) pairs, there is considerably more information to exploit.The mode-parameter fitting routine of the ring-diagrampipeline (whose output is used in Eq. 4) often fits upto p (for 15 ◦ tiles), enabling the use of 21 self- andcross-coupling pairs. The different depths and horizon-tal scales to which these couplings are sensitive provideinformation on the structure of supergranulation. Fu-ture studies that attempt to invert for the depth struc-ture of supergranulation will need to compute all of thesepairs. Furthermore, the information contained in thesecouplings could also be used to measure the Lorentzstresses (e.g. Hanasoge 2017) associated with supergran-ulation. nalyzing supergranular power spectra using helioseismic normal-mode coupling B σ kq coefficients requires careful analysis of sensitivity kernelsfor coupling (derived in Woodard 2006) and inversions atthe surface. This, and the yet-unknown depth profile ofsupergranules, may be addressed using mode-coupling. The authors downloaded HMI data from the GermanData Center at the Max Planck Institute for Solar Sys-tem Research. The Center for Space Science at NYUAbu Dhabi is funded by NYUAD Institute Grant G1502.The HMI data are courtesy of NASA/SDO and the HMIScience Team. Processing of HMI data was performedon the DALMA compute cluster at NYUAD.REFERENCEScoefficients requires careful analysis of sensitivity kernelsfor coupling (derived in Woodard 2006) and inversions atthe surface. This, and the yet-unknown depth profile ofsupergranules, may be addressed using mode-coupling. The authors downloaded HMI data from the GermanData Center at the Max Planck Institute for Solar Sys-tem Research. The Center for Space Science at NYUAbu Dhabi is funded by NYUAD Institute Grant G1502.The HMI data are courtesy of NASA/SDO and the HMIScience Team. Processing of HMI data was performedon the DALMA compute cluster at NYUAD.REFERENCES