Modeling Solar Oscillation Power Spectra. III. Spatiotemporal Spectra of Solar Granulation Velocity Field as Seen in HMI Velocity Measurements

We suggest a physically motivated model of the uncorrelated background, which can be used to improve the accuracy of helioseismic frequency measurements when the background contributes significantly to the formation of spectral lines of acoustic resonances. The basic assumption of our model is that the correlation length of the convective motions is small compared with the horizontal wavelength R ⊙/ℓ of the observations, where ℓ is the degree of the spherical harmonic Y ℓ m (θ, φ). When applied to solar power spectra at frequencies below acoustic resonances, the model reveals a distinct sensitivity to solar rotation: advection of the convective velocity pattern brings spatial correlations in the apparent stochastic velocity field (temporal correlations in the corotating frame induce spatial correlations in the inertial frame). The induced spatiotemporal correlations manifest themselves as an antisymmetric component in the dependence of the convective noise power on azimuthal order m, which allows us to address the solar differential rotation. With 360 days of data obtained by the Helioseismic and Magnetic Imager on board the Solar Dynamics Observatory, we measure three components of the rotation rate as a function of latitude using only ℓ = 300. This result indicates that the model suggests a new way of measuring solar subsurface rotation. This approach can complement traditional measurements based on correlation tracking.


Introduction
The most challenging task in contemporary helioseismology is to reduce systematic errors in estimating solar p-mode frequencies.This problem is most prominent when analyzing long time series, which are expected to have the smallest random errors.A large amount of data accumulated over the decades in dedicated ground-based and space projects calls for significant improvements to the data analysis pipeline to exploit their full diagnostic potential.For a recent account of the available data and its processing, we refer the reader to Larson & Schou (2015, 2018) and Korzennik (2005Korzennik ( , 2023)).
Multiple sources of potential systematic errors come into play when we attempt to measure an oscillation frequency with accuracy better than the width of its resonant line in the observed power spectrum.Systematic offsets are caused by inadequate modeling of the asymmetric line profile, inaccurate treatment of nearby spatial leaks, and incorrect magnitude and/ or frequency gradient of the uncorrelated background noise.Spatial leaks are characterized by the so-called leakage matrix and are subject to approximations regarding instrumental and optical distortions and mode-coupling effects.
This study is focused on global modeling of the uncorrelated background.It is common practice in most mode-fitting procedures to account for the uncorrelated background with a single free parameter for each (n, ℓ) frequency multiplet to ensure numerical stability.When dependence on azimuthal order m is allowed, it is evaluated in a small frequency interval in the vicinity of resonances, an interval typically contaminated by spatial leaks.We are seeking a description of the background in the entire range of (n, ℓ, m) by fitting a single slowly varying function of frequency only.
It is natural to examine noise power in a frequency range below 1 mHz, where global oscillation resonances are buried below the noise level.This study analyzes power spectra obtained from a 360 day long time series of Dopplergrams measured with the Helioseismic and Magnetic Imager (HMI) on board the Solar Dynamics Observatory (SDO).The time series we studied begins on 2019-3-14 and covers nearly a 1 yr period synced with HMIʼs regular 72 day global processing cadence and that is most closely centered on solar activity minimum.Figure 1 shows the observed power at degrees of ℓ = 300 and 100 as a function of m at frequencies around 300, 600, and 900 μHz; the measurements were averaged over ±100 μHz frequency intervals.
We can make two interesting observations: (i) For each of the two values of degree ℓ, the three curves obtained at frequencies that differ by a factor of 3 are essentially the same; the only difference is a nearly uniform vertical shift on the logarithmic scale.This feature indicates that the functional dependence of the noise power on the spatial spectral numbers (ℓ, m) and temporal frequency ω is separable, and (ii) the dependence on m is highly asymmetric.This feature points immediately to the effects of solar rotation, as the instrumentʼs sensitivity does not depend on the sign of m.
With our sign convention, harmonics with positive m are prograde waves, i.e., waves moving in the direction of rotation.In the corotating frame, these waves have smaller frequency and higher noise.
Below is our attempt to understand this behavior in detail.We assume that the noise comes from the turbulent convective Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
velocity field in the solar photosphere.In Section 2, we analyze the spectral measures of this noise, assuming that the correlation length of the convective motions is small compared with the observational wavelength R e /ℓ.We consider in detail the effects of differential rotation.Section 3 describes its measurement from the odd component with respect to m of the noise power in HMI measurements.We also analyze the even component with respect to m, governed by different sensitivity of the instrument to different spatial harmonics of the velocity field.Extension of the leakage matrix computations to include the instrumentʼs response to torsional components of the velocity field, which enter the analysis, is described in the Appendix.Section 4 suggests an initial approximation for the noise power in the frequency interval containing acoustic resonances, which has to be iteratively improved when fitting solar power spectra in frequency measurements.Our results are discussed in Section 5.

