Issue 
A&A
Volume 633, January 2020



Article Number  A109  
Number of page(s)  16  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201936662  
Published online  20 January 2020 
An improved multiridge fitting method for ringdiagram helioseismic analysis
^{1}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: nagashima@mps.mpg.de
^{2}
JILA & Department of Applied Mathematics, University of Colorado, Boulder, CO 803090440, USA
^{3}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, FriedrichHundPlatz 1, 37077 Göttingen, Germany
Received:
10
September
2019
Accepted:
6
November
2019
Context. There is a wide discrepancy in current estimates of the strength of convection flows in the solar interior obtained using different helioseismic methods applied to observations from the Helioseismic and Magnetic Imager onboard the Solar Dynamics Observatory. The cause for these disparities is not known.
Aims. As one step in the effort to resolve this discrepancy, we aim to characterize the multiridge fitting code for ringdiagram helioseismic analysis that is used to obtain flow estimates from local power spectra of solar oscillations.
Methods. We updated the multiridge fitting code developed by Greer et al. (2014, Sol. Phys., 289, 2823) to solve several problems we identified through our inspection of the code. In particular, we changed the (1) merit function to account for the smoothing of the power spectra, (2) model for the power spectrum, and (3) noise estimates. We used Monte Carlo simulations to generate synthetic data and to characterize the noise and bias of the updated code by fitting these synthetic data.
Results. The bias in the output fit parameters, apart from the parameter describing the amplitude of the pmode resonances in the power spectrum, is below what can be measured from the MonteCarlo simulations. The amplitude parameters are underestimated; this is a consequence of choosing to fit the logarithm of the averaged power. We defer fixing this problem as it is well understood and not significant for measuring flows in the solar interior. The scatter in the fit parameters from the MonteCarlo simulations is wellmodeled by the formal error estimates from the code.
Conclusions. We document and demonstrate a reliable multiridge fitting method for ringdiagram analysis. The differences between the updated fitting results and the original results are less than one order of magnitude and therefore we suspect that the changes will not eliminate the aforementioned ordersofmagnitude discrepancy in the amplitude of convective flows in the solar interior.
Key words: Sun: helioseismology / methods: data analysis
© K. Nagashima et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1. Introduction
Understanding solar interior dynamics is crucial to understanding the mechanisms of the solar dynamo. As one example, convection may play an important role in the formation of magnetic flux tubes, as well as in their rise through the convection zone and their tilts at the solar surface (e.g., Brun & Browning 2017). Helioseismology, which uses observation of oscillations on the solar surface, is an important probe of interior dynamics.
Currently there is a major discrepancy between the timedistance (Hanasoge et al. 2012) and ringdiagram (Greer et al. 2015) estimates of the strength of solar subsurface convection at large spatial scales. The measurements from timedistance helioseismology suggest flows ordersofmagnitude weaker than those seen in convection simulations (e.g., in the anelastic spherical harmonic (ASH) convection simulations by Miesch et al. 2008), while the measurements from ringdiagrams are closer to the expectations from simulations.
Timedistance helioseismology (Duvall et al. 1993; Kosovichev & Duvall 1997) is based on measuring and interpreting the travel times of acoustic and surfacegravity wave packets. These travel times are measured from the temporal crosscovariances between the Doppler observations at pairs of points on the solar surface (see Gizon & Birch 2005, for a review). Ringdiagram analysis (Hill 1988; Antia & Basu 2007) measures the Doppler shift of acoustic and surfacegravity oscillation modes in the local power spectra and uses these Doppler shifts to infer local flows in the solar interior. To compute the local power spectra, the solar surface is divided into a number of spatial tiles. Each tile is tracked at a rate close to the solar rotation rate. The power spectra of the solar oscillations (Doppler observations) are computed for each tile. In the threedimensional power spectra, roughly concentric rings with high power are present at each frequency and they correspond to the modes of different radial orders. Flow and wavespeed anomalies in the Sun shift and distort the rings and, hence, one can obtain information about the solar interior from the ring parameters. Flow maps from Doppler observations by the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) onboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) are automatically computed by the SDO/HMI ringdiagram pipeline (Bogart et al. 2011a,b) on a daily basis since the SDO launch in 2010.
The HMI ring pipeline codes separately fit each single ridge (single radial order n) in the power in slices at constant horizontal wavenumber or slices in temporal frequency. Greer et al. (2014) developed an alternative approach based on simultaneously fitting multiple ridges (multiple radial orders) at each horizontal wavenumber. Greer et al. (2015) introduced another innovation: they chose a much denser tile layout than other ringdiagram analyses: 16degree tiles with the separation of 0.25 degree instead of the 7.5 degree spacing (at the equator, the spacing increases at higher latitude to maintain 50% overlap between neighboring tiles) used for 15degree tiles in the HMI ring pipeline. As described in detail in Appendix A.2, the ATLAS code from Greer et al. (2014, 2015) provides seven parameters for each ridge and five parameters for the background power at each wavenumber k. Greer et al. (2015) applied a threedimensional flow inversion to these fit results to estimate the threedimensional flow field in the solar interior. Both the dense packing of tiles and the threedimensional inversions are unique to ATLAS, however, in this paper we focus only on the fitting component of the code.
As one step toward understanding the causes of the abovementioned disagreement between the helioseismic measurements of subsurface convection, here we focus on the ringdiagram analysis described by Greer et al. (2015). We revisit the analysis code (Greer et al. 2014; Greer 2015) that was used in that work and identify several issues through a stepbystep examination of the code. In response to these issues, we have developed an updated method and characterize the updated code by applying it to synthetic data generated from MonteCarlo simulations.
2. Description of the updated code
In this section, we describe our updated ATLAS code. Each of the updates is a response to a problem that we found in our inspection of the original code. In this section we will refer to Appendix A for the details of the original code.
After computing the power spectrum for a particular tile, the processing steps are: (1) remap the power spectrum from Cartesian (k_{x}, k_{y}) to polar (k, θ) coordinates, (2) rebin the power in azimuth θ; the number of grid points in θ is reduced from n_{pix} = 256 to n_{pix} = 64, (3) fit the logarithm of a Lorentzian model to the logarithm of the smoothed power by leastsquares minimization at each k using the Levenberg–Marquardt (Marquardt 1963) technique; the model function has 7n_{r} + 4 parameters at each horizontal wavenumber k, where n_{r} is the number of ridges at the particular value of k, and (4) estimate the covariance matrix of the errors of the fitted parameters by computing the inverse of the Hessian matrix of the cost function. In the following subsections we describe the changes that we introduced to each of these steps.
2.1. Rebinning in azimuth
The computed local power spectra O(k_{x}, k_{y}, ν) and the interpolated spectra in polar coordinates are nonnegative by construction. Following the original code, the interpolated spectra have 256 pixels at each k. At k R_{⊙} ≡ ℓ ∼ 500 (where R_{⊙} is the solar radius), which corresponds to k_{pix} ≡ k/h_{k} = 21 (where h_{k} = 3.37 × 10^{−2} Mm^{−1} is the grid spacing in k) and which we use in most of the MonteCarlo test calculations shown later, this is about twice of the number of grid points in θ from the full resolution, which has n_{pix} ≈ 2πk/h_{k}. For the sake of computational efficiency, it is desirable to reduce the number of points in azimuth. In the updated code, we use a running boxcar smoothing of fourpixel width followed by subsampling by a factor of four to reduce to the number of grid points in θ. This procedure ensures that the resulting smoothed power spectrum O_{s}(k, θ, ν) will be positive. We expect that as flows produce θ variations that are dominantly at azimuthal wavenumber one, it should be possible retain only very low resolution in θ; Sect. 4.1 discusses this in more detail. We retain 64 points in θ for the examples shown in this paper.
The original code used a lowpass Fourier filter to smooth the power spectrum in the θ direction. This procedure produces occasional points where the smoothed power is negative.
2.2. Leastsquares fitting of the logarithm of power
In the updated code, we use a leastsquares fitting to fit the logarithm of the model power to the logarithm of the remapped and smoothed power. Following Greer et al. (2014), the fitting is carried out independently at each horizontal wavenumber k. As the amount of smoothing is increased, the probability distribution function (PDF) of the power, as well as its logarithm, approaches a normal distribution (see Appendix B.4 for more details), and it is therefore appropriate to use a leastsquares fit. The logarithm of the smoothed power has the convenient property that the variance of the logarithm of power (σ_{N} in Eq. (B.4)) depends only on the details of the remapping and smoothing and does not depend on θ or ν (see Appendix B.2 for details).
In our approach, the cost function at a single k is:
where O_{s}(θ, ν) is the observed spectrum at some fixed k after smoothing in θ and P(θ, ν; q) is the model of the spectrum with model parameters q. The summation is taken over all bins j within the fitting range. We note that as the error estimates for lnO_{s} (σ_{N} in Eq. (B.4)) are all the same, we set them all to one in writing the cost function. For the sake of readability, throughout the remainder of this paper we do not introduce or carry notation to denote the value of k; the fitting problems at each k are treated as independent and we will not be comparing fit parameters for different values of k. This is an approximation as the interpolation from (k_{x}, k_{y}) space to (k, θ) space does imply error correlation between the fit parameters at different values of k.
In the updated code, we use the Levenberg–Marquardt technique to solve the minimization problem. In particular, we use mpfit.c (Markwardt 2009), one of the codes in the MINPACK1 leastsquares fitting library (Moré 1978; Moré & Wright 1993). The covariance matrix of the errors associated with the fitted parameters are estimated at the last step of the LevenbergMarquardt procedure. As a practical note, the error estimates obtained by this method must be scaled as we have assumed σ = 1 in Eq. (1); see original PDF, Eq. (B.4). In general, calling the code with an incorrect estimate of σ could cause poor performance of the fitting algorithm. In the current case σ is not far from one (σ ∼ 0.5, see Sect. 3.1) and we do not expect that this is a significant issue here. To implement the σ estimation in the code is a task for the future.
The original ATLAS code used a maximumlikelihood method based on the assumption that the power spectrum in a single bin in (k_{x}, k_{y}, ν) space follows a chisquared distribution with two degrees of freedom; see Appendix B.1 for the original likelihood function. This method does not account for the impact of smoothing on the PDF of the power spectrum.
2.3. The functional form of the power spectrum
The model function for the power spectrum at a single k in our updated code is
where
The fitting parameters for the peaks (n = 0, …n_{r} − 1, where n_{r} is the number of the ridges) are q_{0, n} = ν_{n} is the frequency of the nth peak, q_{1, n} = A_{n} is the amplitude, q_{2, n} = Γ_{n} is the width, (q_{3,n}, q_{4,n}) = u = (u_{x,n}, u_{y,n}) is the horizontal velocity, q_{5, n} = f_{c, n} and q_{6, n} = f_{s, n} are parameters to handle anisotropy. The background is modeled at each k as
where F is again given by Eq. (3) and q_{0, BG} = B_{0} is the amplitude, q_{1, BG} = b is the powerlaw index, q_{2, BG} = f_{c, bg} and q_{3, BG} = f_{s, bg} are the parameters to handle the anisotropy. The number of the parameters in total is 7n_{r} + 4 at each k. For reference, Tables C.1–C.3 show the physical meaning of each of the fitting parameters.
We altered the model function Eq. (2) from the original Eq. (A.1) to make the parameterization more stable. The following subsections describe the motivation for these changes.
2.3.1. Parameters of the anisotropy terms
We replaced F′ defined by Eq. (A.2) with F defined by Eq. (3) and used the same form for the anisotropy terms for the background, although the exact form for the background function is subject to other alterations discussed in Sect. 2.3.2.
Our alteration does not change the space of functions that can be fit with F(θ) but it does not suffer from the issue of indeterminate phase for nearly isotropic power spectra. In the original parameterization, the amplitudes of the anisotropic part of the model function ( and ) are usually much smaller than one. As a consequence, the phase ( or ) does not matter much in the fitting, and, hence, it is unstable. In the particular case of isotropic power spectra with for all n and , the phases and are indeterminate.
2.3.2. Parameters of the background model
We also changed the background parameterizations for the sake of stability. The background model (Eq. (A.3)) originally contained five parameters. In this section we explain our motivation for reducing this to the four parameters shown in Eq. (4).
The background model in Eq. (A.3) is based on the model of Harvey (1985):
where τ is the characteristic timescale of the velocity field in question and σ_{rms} is the rms velocity. In this case the index in the original background model (Eq. (A.3)) would be . Appourchaux et al. (2002) suggested a generalization of the Harvey (1985) model, where can be the range of 2 to 6.
In the ATLAS fittings, however, we typically obtain the index for HMI observations. Also, the rolloff frequency obtained from the fitting of HMI observations, , is quite low (∼1 μHz) and below the frequency resolution for 28.8 h power spectra (9.7 μHz), which is the typical observation length for 15degree tiles in HMI pipeline and 16degree tiles in ATLAS.
This suggests that the background is not related to either the supergranulation (τ ≈ 10^{5} s, or ν ≈ 10 μHz) or granulation time scales (τ ≈ 4 × 10^{2} s, ν ≈ 2.5 mHz), based on Table 1 of Harvey (1985). As these background parameters are not consistent with the original physical model, we might need to reconsider the background model. At this moment, however, this is a task for the future, and we retain this model in the altered form mentioned below.
In the case of , the original background model, Eq. (A.3), can be simplified:
we therefore redefine the background model B(θ, ν) with four parameters q_{i, BG}(i = 0, 1, 2, 3) with the altered anisotropy terms in the form of Eq. (4). As the new background model (q_{i, BG}) has only four parameters instead of the original five, the index i has different meanings in the two models.
3. Performance of the updated code
MonteCarlo simulations are a powerful tool for testing fitting methods. In this section, we use this approach to characterize the performance of the updated fitting code. In Sect. 3.1 we compare the scatter of the fitting results of MonteCarlo simulations with the noise estimated by the updated code. In Sect. 3.2 we measure the bias of the flow estimates for some simple cases.
3.1. Error estimates
We use the approach of Gizon & Birch (2004) to generate realizations of the wavefield. The assumption of this approach is that the real and imaginary parts of the wavefield at each point (k_{x}, k_{y}, ν) are independent Gaussian random variables. In more detail, the method is as follows: (1) create a Lorentzian model (Eq. (2)) for the limit spectrum using a set of input parameters, (2) pick two standard normally distributed random numbers at each grid point (k_{x}, k_{y}, ν) and take sum of the squared numbers divided by two to make a chisquare distribution with two degrees of freedom, and (3) multiply the limit spectrum by these random numbers to obtain one realization of the power spectrum.
After each realization of the power spectrum is generated in (k_{x}, k_{y}, ν) space, it is then remapped and rebinned in the same manner as the local power spectra computed from the observations. To obtain the parameters for the input model power spectrum in the first step above, we used the average over Carrington rotation 2211 of the power spectra for the disc center tile from the SDO/HMI pipeline (these average power spectra were obtained using the Data Record Management System (DRMS) specification hmi.rdvavgpspec_fd15[2211][0][0]). We then used the updated code to fit these average power spectra. The parameters resulting from these fits at ℓ = 328, 492, and 984 (k_{pix} = 14, 21, and 42, respectively) are shown in Tables C.1–C.3, respectively. We later use these parameters to generate MonteCarlo synthetic data.
Figure 1 shows a few slices of the input limit spectrum, namely, Eq. (2) evaluated for the parameters q listed in Tables C.1–C.3. The original observed power spectrum is also shown. Figure 1 shows that the model power spectrum is reasonable compared to the observables, although the details of the peak shapes and the background slopes are not always reproduced by the model and there is still room for improvement with regard to the model function.
Fig. 1.
Slices through the input limit spectrum (black) and observational spectrum averaged over one Carrington rotation (gray), which was used to obtain the input parameters, at θ = 0 and a few different values of k. Vertical dashed lines indicate the fitting ranges and the models are plotted only within these ranges. Lower limits for all k are fixed at 0.4 mHz, while the upper limits depend on the initial guess for each k, and they are the highest peak frequencies plus the widths of the peak. For the cases of ℓ = 328 and 492 (k_{pix} = 14 and 21, respectively), the peaks with the radial order, n, from 0 to 7 are used in the fitting, while for the case of ℓ = 984 (k_{pix} = 42), 0 ≤ n ≤ 3. 
We created 500 realizations of the power spectrum using the above procedure. The fitting code occasionally produces outliers and for these computations, we removed them. Specifically, we removed outliers iteratively: in each iteration we computed the standard deviation of the samples between the tenth and 90th percentile and we discarded points at more than five times the standard deviation away from the mean. We repeated this procedure until no further points were removed. After this outlier removal, the number of valid samples are 489 (ℓ = 328, k_{pix} = 14), 500 (ℓ = 492, k_{pix} = 21), and 496 (ℓ = 984, k_{pix} = 42) out of 500 samples.
Figure 2 shows the average and scatter of the output peak parameters associated with the n = 0…5 ridges for ℓ = 492 (k_{pix} = 21) from fitting 500 MonteCarlo realizations of the power spectrum. Table 1 shows the average and scatter of the output background parameters from the same set of fitting results. From Fig. 2 and Table 1, we see that most of the averages of the parameters estimated by the fitting are within the expected scatter in the mean (, where N_{sample} is the number of MonteCarlo realizations, hence N_{sample} = 500 here). The amplitude is an exception; it is always smaller than the input. This is the result of taking the logarithm of the smoothed power; see Appendix D for details. Also, from Fig. 2 and Table 1, we see that the error estimated by the updated code σ^{code} is consistent with the scatter of the samples σ^{scat}. This shows that the error estimates produced by the updated code are reasonable.
Fig. 2.
Fitting parameters and their errors for the case ℓ = 492 (k_{pix} = 21). Each panel has three parts. In the upper part of each panel, squares indicate the input parameters, . The asterisks indicate the fitting results obtained by the updated code, namely, the mean of MonteCarlo samples, . In the middle part of each panel, the deviation of the mean fitting results from the input, , are illustrated by the circles and the expected scatter of the mean computed from the square root of variance σ^{scat} of the 500 Monte Carlo samples divided by the square root of the sample number, N_{sample} (here it is 500) are shown as the error bars. The dashed horizontal lines are at δq = 0. In the lower part of each panel the scatter of the 500 samples, σ^{scat} and the scaled error estimated by the updated code, σ^{code}, are depicted by the thick gray lines and the thin lines with short horizontal bars on the edge. Errors estimated by the updated code are scaled by σ(α = 1.5, k_{pix}, n_{pix}) = 0.569 as described in Appendix B.2. The error bars on the middle panel of δq_{1, n} are tiny at this scale; the underestimation is relatively large. However, the errors of δq_{1, n} are and, therefore, available from σ^{scat} in the lower panel and N_{sample} = 500; they are ∼0.02 at most. 
Background parameters in the input model for the Monte Carlo simulation and the fitting results at ℓ = 492 (k_{pix} = 21).
We carried out the analogous MonteCarlo simulations for the cases of twice larger k and twothirds k. Although the number of peaks to be fitted is not the same in these cases as the number of prominent peaks is smaller for larger k, the behavior of the fitting results and their error estimates are similar to those in the case shown in Fig. 2 and Table 1.
3.2. Bias of the flow estimates and correlations between the fitting parameters
To measure the bias of the flow estimates and the correlations between the fitting parameters, we make models with a simple flow: isotropic models, except for a specific flow u_{x} for the n = 3 ridge.
The parameters for the model are identical to those from Sect. 3.1, except for the parameters related to the azimuthal angle θ; namely, the peak parameters q_{3, n}, q_{4, n}, q_{5, n}, and q_{6, n} (for all n) and the background parameters q_{2, b} and q_{3, b} are all set to zero. The flow u_{x, n} (q_{3, n}) is nonzero for a single n, n = 3. The limit spectrum is constructed with Eq. (2) and these input parameters, and the realizations of the power spectrum are generated in the same way as those in Sect. 3.1.
In this subsection we compare results from the updated code with results from the original code and also a modified version of the original code. While the original code fits a Lorentzian model to the square root of the power, the modified original code fits a Lorentzian model to the power itself. Although the original and modified original codes assumes the fiveparameter background model (Eq. (A.3)), here we created the limit spectrum using the fourparameter background model (Eq. (4)), which we use in our updated code. As described in Sect. 2.3.2, under the current conditions, our updated fourparameter background model can approximate well the original fiveparameter background model and the original and the modified original codes are applicable, although we need to keep in mind that the meaning of the background parameters in the results from different codes are different. Comparable tests using input parameters obtained by the original code give similar correlation coefficients maps.
3.2.1. Bias of the flow estimate
Figure 3 shows the fitting results for u_{x, n = 3} at ℓ = 492 (k_{pix} = 21) in the form of normalized histograms of 500 Monte–Carlo simulations. The input models have q_{3, n} = u_{x, n} = 0, 80, …400 m s^{−1} for the n = 3 ridge. Before plotting, we removed outliers from the output fit parameters. We removed outliers as described in Sect. 3.1. After the outlier removal procedure, the sample numbers for each methods are 394–466 depending on the value of u_{x, n = 3} (79–93%, modified original), 436–488 (87–98%, original), and 499–500 (100%, updated). The updated code provides almost no outliers, while there are some outliers in the fitting results by the modified original and original codes. The number of outliers depends on u_{x, n}, although there is no clear trend with u_{x, n}; for example, we cannot say that the stronger flow produces more outliers. We did not further investigate the details of the outliers in the fitting results by the original and modified original codes.
Fig. 3.
Fitting results of u_{x, n = 3} at ℓ = 492 (k_{pix} = 21) in the form of normalized histograms by the original code which fit the model to the square root of power (black), the modified original code which fit the power (blue), and the updated code (red). Horizontal lines indicate the average of 500 Monte–Carlo realizations given by the three codes (in the same colors as the histograms, dashdotted black, dashed blue, and solid red). The short vertical lines on the horizontal lines indicates the error estimated by the codes (thin lines with short horizontal bars) and the standard deviations of the 500 fitting results (thick lines) in the same color as the histograms centered at the means (horizontal lines) by the three codes. We note that the fitting results with overly large deviation are omitted. See text for details. 
Figure 3 shows that the original fitting code underestimates the input flow by about 3%. This trend is consistent with what was reported by Greer et al. (2014). The modified original and updated codes produce lessbiased flow estimates. Figure 3 also compares the errors estimated by the code and the scatter in the Monte Carlo simulation. The errors from the updated code are consistent with the scatter of the MonteCarlo results, while the original code overestimates the errors.
3.2.2. Correlations between the fitting parameters
Figure 4 shows the correlation coefficients between the fitting parameters of 500 Monte Carlo samples. For this computation, outliers were removed as described in Sect. 3.1. The original and modified original codes both produce stronger correlations in some of the output parameters than the updated code. In particular, for the original code, these parameters show the strongest correlations: the amplitude A_{n} () and the width Γ_{n′} () at peaks n − n′ ≤ 1, the rolloff frequency ν_{bg} () and the index b () of the background, and the index b () of the background and A_{n} () or Γ_{n} () of highern peaks.
Fig. 4.
Correlation coefficients between the parameters (500 MonteCarlo realizations) for the isotropic model (u_{x, n} = 0 for all n) and for the isotropic model except with u_{x} = 400 m s^{−1} for the n = 3 peak (on the next page) at ℓ = 492 (k_{pix} = 21). Each box indicates the parameters for nth peak (n = 0, 1, …5), q_{i, n} (i = 0, …6, from left to right and from bottom to top in each box) and the background parameters (BG). As we defined in Sect. 2q_{0, n} = ν_{n} is the frequency of the nth peak, q_{1, n} = A_{n} is the amplitude, q_{2, n} = Γ_{n} is the width, q_{4,n}) = u = (u_{x,n}, u_{y,n}) is the horizontal velocity, q_{5, n} = f_{c, n} and q_{6, n} = f_{s, n} are parameters to handle anisotropy. Background parameters are q_{0, BG} = B_{0} is the amplitude, q_{1, BG} = b is the powerlaw index, q_{2, BG} = f_{c, bg} and q_{3, BG} = f_{s, bg} are the parameters to handle the anisotropy of the background. Color scale is shown by the color bar on the right side. We note that the last columns and rows (i = 4) of BG on the top panels are empty because the number of the background parameters in the model in the updated code is four, instead of five. 
Fig. 4.
continued. 
The modified original code and the updated code did not show such strong biases, except for the width and amplitude of each peak, along with the background index and some weaker peaks (smaller and larger n). In terms of correlations between the parameters, the updated code shows an overall improvement in comparison with the original, but the correlation coefficients between u_{x} and u_{y} on the same peaks or on the peaks next to each other are ∼0.1 at most for any u_{x} and any fitting results by all three codes shown in Fig. 4.
3.3. Summary of the MonteCarlo test
The MonteCarlo tests show that the updated code is able to reasonably recover the parameters, apart from the amplitude, that are used to generate the input power spectra. Other than the amplitudes, we were not able to measure a statistically meaningful bias in the fit results using 500 MonteCarlo simulations. The underestimation of the amplitude is the result of taking the logarithm of the rebinned power. Appendix D discusses this issue in detail. Since current ringdiagram flow inversions are carried out using the mode shifts, (q_{3, n}, q_{4, n}) = (u_{x, n}, u_{y, n}), and the other fitting parameters including mode amplitudes are not used, we believe underestimated amplitudes will not substantially impact on further analysis.
4. Discussion
4.1. Further rebinning in azimuth
We explored the idea of rebinning the data further in θ. If we rebin further, the PDF of the resulting power spectrum is closer to Gaussian. Additional rebinning has the additional benefit of reducing the underestimation of the amplitude. Another benefit is that if we can reduce the number of pixels involved in the fit without a significant decrease of the fitting quality, it will reduce calculation costs.
Figure 5 shows the fitting results corresponding to the fourtime further rebinning. Specifically, we rebin the 64pixel azimuth grid into 16 pixels. In this case, valid sample numbers for the three codes are 351–429 (70–86%, modified original), 406–487 (81–97%, original), and 500 (100%, updated). While the updated code had no outliers, the numbers of outliers in the results by the original and modified original codes increase compared to the case without further rebinning (see Sect. 3.2). This also confirms the stability of the updated code.
Fig. 5.
Fitting results of u_{x, n = 3} at ℓ = 492 (k_{pix} = 21), similar to Fig. 3, but for the further rebinned data. In this case, the number of azimuthal pixels, n_{pix} is reduced to 16, instead of the original 64. Legends and colors are similar to those of Fig. 3. Scale factors for the error estimates by the updated code are adjusted accordingly. 
This figure tells us that our updated code is better than the original for this further rebinned case as well, in terms of a more reasonable error estimate and a smaller bias. Moreover, using our updated code, we can rebin further up to 16 pixels at this k (k_{pix} = 21) without a significant increase in the noise or the underestimate of the parameters or outliers in comparison with the original 64pixel case. In the further rebinned case, the correlations between parameters show a trend that is essentially similar to that of Fig. 4. The only exception is in the fitting results obtained by the modified original code; they show less correlation even in the case of u_{x, n = 3} = 400 m s^{−1} in the further rebinned case.
In the original code, there is no scaling factor related to the bin number and the error estimate by the code is needed to be scaled properly. For these 64pixel and 16pixel cases, the scatters are not significantly changed and it is the case as well for the error estimated with the proper scaling. For these plots, again we note that we need a sufficient number of pixels for a good fitting; in this case 8 or 4 bins were too small, because the functional form of the model has terms that vary as cos θ and sin θ but also as cos2θ and sin2θ as parameters.
4.2. Future work
The modifications described in this work are mainly corrections of problems in the original code. The exception to this is the change of the algorithm from the maximum likelihood method based on the chisquare distribution function to the one based on the normal distribution function, namely the leastsquares method, and the fitting not to the squareroot of power but to the logarithm of the power. There are several potential improvements of the analysis. While implementing such further alterations are beyond the scope of this paper, we will briefly discuss some potential future improvements.
One of the open issues is the remapping and rebinning. Currently the power is remapped from Cartesian to polar coordinates and then smoothed in the azimuthal direction. However, it is possible to do the analysis in the original Cartesian system as it is done in the HMI pipeline fitsc module (Bogart et al. 2011a), which fits in slices at constant temporal frequency. The fitting approach shown here would only need to be slightly modified to be carried out in a region with near k. We expect that the main issue would be extending the model to allow for small variations in k. As most parameters are presumably smooth in k, we speculate the fitting small range in the k rather than single k would help the stability of the fits.
In the current code, we approximate the PDF of the remapped and smoothed power spectrum as a normal distribution function. But how we rebin, including the topic mentioned in Sect. 4.1, is still open. In the method shown here, we do not rebin the model function but use the model function calculated on the same grid as the rebinned data. This is an approximation in the limit where the rebinning does not significantly change the shape of the limit spectrum and we expect that it will cause a bias in the case of extreme rebinning. Figures 3 and 5 show that the bias in the fitting results of the updated code is less than about 5 m s^{−1} in the range of u_{x} ≤ 400 m s^{−1} even in the further rebinned case (Fig. 5).
The choice of the model function (currently Eq. (2)) is also an open issue. As shown in Fig. 1, the current model function does not reproduce the detailed structure of the observed power spectrum. For example, the current model function does not take into account the asymmetry of the ridge shape in terms of frequency. Currently, we use the power spectrum of the tile at the disc center only but that is the simplest case and when we investigate the deep convection, we cannot avoid using the tiles on various locations on the disc. In such a case, in order to construct a model function, we need to take into account further effects, such as the effect of the centertolimb variations (Zhao et al. 2012), the lineofsight effect on shape of the power, or the effect of the Postel projection. Also, the effects of differential rotation on eigenfunctions are to be accounted for in future.
5. Conclusions
We identified several problems in the multiridge fitting code ATLAS (Greer et al. 2014) and we updated the code in response. We confirmed that flowestimate biases and error overestimates exist in the fitting results by the original code. The biases that we found are insufficient to resolve the discrepancy presented in Hanasoge et al. (2016).
The updated code is based on a consistent model and an appropriate likelihood function. Monte Carlo tests show the fitting results and their error estimates by the updated code are reasonable and confirm the improvement of the fitting.
The work shown here is limited to the fitting part of the ringdiagram analysis codes. The next step in the ringdiagram analysis is flow inversion using the modeshift parameters (q_{3, n} and q_{4, n} in Eq. (2)). Unlike the HMI ringdiagram pipeline in which the inversion is done at each tile, Greer et al. (2015) used a 3D inversion using multiple tiles. This unique step should be the focus of examination in future works.
This library code minimizes the sum of square of given function, and ℬ_{i} in Eq. (A.5) is given in the original ATLAS code.
Acknowledgments
We thank Rick Bogart for useful discussions. The German Data Center for SDO (GDCSDO), funded by the German Aerospace Center (DLR), provided the IT infrastructure to process the data. The HMI data used are courtesy of NASA/SDO and the HMI science teams. BWH acknowledges NASA support through grants 80NSSC17K0008, 80NSSC18K1125, and 80NSSC19K0267. We acknowledge partial support from ERC Synergy grant WHOLE SUN 810218.
References
 Anderson, E. R., Duvall, Jr., T. L., & Jefferies, S. M. 1990, ApJ, 364, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Antia, H. M., & Basu, S. 2007, Astron. Nachr., 328, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Appourchaux, T., Andersen, B., & Sekii, T. 2002, in From Solar Min to Max: Half a Solar Cycle with SOHO, ed. A. Wilson, ESA SP, 508, 47 [NASA ADS] [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011a, GONGSoHO 24: A New Era of Seismology of the Sun and SolarLike Stars, 271, 012008 [Google Scholar]
 Bogart, R. S., Baldner, C., Basu, S., Haber, D. A., & RabelloSoares, M. C. 2011b, GONGSoHO 24: A New Era of Seismology of the Sun and SolarLike Stars, 271, 012009 [Google Scholar]
 Brun, A. S., & Browning, M. K. 2017, Liv. Rev. Sol. Phys., 14, 4 [Google Scholar]
 Duvall, Jr., T. L., & Harvey, J. W. 1986, in NATO Advanced Science Institutes (ASI) Series C, ed. D. O. Gough, 169, 105 [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S. M., Harvey, J. W., & Pomerantz, M. A. 1993, Nature, 362, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2005, Liv. Rev. Sol. Phys., 2, 6 [Google Scholar]
 Greer, B. J. 2015, Ph.D. thesis, University of Colorado at Boulder [Google Scholar]
 Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Greer, B. J., Hindman, B. W., & Toomre, J. 2014, Sol. Phys., 289, 2823 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annu. Rev. Fluid Mech., 48, 191 [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Nat. Acad. Sci., 109, 11928 [NASA ADS] [CrossRef] [Google Scholar]
 Harvey, J. 1985, in Future Missions in Solar, Heliospheric & Space Plasma Physics, eds. E. Rolfe, & B. Battrick, ESA SP, 235, 199 [NASA ADS] [Google Scholar]
 Hill, F. 1988, ApJ, 333, 996 [NASA ADS] [CrossRef] [Google Scholar]
 Kosovichev, A. G., & Duvall, Jr., T. L. 1997, in SCORe’96 : Solar Convection and Oscillations and their Relationship, eds. F. P. Pijpers, J. ChristensenDalsgaard, & C. S. Rosenthal (Dordrecht: Springer Science Business Media), Astrophys. Space Sci. Libr., 225, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, ASP Conf. Ser., 411, 251 [NASA ADS] [Google Scholar]
 Marquardt, D. 1963, J. Soc. Ind. Appl. Math., 11, 431 [Google Scholar]
 Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Moré, J. J. 1978, in Numerical Analysis, ed. G. A. Watson (Berlin, Heidelberg: Springer), 105 [Google Scholar]
 Moré, J., & Wright, S. 1993, Optimization Software Guide (Society for Industrial and Applied Mathematics) [Google Scholar]
 Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [Google Scholar]
 Woodard, M. F. 1984, Ph.D. thesis, University of California, San Diego [Google Scholar]
 Zhao, J., Nagashima, K., Bogart, R. S., Kosovichev, A. G., & Duvall, Jr., T. L. 2012, ApJ, 749, L5 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Original ATLAS code
A.1. Fitting method of the original ATLAS code
In the original ATLAS code (Greer et al. 2014, 2015), the processing steps after computing the power spectrum for each tracked tile are as follows:

at each frequency ν, remap the log of the square root of the power spectrum from the Cartesian (k_{x}, k_{y}) to polar (k, θ) coordinates,

compute and carry out Fourier interpolation using a boxcar lowpass filter to smooth in the azimuthal (θ) direction. The number of the grid points in azimuth θ is reduced from n_{pix} = 256 to n_{pix} = 64,

fit a Lorentzian model (see Appendix A.2) to the smoothed squareroot power at each k. The fitting is done based on the maximum likelihood method assuming that the PDF of the power is the chisquare distribution with two degrees of freedom (see more details in Appendix B.1), and

estimate error of the fitting parameters separately from the fitting using the Fisher information matrix based on the final fitting results.
A.2. Model for the power spectrum in the original ATLAS code
In the original ATLAS code (Greer et al. 2014, 2015), the power at each k is modeled with the sum of Lorentzians with seven parameters () for each ridge n = 0, …n_{r} − 1, where n_{r} is the number of ridges, and a background model with five parameters () at each k. The total number of fit parameters at each k is 7n_{r} + 5.
The model power spectrum, at a single k, is:
where
and is the frequency of the nth peak, = A_{n} is the amplitude, = Γ_{n} is the width, is the horizontal velocity, and are anisotropy terms in Eqs. (5) and (6) in Greer et al. (2014). The background is modeled at each k as
where is an amplitude, = ν_{bg} is a rolloff frequency, is the powerlaw index, is the amplitude of the anisotropy term, and is the phase of the anisotropy term in Eq. (7) in Greer et al. (2014).
A.3. Problems in the original ATLAS code
First, the Lorentzian model for the power in the original code was fit to the square root of the observed power. This is inconsistent with what is stated in Greer et al. (2014, 2015) and it is also an inconsistency between the model and the observable. Therefore, we made them consistent in the updated code (see Sect. 2.2).
Second, the cost function to be minimized in the original code was not a good approximation. The cost function to be minimized in the maximum likelihood method based on the chisquaredistribution with two degrees of freedom as PDF is
where O_{i} is the observed spectrum and P_{i} is the model, as given in Sect. A.2. The sum is taken over all grid points (θ_{i}, ν_{i}) in the fitting range.
Instead, the original code minimizes
using the mpfit.c code^{1}. This is conceptually inconsistent with the description in Greer (2015), although the likelihood function itself is not explicitly stated there. While it can be shown that the results of minimizing of the exact cost function (Eq. (B.3)) and the wrong one (Eq. (A.5)) are identical in the limit of linear perturbations, there is no reasonable computational or physical reason to take the extra square in the calculations. Moreover, it is not useful for the error estimates. Therefore, we decided to correct this issue; Sect. 2.2 describes our changes.
Third, in the original ATLAS code, error estimates are obtained by an independent calculation of the Fischer information using the final fitting results and using the chisquare distribution with two degrees of freedom as the likelihood function. It is not consistent to compute the error estimates with a likelihood function that is different than what was used in the fitting itself. In Sect. 2.2, we also describe a consistent approach to the computation of error estimates in our updated code.
Fourth, several parameters are unstable. Therefore, we have changed the parameterizations as we discussed in Sects. 2.3.1 and 2.3.2.
Appendix B: Probability distribution function (PDF) and maximum likelihood method based on the PDF
B.1. PDF of the raw power spectrum
Woodard (1984) demonstrated that the PDF of a single observed power spectrum divided by the expectation value of the spectrum is the chisquare distribution with two degrees of freedom. On this basis, Duvall & Harvey (1986) and Anderson et al. (1990) introduced the probability density at a given grid point in the Fourier space (k_{x}, k_{y}, ν)_{i}
where O_{i} is the observed spectrum, P_{i} is the model for the limit spectrum, and q is the model parameters. The joint probability density for the experimental outcome at horizontal independent wavevectors and frequencies is given by
This is the likelihood function, and the model parameters which maximize this ℒ are the targets in the maximum likelihood method. In practice −lnℒ is minimized to use standard minimization procedures. To make the computation simpler, lnO_{i} is subtracted from −lnℒ. The minimization of −lnℒ − lnO_{i} in terms of q is identical to that of −lnℒ, because O_{i} is not a function of the model parameter q. In conclusion,
is minimized.
B.2. PDF of logarithm of averaged power
The PDF of a wellaveraged power spectrum is also a chisquare distribution but with many more than the two degrees of freedom of the original and given the central limit theorem, it can be approximated by a normal distribution. To carry out the leastsquares fitting based on the normal distribution function and estimate errors of the fitting results, we need the variance of the normal distribution function. We therefore take advantage of the fact that the logarithm of the averaged spectrum obeys a normal distribution function with a constant variance. In Appendix B.3, we show that if we have a spectrum obtained by averaging over N spectra, each normally distributed with the same expectation value M, the logarithm of averaged spectra, y, obeys the normal distribution function :
where , therefore, σ_{N} is a constant over y. We note that the function f(y) is linearized around the mean to derive Eq. (B.4); see Appendix B.3 for details.
In the present case, we fit one k at a time. To do this, the spectra are, at each frequency, linearly interpolated in k_{x} and k_{y} to a circle with radius k_{pix}, where k_{pix} is k in units of bins, and smoothed to n_{pix} azimuths. On this basis, the question is if we can still use a leastsquares fit with a diagonal covariance and a single σ_{N} and, if so, which value of σ_{N} (or equivalently N) in Eq. (B.4) should be used. In particular, it needs to be considered that more than 2πk_{pix} are used for the interpolation and that the averaged values will be correlated.
In the limit of averaging over the entire circle, namely n_{pix} = 1, and in the case with k_{pix} ≫ 1, it might be reasonable to assume that N = 2πk_{pix}α, with α accounting for the fact that the averaging is essentially over an annulus around k_{pix}. Indeed, a MonteCarlo test shows that α ≈ 1.5 gives a very good estimate of the error on the average, which will be shown in Fig. B.1 later in this subsection.
Fig. B.1.
Standard deviation of the logarithm of chisquare distributed white noise rebinned on the annulus with the radius of k_{pix} pixels, σ^{scat}, is plotted against the squareroot of the pixel number on the annulus after rebinning, , at three k_{pix} shown with different symbols in different colors. In the case of n_{pix} ≪ 2πk_{pix}, σ^{scat} = σ_{α}(α = 1.5), which is shown with the dotted lines, is a valid approximation, while in the limit of n_{pix} ≳ 2πk_{pix}, σ^{scat} deviates from σ_{α} and approaches a constant of 2/3. 
For n_{pix} ≫ 1 the variances, as well as the offdiagonal elements of the covariance matrix of the interpolated datapoints will, in general, depend on azimuth. Assuming that the fitted function may be linearized in the fitted parameters, a linear fit implies that the fitted parameters are given by a linear combination of the observed values. Assuming that the functions are smooth (as in the present case where they are loworder harmonic functions), the coefficients in the linear combination are also smooth. From this it follows that if the variations in the properties of the covariance matrix average on the scale of the variations in the fitted functions, then the same scaling factor may be used and that
which a MonteCarlo test again confirms.
To validate the arguments above, a simple MonteCarlo test is carried out to measure the scatter (standard deviation, σ^{scat}) of the logarithm of averaged white noise. Whitenoise fields with the chisquare distribution with two degree of freedom in a threedimensional Cartesian Fourier space (k_{x}, k_{y}, ν) = (384, 384, 1152) [pix] are remapped at several annuli with specific radii k_{pix} in the same way as the data in the updated analysis code. At first, there are 256 pixels on each annulus after remapping, we then rebin them into n_{pix} = 4, 8, 16, 32, 64, and 128 pixels. The standard deviation of the logarithm of the rebinned white noise, σ^{scat}, are computed and plotted against the square root of the number of datapoints after rebin on the annuli, n_{pix} = 14, 21, and 42 in Fig. B.1. As we mentioned above, in the case of n_{pix} ≪ 2πk_{pix}, σ^{scat} ≃ σ_{α}(α = 1.5), where . In the case of n_{pix} ≳ 2πk_{pix} the deviation of σ^{scat} from σ_{α} is not negligible, and σ^{scat} approaches a constant 2/3; this constant can be derived from the bilinear interpolations of chisquare distributed random variables in two dimensions. Also we note that we have not included the apodization effect in the data to create the noise field in this test calculation. Without apodization σ^{scat} might be overestimated; in any case, this does not affect the parameter estimates but only the error estimates.
In our updated code, σ_{N} is set to 1 (see the cost function Eq. (1)), as we mentioned in Sect. 2.2, and after the fitting is done the error provided by the code is scaled by σ_{α}. Although calling the minimization code with an incorrect error estimate can lead to worse convergence, the difference in this case should be negligible. Implementing σ_{α} in the code is a task for the future.
B.3. Derivation of PDF of the logarithm of averaged power
In this section we briefly summarize how to derive the PDF of the logarithm of the power.
Suppose x is the average of N spectra whose averages are M. If N is large enough, then the PDF of this spectrum is the normal distribution
and . Here we define
and when x is close to M,
Therefore, the inverse function of g is given by
and
Hence, the PDF of y = lnx is given by
where . Therefore, the logarithm of the averaged power has a constant error, which is independent of the original σ and depends only on the number of the averaged spectra, N. For comparison, approximated PDF and the distribution functions of the MonteCarlo samples are given in Appendix B.4.
B.4. Comparison of the distribution functions with Gaussian and chisquare distribution functions
In this subsection we show that Gaussian function is a good approximation for the distribution functions of the logarithm of the synthesized power. We also discuss the degrees of freedom of the chisquare distribution functions of the unwrapped and smoothed power here.
Figure B.2 shows the distribution function of the power and the logarithm of the power at k_{pix} = 21 (ℓ = 492) from the Monte Carlo simulations. The Gaussian function centered at zero whose width is the standard deviation of the samples are overplotted in the figure. These plots show that the Gaussian function is a good approximation for the distribution function of the logarithm of the unwrapped and smoothed power near the mean value, though the tails of the distribution function are longer than those of the Gaussian.
Fig. B.2.
Normalized frequency distribution functions of the normalized synthetic power at k_{pix} = 21 (ℓ = 492, black solid lines): top panels are for the subsampled power and the second panels are the unwrapped power without further smoothing, while the lower panels are for the rebinned powers. Left panels: (P − ⟨P⟩)/⟨P⟩, and right panels: lnP − ln⟨P⟩, where ⟨⟩ indicates the expectation value (model). The input model of the power is isotropized (q_{i, n} = 0 for 3 ≤ i ≤ 6 for all n and q_{i, b} = 0 for i = 2, 3) but otherwise it is the same as the model with the parameters given in Table C.2. For this plot, we use 50 realizations for the top two rows and 500 realizations for the rebinned data (lower three rows). The black dashed lines on the lower four panels are the Gaussian function centered at zero whose width is the standard deviation of the samples (shown on the panels). This validates our approximation of the distribution function of the logarithm of the power by the Gaussian function. The thick red dashed curves on the left panels are the chisquare distribution function of various degrees of freedom. The degree of freedom (d.o.f.) for each distribution is determined from the standard deviation of each set of normalized power. See the text for details. 
In the left panels of Fig. B.2, the distribution function of the power is also compared with the chisquare distribution functions with various degrees of freedom: the chisquare distribution functions shown here are
where X = (P − ⟨P⟩)/⟨P⟩, ⟨⟩ indicates the expectation value, and d is the number of degrees of freedom (d.o.f.). Since in this case ⟨X⟩ = 0 and the variance of X is 2/d, we determine d.o.f. from the variance and show the corresponding chisquare distribution function in the figure. As the original synthesized power is created from the model and the twod.o.f. chisquare distributed random numbers, subsampled power is purely two d.o.f. The unwrapped power has a fourd.o.f. chisquare distribution; this indicates that unwrapped power is effectively the subsampled power rebinned over two pixels. Because in this case the radius of annulus is k_{pix} = 21 and hence roughly 2πk_{pix} independent pixels are on the annulus, the unwrapped power with 256 pixels on the annulus is roughly twice oversampled. Therefore, the d.o.f. of the power unwrapped and rebinned over two pixels in azimuth (n_{pix} = 128, third panel) is five. This is almost unchanged from the original (second panel, d.o.f. = 4). This is also the case when we further rebin the data over four pixels (fourth panel); d.o.f is only eight and this is smaller than 16, which we might have expected for the fourpixel rebinned fourd.o.f. distributed independent variables (i.e., unwrapped power). But we should also note here that 12 is the d.o.f. which we expect on the basis of the effective rebinning numbers we discussed in Appendix B.2 (Eq. (B.5)), N ∼ 3, where α = 1.5. It is also larger than the dof of the distribution function (8); this is partly because the approximation (Eq. (B.5)) is not sufficiently good for n_{pix} = 64 case as was shown in Fig. B.1. When we see the further rebinned case with n_{pix} = 16, namely the unwrapped power averaged over 16 pixels (the lowermost panel), the d.o.f. is 25 and close to what we expect from Eq. (B.5), N ∼ 6: in this case the d.o.f. is 6 × 4 = 24. For this case the approximation by Eq. (B.5) with α = 1.5 is good according to Fig. B.1. This closeness of the distribution function of the remapped and smoothed power to the chisquare distribution with some degrees of freedom suggest that one might expect improvement to the fitting method by using the maximum likelihood method based on the chisquare distribution with specific degrees of freedom.
Appendix C: Input parameters for smaller and larger wavenumbers
The parameters at several ℓ obtained from the fitting of the observed spectra averaged over one Carrington rotation 2211 for the disc center tile from SDO/HMI pipeline are shown in Tables C.1 (ℓ = 328), C.2 (ℓ = 492), and C.3 (ℓ = 984). These average power spectra were obtained from DRMS specification hmi.rdvavgpspec_fd15[2211][0][0] and the fitting was done using the updated code. These parameters are used as input parameters to generate MonteCarlo synthetic data in Sect. 3.
Input parameter set at ℓ = 328 (k_{pix} = 14) obtained from the fitting the average power spectrum over one Carrington rotation to the model by the updated code.
Input parameters at ℓ = 492 (k_{pix} = 21) similar to Table C.1.
Input parameter set at ℓ = 984 (k_{pix} = 42), similar to Table C.1.
Appendix D: Underestimated amplitude – Logarithm of average and average of logarithm
Figure 2 shows that the amplitudes (q_{1, n} for all n and q_{0, b} in Eq. (2)) measured by the updated code are underestimated. This is the result of taking the logarithm of the rebinned power and the effect can be reproduced in the following simple test.
To see why this is the case, we define x_{i} (i = 0, 1, …n − 1) as a series of random variables whose distribution function is the chisquare distribution with two degrees of freedom. The expectation values and the variance of x are ⟨x⟩ = 2 and , respectively, and y_{j}(j = 0, 1, …n/a − 1) is smoothed x_{i} over each apixel range:
In this case ⟨y⟩=⟨x⟩ = 2, and ln⟨y⟩=ln(2), but ⟨ln(y)⟩ ≤ ln⟨y⟩. The more the array is rebinned (the larger a becomes), the less the scatter of y becomes: ⟨ln(y)⟩ approaches ln⟨y⟩.
This trend is alleviated if the spectra are more strongly averaged. Figure D.1 shows a simple test calculation to illustrate how the average of logarithm is reduced from the logarithm of the average; the ratio of ⟨ln(y)⟩ to ln⟨y⟩ is plotted against the pixel number ratio of the rebinned data to the original. In this simple test calculation, we make 500 realizations of sets of x_{i} (i = 0, 1, …n − 1, where n = 256), calculate y with a = 1, 2, 4, 8, …256, and compute how much smaller the expectation value of the logarithm of y than the logarithm of the expectation value is. The linear regression of the data points are also shown in the plot. For comparison the amplitude fitting results shown in Fig. D.2 are replotted here with squares. This explains how the logarithm of amplitude is reduced. In the case of no rebinning (X = n_{rebin}/n_{original} = 1) or only 2pixel rebin (X = 1/2) the fitting results are deviated significantly from this simple test calculation, but this comes as no surprise because in these cases, the assumption of the wellaveraged power as data is not appropriate.
Fig. D.1.
exp⟨ln(y)⟩/⟨y⟩, where y is the aelement average of chisquare distributed random variables, as a function of the pixel number ratio of the rebinned data to the original (asterisks). The average ratio of the amplitude fitting results to the input shown in Fig. D.2 is depicted by squares. The dashed line and the equation on the panel are the linear regression of the data points. 
Figure D.2 shows the dependence of the amplitude measurements on the amount of rebinning. The input model is identical to the one used in the calculation in Sect. 3.1. The default rebinning is from 256 pixels to 64 pixels and in this case, the output amplitude is about 88% of the input amplitude but with further rebinning to 16 pixels, it increases up to 96%.
Fig. D.2.
Dependence of the peak amplitude parameters at ℓ = 492 (k_{pix} = 21) measured by the updated code on the rebinning strength. Left panel: input (black dotted lines with squares) and measurements with three different rebinning with error bars (scatter of 500 realizations): no rebinning (256 pixels on the azimuthal grid, blue), original rebinning (64 pixels, green), and extra rebinning (16 pixels, red). Right panel: ratio of the amplitude to the input. The numbers are the grid number and the average ratio of six peaks. The input parameters are identical to those in Fig. 2. 
All Tables
Background parameters in the input model for the Monte Carlo simulation and the fitting results at ℓ = 492 (k_{pix} = 21).
Input parameter set at ℓ = 328 (k_{pix} = 14) obtained from the fitting the average power spectrum over one Carrington rotation to the model by the updated code.
All Figures
Fig. 1.
Slices through the input limit spectrum (black) and observational spectrum averaged over one Carrington rotation (gray), which was used to obtain the input parameters, at θ = 0 and a few different values of k. Vertical dashed lines indicate the fitting ranges and the models are plotted only within these ranges. Lower limits for all k are fixed at 0.4 mHz, while the upper limits depend on the initial guess for each k, and they are the highest peak frequencies plus the widths of the peak. For the cases of ℓ = 328 and 492 (k_{pix} = 14 and 21, respectively), the peaks with the radial order, n, from 0 to 7 are used in the fitting, while for the case of ℓ = 984 (k_{pix} = 42), 0 ≤ n ≤ 3. 

In the text 
Fig. 2.
Fitting parameters and their errors for the case ℓ = 492 (k_{pix} = 21). Each panel has three parts. In the upper part of each panel, squares indicate the input parameters, . The asterisks indicate the fitting results obtained by the updated code, namely, the mean of MonteCarlo samples, . In the middle part of each panel, the deviation of the mean fitting results from the input, , are illustrated by the circles and the expected scatter of the mean computed from the square root of variance σ^{scat} of the 500 Monte Carlo samples divided by the square root of the sample number, N_{sample} (here it is 500) are shown as the error bars. The dashed horizontal lines are at δq = 0. In the lower part of each panel the scatter of the 500 samples, σ^{scat} and the scaled error estimated by the updated code, σ^{code}, are depicted by the thick gray lines and the thin lines with short horizontal bars on the edge. Errors estimated by the updated code are scaled by σ(α = 1.5, k_{pix}, n_{pix}) = 0.569 as described in Appendix B.2. The error bars on the middle panel of δq_{1, n} are tiny at this scale; the underestimation is relatively large. However, the errors of δq_{1, n} are and, therefore, available from σ^{scat} in the lower panel and N_{sample} = 500; they are ∼0.02 at most. 

In the text 
Fig. 3.
Fitting results of u_{x, n = 3} at ℓ = 492 (k_{pix} = 21) in the form of normalized histograms by the original code which fit the model to the square root of power (black), the modified original code which fit the power (blue), and the updated code (red). Horizontal lines indicate the average of 500 Monte–Carlo realizations given by the three codes (in the same colors as the histograms, dashdotted black, dashed blue, and solid red). The short vertical lines on the horizontal lines indicates the error estimated by the codes (thin lines with short horizontal bars) and the standard deviations of the 500 fitting results (thick lines) in the same color as the histograms centered at the means (horizontal lines) by the three codes. We note that the fitting results with overly large deviation are omitted. See text for details. 

In the text 
Fig. 4.
Correlation coefficients between the parameters (500 MonteCarlo realizations) for the isotropic model (u_{x, n} = 0 for all n) and for the isotropic model except with u_{x} = 400 m s^{−1} for the n = 3 peak (on the next page) at ℓ = 492 (k_{pix} = 21). Each box indicates the parameters for nth peak (n = 0, 1, …5), q_{i, n} (i = 0, …6, from left to right and from bottom to top in each box) and the background parameters (BG). As we defined in Sect. 2q_{0, n} = ν_{n} is the frequency of the nth peak, q_{1, n} = A_{n} is the amplitude, q_{2, n} = Γ_{n} is the width, q_{4,n}) = u = (u_{x,n}, u_{y,n}) is the horizontal velocity, q_{5, n} = f_{c, n} and q_{6, n} = f_{s, n} are parameters to handle anisotropy. Background parameters are q_{0, BG} = B_{0} is the amplitude, q_{1, BG} = b is the powerlaw index, q_{2, BG} = f_{c, bg} and q_{3, BG} = f_{s, bg} are the parameters to handle the anisotropy of the background. Color scale is shown by the color bar on the right side. We note that the last columns and rows (i = 4) of BG on the top panels are empty because the number of the background parameters in the model in the updated code is four, instead of five. 

In the text 
Fig. 4.
continued. 

In the text 
Fig. 5.
Fitting results of u_{x, n = 3} at ℓ = 492 (k_{pix} = 21), similar to Fig. 3, but for the further rebinned data. In this case, the number of azimuthal pixels, n_{pix} is reduced to 16, instead of the original 64. Legends and colors are similar to those of Fig. 3. Scale factors for the error estimates by the updated code are adjusted accordingly. 

In the text 
Fig. B.1.
Standard deviation of the logarithm of chisquare distributed white noise rebinned on the annulus with the radius of k_{pix} pixels, σ^{scat}, is plotted against the squareroot of the pixel number on the annulus after rebinning, , at three k_{pix} shown with different symbols in different colors. In the case of n_{pix} ≪ 2πk_{pix}, σ^{scat} = σ_{α}(α = 1.5), which is shown with the dotted lines, is a valid approximation, while in the limit of n_{pix} ≳ 2πk_{pix}, σ^{scat} deviates from σ_{α} and approaches a constant of 2/3. 

In the text 
Fig. B.2.
Normalized frequency distribution functions of the normalized synthetic power at k_{pix} = 21 (ℓ = 492, black solid lines): top panels are for the subsampled power and the second panels are the unwrapped power without further smoothing, while the lower panels are for the rebinned powers. Left panels: (P − ⟨P⟩)/⟨P⟩, and right panels: lnP − ln⟨P⟩, where ⟨⟩ indicates the expectation value (model). The input model of the power is isotropized (q_{i, n} = 0 for 3 ≤ i ≤ 6 for all n and q_{i, b} = 0 for i = 2, 3) but otherwise it is the same as the model with the parameters given in Table C.2. For this plot, we use 50 realizations for the top two rows and 500 realizations for the rebinned data (lower three rows). The black dashed lines on the lower four panels are the Gaussian function centered at zero whose width is the standard deviation of the samples (shown on the panels). This validates our approximation of the distribution function of the logarithm of the power by the Gaussian function. The thick red dashed curves on the left panels are the chisquare distribution function of various degrees of freedom. The degree of freedom (d.o.f.) for each distribution is determined from the standard deviation of each set of normalized power. See the text for details. 

In the text 
Fig. D.1.
exp⟨ln(y)⟩/⟨y⟩, where y is the aelement average of chisquare distributed random variables, as a function of the pixel number ratio of the rebinned data to the original (asterisks). The average ratio of the amplitude fitting results to the input shown in Fig. D.2 is depicted by squares. The dashed line and the equation on the panel are the linear regression of the data points. 

In the text 
Fig. D.2.
Dependence of the peak amplitude parameters at ℓ = 492 (k_{pix} = 21) measured by the updated code on the rebinning strength. Left panel: input (black dotted lines with squares) and measurements with three different rebinning with error bars (scatter of 500 realizations): no rebinning (256 pixels on the azimuthal grid, blue), original rebinning (64 pixels, green), and extra rebinning (16 pixels, red). Right panel: ratio of the amplitude to the input. The numbers are the grid number and the average ratio of six peaks. The input parameters are identical to those in Fig. 2. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.