Damping of Propagating Kink Waves in the Solar Corona

Before examining the velocity signals from CoMP, it is beneficial to understand the geometry of the loop system that will be considered. To provide some insight, the potential field source surface (PFSS—Schrijver & DeRosa 2003) extrapolation package available in SolarSoft is used to provide an indication of the local magnetic field structure in the corona (Figure 1: left-hand panel). To determine the validity of the obtained field extrapolations, several attempts to generate extrapolations in the neighborhood of the footpoints are undertaken and we find that the given PFSS loops are indeed unique. Further extrapolations were undertaken to examine the solution for constant latitudinal points and to ascertain that the loops we obtained from the initial extrapolations are the best representation for the observed CoMP loops. The extrapolated field lines obtained after these initial checks are shown (Figure 1: left-hand panel) and visually represent the coronal structures well. There is also a close agreement with the direction of wave propagation determined from CoMP, which is believed to follow the magnetic field lines (Figure 1 center and right-hand panels, see Section 3.2 for further details). The projection of the field lines onto the magnetogram is shown in Figure 2.

Zoom In Zoom Out Reset image size

The CoMP Doppler velocity image sequence shows coherent fluctuations propagating through the corona, which are interpreted as propagating kink waves. We begin by determining the direction of the propagation of the waves. The coherence between the velocity time-series of each pixel and its neighboring pixels is calculated using an FFT-based method (McIntosh et al. 2008; Tomczyk & McIntosh 2009) and correlation maps are derived. Selecting pixels in the neighborhood where the coherence value is greater than 0.5 defines a coherence island. This coherence island has a distinct direction, following the apparent trajectory of the propagating waves. The direction of wave propagation is then taken to be aligned with the island, which is determined by fitting a line that minimizes the sum of perpendicular distances from the points to the line. This is performed for each pixel in the field-of-view, enabling us to create a wave angle map, which is displayed in the right-hand panel of Figure 1. The shown angle gives the direction of propagation measured counterclockwise from a due east direction. Given that the kink mode propagates along the magnetic field, this angle should also represent the magnetic field orientation in the plane-of-sky (POS). This method does indeed show excellent agreement with polarimetric measurements of the POS direction of the magnetic field (Tomczyk & McIntosh 2009).

The wave angle map is used to determine the path of the wave propagation through the corona, enabling the kink wave packets to be followed and allowing us to determine how they evolve as they propagate. We select five different wave paths with increasing lengths (center panel of Figure 1), where the selected paths are assumed to follow the quiescent coronal loops and, to satisfy the restrictions of Equation (4), are assumed to represent half the total loop length (this assumption is discussed further in Section 4). The velocity signal along the wave paths is extracted to create time–distance maps, where cubic interpolation is used to map the velocities from the selected wave paths onto (x, t) space.¹ For each wave path shown in Figure 1, we also extract the neighboring five wave paths on either side of the original wave path. Each additional path is calculated using the normal vector to the original, and are separated by one pixel in the perpendicular direction.

These velocity time–distance maps are composed of both the inward and outward propagating kink waves. Taking a Fourier transform of the velocity time–distance maps enables us to produce the k–ω spectra for velocity power, as shown in Figure 3. The wave power is separated for the inward and outward components of the wave propagation. It is evident from all k–ω spectra that the outward wave power dominates over the inward wave power. Filtered time–distance diagrams are created by taking the inverse Fourier transform of the inward and outward halves of the k–ω spectra separately (Tomczyk et al. 2007; Tomczyk & McIntosh 2009; Morton et al. 2015). The filtered time-series are used to obtain the wave propagation speed along the wave path for both outward and inward propagating waves. The time-series at the center of the wave path are cross-correlated to the neighboring time-series along the path. The lag of the cross-correlation is determined by fitting a parabola to the peak of the correlation function. The propagation speed is then calculated by fitting the slope of the observed lags as a function of the position along the wave path.

Zoom In Zoom Out Reset image size