Spectral Measures of the Granulation Velocity Field in the Spatial and Temporal Domains
We work in a spherical coordinate system (r, θ, j) aligned with the solar rotation axis, and expand the time-dependent surface velocity field v(θ, j, t) in vector spherical harmonics as where ∇ 1 is the angular part of gradient operator, sin , and hats designate unit vectors.We make a Fourier transform of the time series of length T, with T large compared to the correlation timescale of the convective motions Using orthogonality properties of vector spherical harmonics where These expressions are obtained by taking the scalar product of both sides of Equation (1) with rY , , and r Y , , respectively, integrating over the angles θ and j and taking the Fourier transform.
We assume v(θ, j, t) to be a particular realization of a stationary stochastic process with zero mean.Quantities in the left-hand sides of Equations (4)-( 6), represented by stochastic integrals in the right-hand sides, are thus random variables with zero mean,  We associate v(θ, j, t) with the turbulent velocity field of convective motions imposed on a stationary large-scale background flow produced by differential rotation and meridional circulation.The basic assumption of our model is that for an observer moving with the background flow, the convective velocities do not correlate in space.For the granulation velocity field, this assumption can only be valid in a limited range of the harmonic degree ℓ, that is, where the typical size of a granule is small compared with the horizontal wavelength R e /ℓ.Here, we consider power spectra in the range of 0 ℓ 300 and limit our analysis to the background flow produced by differential rotation only: the effects of meridional circulation require a different treatment and will be left for further studies.
To make derivations more transparent, let us first consider a nonrotating model before generalizing the results to include the effects of solid-body rotation and then differential rotation.
(i) Nonrotating Sun.Changing the order of integration, we rewrite Equation (4) as , , ò q j w q j = w , and consider the covariance Considering the first integral on the right-hand side as a sum of integrals over small angular areas Δϖ i and likewise for the second integral indexed by j, we notice that the result is only nonzero for diagonal elements i = j; the expectation value of the cross-terms is zero because the v r ˜values are not correlated in space.We also know that the variance has an additive property.Therefore, if we replace the entire integration domain 4π on the right-hand side of Equation (8) by a small angular element Δϖ, the result will be proportional to Δϖ.We thus have where is positive spectral measure of the variance of vertical velocities, which we assume to be uniform over the solar surface.
We proceed in a similar manner with the contribution from the horizontal components of the velocity field, v θ (θ, j, t) and v j (θ, j, t), which have corresponding Fourier transforms v , , ˜( ) q j w q and v , , ˜( ) q j w j .We assume that horizontal velocities are isotropic in the azimuthal direction and, therefore, the two orthogonal horizontal components do not correlate with each other and where h 2 ( ) s w is the spectral measure of absolute values of horizontal velocities.The result is where we note that at ℓ = 0, the horizontal components are identically zero.
We expect no correlation between U, V, and W because of the orthogonality of corresponding velocity components and symmetry considerations.To show this, let v be the velocity vector at a particular point on the solar surface, with v v rY , . From geometrical considerations, the joint probability density function p(v r , v h ) is symmetric in v h , i.e., for any value of v r , two events with v h of the same magnitude but of opposite sign have the same probability.The expectation value of their product E[v r v h ] = 0, and hence, there is no correlation between U ℓm and V ℓ m ¢ ¢ .The same arguments apply to correlations between U and W and between V and W.
We now extend the analysis to include the effects of rotation.
(ii) Solid-body rotation.In this scenario, the Sun rotates with uniform angular velocity Ω in the observational frame.The effect of Coriolis forces on the velocity field at the scale of solar granulation is expected to be smeared away by spatial averaging, so we can assume that in the corotating frame, the observable statistical properties of convective motions are not affected by rotation.When the convective velocity field is observed in another reference frame, the only change is due to advection: in the spherical harmonic decomposition, a component of azimuthal order m will have its temporal frequency shifted by mΩ.The net result is that r 2 ( ) s w in Equation (9) has to be replaced with m r 2 ( ) s w -W , and similarly h 2 ( ) s w in Equation (12).(iii) Differential rotation.We now allow the rotation to change with latitude.When the rotation is uniform, the variance of U ℓm (ω) can be written as Since contributions to the variance coming from different latitudes simply add up, the same expression will be valid when Ω on the right-hand side is allowed to depend on latitude, meaning we can divide the spherical surface into thin latitudinal belts with each one in its own corotating frame.We will assume now that the rotation is slow, limiting the analysis to terms linear in Ω, and hence, where P z ℓ m ( ) are associated Legendre polynomials.Following an approach that is standard in solar seismology, we represent Ω(z) by an expansion where P s (z) are Legendre polynomials.Note that only even components of Ω(z) enter our result, as P z ℓ m 2 [ ( )] is an even function of z.
The required angular integrals are where s = 2k − 1, sℓ m g are odd polynomials of degree s in m defined in Ritzwoller & Lavely (1991), and  - are polynomials currently used in solar seismology to describe frequency splittings of solar oscillations, following normalization defined in Schou et al. (1994).Equation (17) can be derived by expanding dP s (z)/dz in P i (z), i < s, and evaluating integrals of triple products of Legendre polynomials.Convenient recurrence relations for evaluating  - can be found in Vorontsov (2007).We thus have Introducing a coefficients, commonly used in solar seismology, we have The relation between the expansion coefficients Ω s and a s is provided by Equation (17); in particular, Variances of V ℓm (ω) and W ℓm (ω) (Equations ( 5) and (6)) are transformed by the effects of differential rotation in precisely the same way.
We have an interesting observation: under the effects of differential rotation, each spectral component of velocity variances split in its observed frequency in precisely the same way as an undistorted frequency of solar oscillations would split under the effects of the same differential rotation if the influence of Coriolis forces can be discarded (leaving the effects of advection only) and differential rotation does not change with depth.
We also note that the possible inaccuracy introduced by linearization in the rotation rate (Equations ( 14) and ( 15)) can only affect the response to differential components.The response to the dominant Ω 1 component is treated correctly whatever its magnitude because