Finally, the wave power as a function of frequency for the inward and outward components is calculated by summing the spectra in the k-direction. For each loop, the inward and outward spectra are averaged over the neighboring wave paths to suppress the variability (see Section 3.4 and Appendix A for further discussion). From this one-dimensional averaged wave power, the ratio of the outward and inward power, $\langle P(f){\rangle }_{\mathrm{ratio}}$, is determined. For each of the coronal loops that we have studied, the ratio of power spectra displays an increase in magnitude as frequency increases (Figure 3). This signature demonstrates the frequency dependence of the change in outward and inward power, indicating that a frequency-dependent process is in action to attenuate the waves (e.g., resonant absorption).

As discussed in the introduction, to estimate the quality factor from the obtained power ratio, the model power ratio given by Equation (4) should be fitted to the data in a robust manner. Therefore, it is necessary to discuss the statistics of the power ratio.

In Verth et al. (2010), the ratio of the outward and inward spectra were fitted with the model given by Equation (4) using a least-squares minimization. However, as we will show, the assumption that the power ratio values at each frequency ordinate are normally distributed (implicit in least-squares) is incorrect and leads to a poor estimate of model parameters and their uncertainties. Here, we present a new method for the MLE of model parameters from the ratio of two power spectra obtained via a DFT.

The power spectra, I (f_i) at each frequency ordinate, f_i; i = 0, 1, 2, 3, ..., n, from the DFT are distributed about the true power value, P(f_i) as

$$\begin{eqnarray}&&I({f}_{i})=P({f}_{i})\displaystyle \frac{{\chi }_{2}^{2}}{2}.\end{eqnarray} \tag{ 5 }$$

Here ${\chi }_{2}^{2}$ represents a random variable from the chi(χ)-squared distribution with two degrees of freedom, distributed as

$$\begin{eqnarray}&&{\chi }_{2}^{2}=\displaystyle \frac{1}{2}\exp \left(-\displaystyle \frac{x}{2}\right)\end{eqnarray} \tag{ 6 }$$

(see e.g., Vaughan 2005). Suppose now that we are interested in taking the ratio of the values x and y, drawn from two independent ${\chi }_{2}^{2}$ distributions X and Y. The associated probability distribution function (PDF) is Z = X/Y and the distribution of Z is then given by

$$\begin{eqnarray}&&{\psi }_{z}={\int }_{0}^{\infty }y{\psi }_{{xy}}({zy},y){dy}.\end{eqnarray} \tag{ 7 }$$

Given that x and y are independent, ψ_xy is given by

$$\begin{eqnarray}&&{\psi }_{{xy}}=\displaystyle \frac{1}{4}\exp \left(-\displaystyle \frac{x+y}{2}\right).\end{eqnarray} \tag{ 8 }$$

Hence,

$$\begin{eqnarray}&&{\psi }_{z}=\displaystyle \frac{1}{{(1+z)}^{2}},\end{eqnarray} \tag{ 9 }$$

and the distribution of the ratio of any two given power spectra, z (i.e., ratio of ${\chi }_{2}^{2}$ distributions) is given by the log-logistic distribution (Equation (9)). For a non-normalized random variable, r, one can obtain the probability distribution by change of variable, introducing

$$\begin{eqnarray}&&z=\displaystyle \frac{r}{s}\end{eqnarray} \tag{ 10 }$$

where s is the appropriate normalizing factor. The resulting PDF is given by

$$\begin{eqnarray}&&g(r)=\psi \left(\displaystyle \frac{r}{s}\right)\displaystyle \frac{{dz}}{{dr}}.\end{eqnarray} \tag{ 11 }$$

Hence,

$$\begin{eqnarray}&&g(r)=\displaystyle \frac{1}{s}\displaystyle \frac{1}{{\left(\tfrac{r}{s}+1\right)}^{2}}.\end{eqnarray} \tag{ 12 }$$

For the power spectra ordinates, we know that $2I/P$ is ${\chi }_{2}^{2}$. Hence, if $x=2{I}_{1}/{P}_{1}$ and $y=2{I}_{2}/{P}_{2}$ then, $z={I}_{1}{P}_{2}/{I}_{2}{P}_{1},$ and $r={I}_{1}/{I}_{2}$, $s={P}_{1}/{P}_{2}$. Thus, the PDF of the ratio of the power spectra ordinates is calculated to be

$$\begin{eqnarray}&&g({R}_{i})=\displaystyle \frac{1}{{S}_{i}}\displaystyle \frac{1}{{\left(\tfrac{{R}_{i}}{{S}_{i}}+1\right)}^{2}},\end{eqnarray} \tag{ 13 }$$

where ${R}_{i}={I}_{1i}/{I}_{2i}$ , is the power spectra ratio and ${S}_{i}={P}_{1i}/{P}_{2i}$ is the true ratio of the spectral power.

In this study, several power spectra are summed, which changes the distribution by altering the number of degrees of freedom. The ratio of two ${\chi }_{\nu }^{2}$ distributed variables can be shown to be distributed following the F-distribution, given by

$$\begin{eqnarray}&&F(z;\nu ,\varphi )=\displaystyle \frac{1}{\beta \left(\tfrac{\nu }{2},\tfrac{\varphi }{2}\right)}{\left(\displaystyle \frac{\nu }{\varphi }\right)}^{\tfrac{\nu }{2}}{z}^{\tfrac{\nu }{2}-1}{\left(1+z\displaystyle \frac{\nu }{\varphi }\right)}^{-\tfrac{\nu +\varphi }{2}}.\end{eqnarray} \tag{ 14 }$$

where ν and φ are the degrees of freedom (number of parameters), and β is the beta function. The log-logistic distribution is recovered for ν = φ = 2. For ν = φ, the F-distribution simplifies to

$$\begin{eqnarray}&&F(z;\nu ,\nu )=\displaystyle \frac{1}{\beta \left(\tfrac{\nu }{2},\tfrac{\nu }{2}\right)}{z}^{\tfrac{\nu }{2}-1}{\left(1+z\right)}^{-\nu }.\end{eqnarray} \tag{ 15 }$$

The F-distribution is an asymmetric distribution with a minimum value 0 and no maximum value. In Figure 4 we show the nature of the distribution for various values of the degrees of freedom ν and φ. There is a different F-distribution for each combination of these two degrees of freedom. The distribution is heavily right-skewed for smaller values of ν and φ, which means there is a long tail and an increased chance of more extreme large values. As the degrees of freedom increase, the F-distribution is more localized.

Zoom In Zoom Out Reset image size

As previously, substituting in the normalized variables gives

$$\begin{eqnarray}&&{g}_{\nu ,\nu }\left(\displaystyle \frac{{R}_{i}}{{S}_{i}}\right)=\displaystyle \frac{1}{{S}_{i}}\displaystyle \frac{1}{\beta \left(\tfrac{\nu }{2},\tfrac{\nu }{2}\right)}{\left(\displaystyle \frac{{R}_{i}}{{S}_{i}}\right)}^{\tfrac{\nu }{2}-1}{\left(1+\displaystyle \frac{{R}_{i}}{{S}_{i}}\right)}^{-\nu }.\end{eqnarray} \tag{ 16 }$$

Assuming a model S(θ) for the true power ratio, with unknown parameters θ, the joint probability density of observing N periodogram ratio points R_i is given by the likelihood function, ${ \mathcal L }$, where

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{i=1}^{n}p({R}_{i}| {S}_{i})=\displaystyle \prod _{i=1}^{n}\displaystyle \frac{1}{{S}_{i}}F\left(\displaystyle \frac{{R}_{i}}{{S}_{i}};\nu ,\nu \right).\end{eqnarray} \tag{ 17 }$$

Maximizing the likelihood is equivalent to minimizing the negative of the log of the likelihood function, namely