Solar Convective Velocity Field as Seen in SDO HMI
Power Spectra Instrumental response to different velocity field components does not depend on the sign of the azimuthal order m.Parameters of differential rotation (Equations ( 19) and (20)) can thus be calculated by shifting in frequency the power spectra of individual orders m to eliminate the odd (in m) component of the noise power.This procedure was implemented iteratively to account for a finite frequency window (±100 μHz in our measurement).The result obtained at ℓ = 300 at frequencies around 900 μHz is a 1 = 390.0± 0.9 nHz (synodic), a 3 = 20.0 ± 1.3 nHz, and a 5 = 2.4 ± 1.7 nHz.Corresponding coefficients of the polynomial expansion of the angular rotation rate (Equation ( 16)) are thus Ω 1 /(2π) = a 1 , Ω 3 /(2π) = −13.4± 0.9 nHz, and Ω 5 /(2π) = 1.3 ± 0.9 nHz.To evaluate the quality of this fit, we use a merit function defined as the RMS value of the residuals weighted by the observational uncertainties, which are evaluated under the standard assumption that observational power in an individual frequency channel has a χ 2 distribution with 2 degrees of freedom.Ideally, the value of this merit function should be close to 1.In our measurement, it is 1.073, which indicates that the targeted odd (in m) component is successfully eliminated in the derotated power spectra.
This result should be compared with other available measurements.Helioseismic measurements of solar internal rotation lose their accuracy in the subsurface layers, where the rotation varies rapidly with depth, and global modes lose their resolving power.However, the measurements reduced to solar activity minimum (Figure 8 of Vorontsov et al. 2002) indicate the surface values of Ω 1 /(2π) ; 435 nHz (sidereal, or synodic plus 31.6 nHz), Ω 3 /(2π) ; −13 nHz and Ω 5 /(2π) ; 1 nHz.The mean rotation rate Ω 1 /(2π) inferred from the convective noise is thus about 13 nHz slower; Ω 3 and Ω 5 appear to be in perfect agreement.
A classical result of measuring solar differential rotation using correlation tracking (Snodgrass & Ulrich 1990) sidereal, which translates to Ω 1 /(2π) = 468.0nHz sidereal, Ω 3 /(2π) = −5.15nHz and Ω 5 /(2π) = −1.46nHz.The difference between our measurement and this result is much bigger.One realistic scenario is that the measurements refer to different effective depths below the visible solar surface.Still, the result of Snodgrass & Ulrich (1990) is hard to reconcile with the results of our earlier helioseismic measurements (e.g., Vorontsov et al. 2002), where Ω 3 /(2π) ; −14 ± 1 nHz was found to be nearly constant with depth over the entire convective envelope, and Ω 1 /(2π) was found to increase with depth to a maximum value of about 449 nHz (sidereal) at a depth of about 6% of the solar radius.
Our measurement of the rotation of the solar granulation pattern requires a certain level of data quality.It benefits from going to a higher degree ℓ (wider range of azimuthal orders m), from observations with better spatial resolution (spatial leaks are not accounted for in the rotation measurement), and from observations of longer duration (better signal-to-noise ratio).When using 360 days of SDO HMI data, the measurement of relatively small differential components of the rotation rate loses stability at degree ℓ less than about 200, leaving the possibility of evaluating mean rotation only.This precludes the use of data from the Global Oscillation Network Group because they only provide spherical harmonic time series up to ℓ = 200.
Power spectra up to ℓ = 300 are indeed available from the Structure Program of the Michelson Doppler Imager (MDI) on board the Solar and Heliospheric Observatory (SOHO), but even these are unsuitable for measuring the differential components of the rotation due to contamination by spatial leaks resulting from insufficient spatial resolution and gaussian smoothing.Our detailed analysis of these power spectra in the entire range of degree 0 ℓ 300 shows that when ℓ increases, the spatial leaks induce a rapidly growing negative bias to the measurement of the mean rotation rate with a 1 starting from the degree ℓ of about 150.
We have attempted a measurement identical to that described above at ℓ = 300 but using 63 days of data taken by the Dynamics Program of SOHO MDI in 1996, which has a much higher resolution than the Structure Program data.The result is a 1 = 374.3± 2.2 nHz and a 3 = 20.1 ± 3.1 nHz, where a 5 could not be determined due to the shorter length of observation.While the a 3 coefficient agrees with the HMI measurement, the a 1 coefficient appears to be about 16 nHz smaller.We conjecture that the difference comes from contamination of the MDI power spectra with bigger spatial leaks due to a smaller spatial resolution of the instrument.Some of the difference may also be attributed to the different heights of formation of the spectral lines used by the two instruments: SOHO MDI was observing the Sun slightly higher in the atmosphere (Fleck et al. 2011).
With the odd (in m) component successfully eliminated in the properly de-rotated observational power spectra, we now analyze the remaining even component.For the same measurement at ℓ = 300 and frequencies of 900 ± 100 μHz, this component is shown by a thin line in Figure 2, where the remaining even component of the observed power of w , which is the m-averaged value of B ℓm 2 ( ) w .The contribution of the convective velocity field to the observational power spectra comes through multiple U ℓm , V ℓm, and W ℓm components (Equations ( 4)-( 6).As these components do not correlate with each other, we have, for the de-rotated power spectra, where we introduce the notation R, H, and T to designate separate leakage matrices which specify sensitivity coefficients of the instrument to radial components of the velocity field, horizontal components of the poloidal vector fields, and components of the toroidal fields (Equation ( 1)), respectively.
To make sure that a sufficient amount of spectral leaks are accounted for, the leakage matrices were computed with ℓ¢ in the range of ℓ ± 30 and m¢ in the range of m ± 30.Computation of the leakage matrices followed a semi-analytic approach described in Vorontsov & Jefferies (2005), which was generalized to include the instrumentʼs response to toroidal velocity fields; details can be found in the Appendix.To account for a finite resolution of the instrument in the CCD plane, the leakage matrix analysis involves the convolution of the images with a 2D Gaussian point-spread function (PSF).
When working with high-resolution HMI data in the intermediate range of degree ℓ 300, the width of the PSF was set to zero, meaning infinite resolution, or a PSF described by 2D Dirac δ-function.The solar B angle was set to 5.11°, the rms value of its annual variation (for a small B angle, its effect on the leakage matrix elements is quadratic in its magnitude).
The observed power B ℓm 2 ( ) w considered as a function of m at ℓ = 300, is fitted by a linear combination of the two terms in the right-hand side of Equation ( 23 s s =  indicates that horizontal velocities in the turbulent flow are about three times bigger than vertical velocities.The fit quality remains adequate when the same analysis is applied to data at a smaller degree of ℓ.An interesting observation is that the measured ratio h r 2 2 s s increases monotonically to 18.2 ± 0.5 at ℓ = 100 and 20.5 ± 2.9 at ℓ = 5.At a degree of ℓ < 5, this measurement loses stability due to an insufficient number of the available orders m.This finding may indicate that bigger convective cells have a bigger average ratio of horizontal to vertical velocities. Another finding is that the observed m-averaged value B 2 ( ) w stays nearly constant in the entire range of degrees: it drops monotonically when going from ℓ = 0 to 300, but only by about 15%.For comparison, in medium-ℓ SOHO MDI measurements of much smaller spatial resolution, this variation amounts to 2 orders of magnitude.This behavior indicates that in the range of degrees of ℓ 300, the spatial resolution of the HMI instrument is indeed almost perfect.To clarify this point, our analysis can be made independent of the leakage matrix computations-assuming, of course, that the spatial resolution of the instrument is perfect.
In addition to the coordinate system (θ, j) with the z-axis aligned with the solar rotation axis, we introduce another coordinate system , ( ) q j ¢ ¢ with the z¢-axis (from which q¢ is counted) directed from the Sun toward the observer.Considering the projection of the turbulent velocity field v(θ, j, t) (Equation ( 1)) on the CCD plane directly, without its decomposition in vector spherical harmonics, we have a result that does not depend on the target degree ℓ.By expanding cos sin 2 2 ( ) q q ¢P ¢ and sin sin 2 2 ( ) q q ¢P ¢ in Equation (25) in spherical harmonics and transforming the result to (θ, j)-coordinates, it is also possible to evaluate the right-hand side at individual m-values.We skip the details of this analysis, as its principal motivation was to check the accuracy of our leakage matrix computations.The numerical results of the two approaches turned out to be the same.
The slight variation of the apparent values of B 2 ( ) w with degree ℓ indicates that the leakage matrices can be improved by setting the PSF width to a small but nonzero value.We conclude that the measurements of the solar noise can be used to calibrate the effective PSF of the instrument.This option may be particularly interesting for analyzing data obtained with the SOHO MDI instrument.

Temporal Domain
We can de-rotate the power spectra by frequency shifting them according to the results of the differential rotation measurement.The m-average of the de-rotated spectra in the entire frequency domain of SDO HMI data at ℓ = 300 is shown in Figure 3.At frequencies less than about 200 μHz, the variation of the observed power with ℓ and m cannot be explained by our model, which loses its ability to fit the data with any reasonable accuracy.In this spatiotemporal domain, our assumption of negligibly small correlation length is violated by supergranularscale convective motions, as seen in Doppler-velocity power as a function of m and ω, shown in a grayscale in Figure 4.The well-defined ridge at m > 0 (prograde waves) is produced by the rotation of the solar supergranulation pattern.A small but noticeable curvature of the ridge is due to faster rotation of the equatorial regions.At frequencies less than about 50 μHz, the observed power drops rapidly because of the de-trending applied to the time series of spherical harmonic components.
In this study, the data analysis was limited to frequencies below the oscillation resonances.We can hope that the dependence of the noise power on ℓ and m (the B m ℓ 2 ( )) measured in this frequency range will stay the same at higher frequencies; this assumption, of course, remains to be verified by addressing residuals of spectral fitting procedures.We note that at frequencies from about 3 mHz and higher, the measurement of the background component is difficult because the signal-to-noise ratio of acoustic resonances becomes very high.At these frequencies, it is now the uncorrelated background that gets buried below the resonant power.
We suggest a simple model for the frequency dependence of the background noise B 2 ( ) w to be used as an initial guess in the mode-fitting procedures.Imagine a convective eddy emerging on the solar surface from below at time t = 0. Let the observed   We adjust a linear combination of these seismic events with two different values of τ to approximate the expected variation of the uncorrelated background in the entire frequency range; the result is shown in Figure 3 by two dashed lines for the two separate components and by the thick gray line for their sum.The fitted values of τ, about 6 minutes and 1 minute, are of the order of the lifetimes of solar granules and shorter.A relative excess of observational power at the highest frequencies may be due to an aliasing signal coming from frequencies higher than the Nyquist frequency.

Discussion
Implementation of our model in frequency measurements is relatively straightforward.At each degree, ℓ, the even functions B m ℓ 2 ( ) are measured from the de-rotated power spectra around some frequency below all the detectable resonances.The initial approximation for the frequency dependence of the uncorrelated background B 2 ( ) w is then improved by fitting individual multiplets in the power spectra.When the instrumentʼs resolution is imperfect, as with SOHO MDI measurements or those of a high degree ℓ from SDO HMI, the minor sensitivity to modes of higher degree ℓ will be captured in B m ℓ 2 ( ).The suggested measurement of differential rotation from power spectra at frequencies below acoustic resonances should be undertaken with more data of lower and higher degree ℓ analyzed in different frequency intervals.It would be interesting to extend these measurements to data sets obtained at different times to explore temporal variations of the subsurface differential rotation, also known as torsional oscillations.
In our limited exercise with observational data, we have another finding that deserves more extensive data analysis.The inferred ratio of magnitudes of horizontal and vertical components of convective velocities h r 2 2 s s clearly tends to get bigger when ℓ gets smaller; it indicates that in bigger convective cells, horizontal velocities become more dominant.
Our model becomes inconsistent with observations at frequencies below 200 μHz, since our basic assumption of small correlation length breaks down when observations are affected by signals from supergranular convective cells.Here, we enter the spatiotemporal domain targeted by Beck & Schou (2000) in their encouraging measurements of differential rotation of solar supergranulation pattern from Dopplergrams provided by the SOHO MDI instrument.An approach that is more sophisticated than ours is needed to deal with turbulent velocity correlations simultaneously in both space and time.
result is shown in Figure2by a thick gray line.Visual inspection of the fit quality and the value of the corresponding merit function indicates that the approximation of the measured function of m by a linear combination of two functions coming exclusively from the leakage matrix analysis is perfectly adequate.The inferred ratio 9

Figure 2 .
Figure2.Even component in the observational noise power as a function of azimuthal order m at l = 300 measured at frequencies around 900 μHz (thin line).The thick gray line shows its approximation obtained by fitting synthetic power.
is the Fourier transform of v(θ, j, t) at frequency shifted by advection effects, and sin ( ) q P ¢ is an apodization function, sin q¢ being the radial coordinate in the image plane in units of the apparent solar radius.Evaluating the measure of the stochastic signal in the way described in Section 2 line-of-sight projection effects.Using an addition theorem for spherical harmonics, the m-averaged value is

Figure 3 .
Figure 3. De-rotated and m-averaged SDO HMI power spectrum at l = 300 (thin line).Dashed lines show two simple models (see text); their sum is shown by the thick gray line.