$$\begin{eqnarray}\begin{array}{rcl}-2\mathrm{ln}{ \mathcal L } & = & 2\displaystyle \sum _{i=1}^{n}\left[\mathrm{ln}{S}_{i}+\mathrm{ln}\beta \left(\displaystyle \frac{\nu }{2},\displaystyle \frac{\nu }{2}\right)\right.\\ & & +\ \left.\left(1-\displaystyle \frac{\nu }{2}\right)\mathrm{ln}\left(\displaystyle \frac{{R}_{i}}{{S}_{i}}\right)+\nu \mathrm{ln}\left(1+\displaystyle \frac{{R}_{i}}{{S}_{i}}\right)\right].\end{array}\end{eqnarray} \tag{ 18 }$$

We have verified that this likelihood function provides consistent estimators for the model parameters θ (see Appendix A).

The observed power ratios shown in Figure 3 are then used to estimate the model parameters for the power ratio, i.e., the power ratio scaling factor, P_out/P_in and the factor in the exponential, $2L/{v}_{\mathrm{ph}}\xi$ given in Equation (4). We use the Powell method for minimization, making use of the IDL POWELL function (e.g., Barret & Vaughan 2012).

The associated confidence intervals on the model parameters can be estimated by utilizing the Fisher Matrix $({ \mathcal F })$. The components of ${{ \mathcal F }}_{{ij}}$ are defined as the expected value of the Hessian $({ \mathcal H })$

$$\begin{eqnarray}&&{{ \mathcal F }}_{{ij}}=\left\langle -\displaystyle \frac{{\partial }^{2}\mathrm{ln}{ \mathcal L }}{\partial {\theta }_{i}\partial {\theta }_{j}}\right\rangle ,\end{eqnarray} \tag{ 19 }$$

where θ represents the model parameters (Pawitan 2001; Bevington & Robinson 2003). The Fisher matrix is a N × N matrix for N model parameters. The inverse of the Fisher Matrix gives the covariance matrix, the diagonal elements of which give the standard error squared on each model parameter, σ². The off-diagonal matrix elements provide the covariances between parameters. The Fisher Matrix only gives reliable uncertainties when the likelihood surface can be approximated by a multi-dimensional Gaussian. Here we will give the values obtained from the covariance matrix as the estimated parameter uncertainties, and have checked that they are in close agreement with more involved methods of calculating confidence levels; e.g., Wilks confidence intervals (Bevington & Robinson 2003). At best, the given uncertainties and confidence intervals should be taken as a lower limit.

We use the standard errors to calculate the point-wise Wald 95% confidence intervals (Bevington & Robinson 2003) for the model. The likelihood surface and covariance matrix suggest covariance between the model parameters and this is included in the confidence interval calculation. For the measured power ratios given in Figure 3, the likelihood surfaces are close to a bivariate Gaussian, thus the corresponding confidence bands calculated are reliable. We note that in the case of the ratio of two single (i.e., non-averaged) power spectra (ν = 2), the likelihood function is irregular and the Fisher Matrix will likely provide a poor coverage of the confidence intervals (see examples in Appendix A).

Potential field extrapolations are undertaken with PFSS to determine the geometry for the quiescent coronal loop system shown in Figure 1. In particular, we are keen to examine whether the wave paths determined from following the Alfvénic fluctuations are situated in the POS, which has been the implicit assumption in previous analyses (Tomczyk & McIntosh 2009; Verth et al. 2010). This assumption has an impact on the measured propagation speeds and the lengths of the loops, both of which are important quantities for determining the equilibrium parameter ξ from the data (see Equation (4)). We note that the plotted magnetic field lines in Figures 1 and 2 are not supposed to represent the specific coronal loops along which we believe the waves are propagating. However, we expect that the extrapolated field will represent the general behavior of the magnetic field in the region, and, as such, will describe the oscillating loops. Moreover, as mentioned earlier, the spatial resolution of CoMP essentially precludes the identification of the individual coronal structures. The extrapolated field demonstrates that the loops are approximately situated in the POS, with a maximum angle between the loops and POS found to be 20° (see Figure 2).

Furthermore, it can be inferred from the extrapolated field lines plotted in the left-hand panel of Figure 1 that the geometry of the coronal loops is not symmetric about the apexes. Given that the model used for fitting the wave damping is derived under the assumption that both the outward and inward waves have propagated along half of the loop (Equation (4)), this will likely affect our estimates for ξ (discussed further in Section 4.3). In Appendix B, we give a more general model for the exponential damping that can be fitted to the data when measuring over a segment of the loop. Although knowledge of total loop length and the segment length are required, for this dataset there   are no stereoscopic data available that would help us to achieve this. Moreover, we are hesitant in trying to determine any one-to-one correspondence between the extrapolated field lines and the wave angle guided tracks. Given that the main purpose of this work is to present a more appropriate method for fitting the observed power ratio and demonstrate that the least-squares method gives incorrect model parameters, this limitation does not invalidate our aim. Hence, the general formula is provided for future work and we ask that these limitations are kept in mind as we proceed with the analysis.

Using two-dimensional DFTs, we determine the inward and outward components of the wave power corresponding to the Alfvénic waves propagating along five wave paths of increasing length (the wave paths are shown in Figure 1 center panel). The k–ω diagrams for three of the wave paths are displayed in Figure 3 (left-hand column) and they provide an indication of the relative strength of the outward and inward propagating Alfvénic waves in the segment of loop under consideration. The k–ω diagrams have the distinct ridges that have been reported in previous observations, corresponding to the near dispersion-less kink mode, where the negative frequencies correspond to outward waves and positive are inward waves. Given that the spatial frequency resolution is lower for the shorter loops (Figure 3 top left) compared to the longer loops (Figure 3 bottom left), the k–ω diagrams are less well resolved for the shorter loops. Nevertheless, it can be noticed that as the length of the loop increases, the relative power in the outward propagating Alfvénic waves to the inward propagating waves increases. Assuming that the Alfvénic waves entering the corona at both footpoints of the loops have the same power spectra, then this potentially has a trivial explanation: for longer loops, the inward propagating waves will have traveled further distances and they should be expected to have been damped to a greater degree, as suggested by Equation (4). The plots in the right-hand columns of Figure 3 are obtained by collapsing the spectra in the wave number direction and taking the ratio of the outward to inward spectra. The power ratio shows an apparent upward trend as a function of frequency, which indicates wave damping, with the relative magnitudes of the power ratio supporting the visual impression from the k–ω diagrams, which indicates greater wave damping for the longer loops. Following Verth et al. (2010), we only show the ratio of power spectra up to 4 mHz. This is largely because the signal drops below the noise level for the inward propagating waves beyond this frequency and leads to a turnover in the power spectra.

Using the derived likelihood function (Equation (18)), we are able to fit the power ratio model (Equation (5)) to the data points shown in Figure 3. The maximum likelihood model parameters are used to define the model power ratio curve (red-solid line right-hand column of Figure 3). The values in the covariance matrix enable us to generate the point-wise Wald confidence bands at 95% via bootstrapping. The confidence bands demonstrate that, in each case, there is a clear trend in the power ratio as a function of frequency. This supports the idea that frequency-dependent wave damping is in action along each wave path (Verth et al. 2010).

Given that previous work has employed the least-squares method for fitting the power ratio model, we also demonstrate the differences between the parameter estimates from least-squares and MLE methods. In Figure 3 (right-hand column) the model curves obtained from the least-squares (black solid line) are overplotted and demonstrate that they underestimate the amount of damping present, i.e., corresponding to flatter curves, when compared to the MLE method.

To estimate the equilibrium parameters (quality factor), ξ, for each selected wave path, we also require the length of the wave path used, assumed to be the half-loop length (L), and the propagation speed of the Alfvénic waves (υ_ph). The values are summarized in Table 1. The measured propagation speeds of ≈680 km s⁻¹ are consistent with the values obtained in previous studies (Tomczyk & McIntosh 2009; Morton et al. 2015). This value is averaged over the outward and inward wave propagation speeds because we are not considering the potential influence on damping from flows along the loop. The presence of flows leads to modification of the TGV relation (Soler et al. 2011).² Furthermore, studies of the interaction between the flows (namely the solar wind) and Alfvén waves suggest wave action conservation is important, which can result in dissipation-less waves undergoing apparent damping (Jacques 1977; Heinemann & Olbert 1980; McKenzie 1994; Cranmer et al. 2007; Li & Li 2007; Chandran et al. 2015). In the case of coronal loops, estimating flows is not a trivial endeavor; although the corona is likely to be in a state of thermal non-equilibrium and flows are expected to be present throughout. Several studies have tried to quantify the flow speeds, largely in active region loops, which are typically of the order of 10–50 km s⁻¹ (Reale 2010). Moreover, speeds of 74–123 km s⁻¹ have also been found in a single event (Ofman & Wang 2008). These studies suggest the axial flow speed is potentially small compared to the local Alfvén speed, and thus we expect this would have little effect on our results. However, further examination of flows in coronal loops is clearly required to assess their impact.

The increase in loop length between the wave paths corresponds to loops reaching higher altitudes in the corona. In Figure 5, we show the measured values of ξ as a function of loop length. Our measurements suggest that for the longer loops that reach higher up in the corona, the quality factors increases and, hence, the damping length increases. This suggests that the Alfvénic waves are subject to a reduced rate of damping. This is in contrast to the k–ω diagrams and power ratios, which show a greater difference between the outward and inward wave power. This is naturally explained by the fact that the inward waves have propagated further along the longer loops and have been damped to a greater degree than those in the shorter loops, despite the apparent reduced rate of damping in the longer loops. This result is present in both the MLE and least-squares fitting; however, the least-squares approach tends to overestimate the fitted values of power ratio and equilibrium parameter.

Zoom In Zoom Out Reset image size

Given the aforementioned problems with this dataset, related to identifying whether the selected wave paths are truly half the loop length, we are cautious in our interpretation of this variation in quality factor with loop length. A physical explanation for the decrease in damping rate can be made in terms of the density ratio between the internal and external plasmas. If we assume that the coronal loops are subject to similar rates of heating, and the rate of chromospheric evaporation is similar, then the average density of the longer loops is likely to be less than those of shorter loops. Hence, compared to the ambient plasma the density ratio (ρ_i/ρ_e) for longer loops is, on average, less than for the shorter loops. Equation (2) then implies that the equilibrium parameter will increase as the density ratio decreases, which matches the observed behavior.

Moreover, the fact that we have potentially not measured the wave path along half a loop will change the model that should be fitted to the power ratio (Equation (22)) and this will alter the measured values of the parameters. Considering the magnetic field extrapolation, there is the possibility that we have measured a loop segment that is less than half the loop length. In such a scenario, the average power over this shorter segment, compared to a half-loop segment, will be greater for outward waves (because the wave amplitudes averaged over have been damped less over this distance) and less for the inward waves (because the wave amplitudes averaged over will have been subject to greater damping). Hence, the power ratio will be artificially enhanced, giving the appearance of greater damping. This would lead to an underestimate of ξ compared to its true value. Hence, the observed effect of increasing ξ with height would be more pronounced.

Finally, it is also worth commenting on the measured value of the factor P_out/P_in, which represents the power of the waves input into the corona at each footpoint of the loop. The power ratio that we have obtained is almost equal to unity in the case of least-squares estimation, which is consistent to the previous study of Verth et al. (2010) and was interpreted as the same wave power being generated at both of the footpoints. In the case of MLE estimation, we obtain that the power ratio is less than unity, which implies that the wave power generated at the footpoint associated with inward waves is larger than the other. The current level of uncertainties associated with our measurements does not permit us to rule out that the input power is equal at both footpoints. However, it would not be surprising if the magnitude of the wave power is different at both footpoints, given that the physical conditions at the wave source region are likely to be dissimilar.