Proper Orthogonal and Dynamic Mode Decomposition of Sunspot Data

High resolution solar observations show the complex structure of the magnetohydrodynamic (MHD) wave motion. We apply the techniques of POD and DMD to identify the dominant MHD wave modes in a sunspot using the intensity time series. The POD technique was used to find modes that are spatially orthogonal, whereas the DMD technique identifies temporal orthogonality. Here we show that the combined POD and DMD approaches can successfully identify both sausage and kink modes in a sunspot umbra with an approximately circular cross-sectional shape.


Introduction
Analysis of oscillations in sunspot data began in the late 1960's, see e.g., Beckers and Tallant [1]. These authors determined the observational parameters of umbral flashes, a phenomenon that shows oscillatory behaviour in a sunspot. In the early 1970's several studies looked at Doppler velocity oscillations in sunspots. Bhatnagar [2] determined Doppler velocity oscillations with a period of the order of 180 − 220s. Later Beckers and Schultz [3] observed a peak period of around 180s. Measuring intensity oscillations directly from time lapse filtergram movies, Bhatnagar and Tanaka [4] detected periodicities of the order of 170 ± 40s. Later on, Moore [5] detected Doppler velocity oscillations of 120 − 180s and 240 − 300s in umbral and penumbral regions, respectively To the present day the study of oscillations in sunspots has mainly been carried out by Fourier transforming data to provide, e.g. power spectra, either on a pixel by pixel basis or integrating across a particular Region of Interest (ROI). Although such analysis can provide valuable information, for the identification of coherent structures, e.g. magnetohydrodynamic (MHD) wave modes, in the temporal and spatial domain across a particular ROI, the basic Fourier transform approach has its limitations. Despite this, one can fine tune a Fourier filter in the spatial and temporal domains to try and identify particular MHD wave modes, as was presented by Jess et al. [6] (hereafter J17) in order to detect a slow kink body mode in a sunspot umbra. In the present work we aim to apply the more advanced techniques of Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD) to identify low order MHD wave modes as coherent oscillations across the sunspot umbra, both in the spatial and temporal domains, using the same sunspot data as [6][7][8]. In the more general solar context, POD has previously been applied to decompose the Doppler velocity of the entire solar disc [9][10][11] and numerical convection data [12].

Observations
The dataset we will analyse has been previously used for studies of running penumbral waves [7], connections between photospheric and coronal magnetic fields [8] and in the detection of an umbral kink mode [6]. The portion of the complete multi-wavelength dataset used in the present study consists of a 75-minute observing sequence of Hα images acquired by the Hydrogen-Alpha Rapid Dynamics camera (HARDcam; [13]). The Hα time series, which observed the approximately circular sunspot present within the active region NOAA 11366, were obtained during excellent seeing conditions between 16:10 -17:25 UT on 10 December 2011, with the Dunn Solar Telescope (DST) at Sacramento Peak, New Mexico. The sunspot under investigation was located at N17.9W22.5 in the conventional heliographic co-ordinate system, or (356 , 305 ) in heliocentric co-ordinates. The filter employed had a full-width at half-maximum of 0.25 Å, which was centered on the Hα line core at 6562.8 Å. A platescale of 0.138 per pixel was used to provide a field-of-view size equal to 71 × 71 . On-site high-order adaptive optics [14], post-facto (speckle) image reconstruction techniques [15] and image destretching relative to simultaneous broadband continuum images [16] were implemented to improve the final data products, providing a cadence of 1.78 s. A sample Hα image of the sunspot is displayed in the left panel of Figure 1.

Modal decomposition techniques
Following the approach by Higham et al. [17], we are going to employ two techniques to identify low order MHD modes from the intensity time series. The POD method can be used to identify MHD wave modes by imposing the criteria that modes are spatially orthogonal. The second method, DMD, assumes a temporal orthogonality of modes. Hence, if observed MHD wave modes do not have identical frequencies, and this difference is resolved in the frequency domain, then DMD offers an optimal methodology to identify such modes.  The PSD shows peaks between frequencies 4.3 mHz and 6.5 mHz (corresponding to periods of 153 -232s).

(a) Proper Orthogonal Decomposition (POD)
The POD technique was developed by Pearson [18] as an analogue of the principal axis theorem in mechanics. POD was introduced as a mathematical technique in fluid dynamics by Lumley [19] to identify coherent structures in turbulent flow-fields. In the literature POD takes a variety of names depending on the field of application, such as principal component analysis (PCA) and Hotelling analysis. Since POD will produce as many modes as there are time snaphots in a dataset, the challenging part of this type of analysis is to identify which of the POD modes actually have a physical meaning. Hence, for identification of MHD wave modes in the umbral regions of sunspots care should be taken to compare POD modes with what we should expect from theoretical models, e.g. the MHD wave modes of cylindrical magnetic flux tubes. Let us consider the sequence of ROI intensity snapshots of the sunspot of spatial size X × Y and a time domain of size T . Each of these snapshots can be reorganized in a column matrix W ∈ R N ×T , where N = XY and N T , where each column of W will be defined as w i with i = 1...T such that There are several approaches to applying POD to a dataset. The classical POD method [20] is performed by computing the eigenvalues and the eigenvectors of the covariance matrix of the dataset.
Another approach is to obtain the POD of W using the optimum low rank approximation [20] and this is known as the Singular Value Decomposition (SVD). Applying the SVD, we obtain This decomposition gives the spatial structure of each mode in the columns of the matrix Φ ∈ R N ×T , i.e. φ i with i = 1...T and these modes are orthogonal to each other. The temporal evolution of the POD modes are given by the columns of the matrix C ∈ R T ×T . The particular spatial and temporal output of the POD presented here is the product of the N two dimensional spatially orthogonal eigenfunctions with their associated one dimensional time coefficients. Since POD places no restriction on the time coefficients, these can be periodic or aperiodic and the amplitude can also vary with time. The matrix S ∈ R T ×T is a diagonal matrix, and the modes are generally ranked according to their contribution to the total variance of the snapshot series. This contribution is given by the diagonal elements of matrix, λ, by means of the vector λ = diag(S) 2 /(N − 1).

(b) Dynamic Mode Decomposition (DMD)
The DMD technique, first introduced by Schmid [21], is a data-driven algorithm that can extract the dynamic information of the flow generated by numerical simulations or in a measured physical experiment [22]. DMD modes represent the spatial structure of the mode where the associated eigenvalues give information about the oscillation frequencies of the modes. DMD is a widely used technique in the field of fluid mechanics, e.g., jet flows [23,24], bluff body flows [25] and visco-elastic fluid flows [26]). It can therefore also extract information about the coherent spatial structure of observed MHD wave modes if the modes have distinct frequencies as we show in this study.
To apply the DMD technique the time snapshots have to be organized in columns analogously to POD, but in two matrices defined as where W B is shifted by a snapshot of W A such that τ = (T − 1). The matrices W A and W B are related by a linear operator A ∈ C N ×N as DMD is based on approximating the eigenvalues and eigenvectors of the linear operator A without actually computing them exactly since for most practical applications the size of A is too large. Using SVD the matrix W A is decomposed as   From this we defineÃ =Φ * AΦ, whereÃ ∈ C τ ×τ is the optimal low-dimensional representation of A, (note that τ N ) so that we can calculate the complex eigenvalues, µ i , and associated eigenvectors, z i , ofÃ, where i = 1...τ . To obtain the spatial structure of the DMD modes, we follow the method developed by Jovanovic et al. [24] by calculating a Vandermonde matrix, where i = 1...τ and j = 1...τ . After this operation is completed, the spatial structure of the DMD modes are obtained from 9) and the distinct frequencies associated with each these modes are where fs is the sampling frequency.

Method, results and MHD wave mode identification
Our goal is to use POD and DMD in combination to identify coherent oscillations across the sunspot's umbra and compare these modes with the MHD wave modes of a cylindrical magnetic flux tube predicted from theory. The particular ROI of the sunspot umbra to be studied is shown by the green box in Figure  1. Firstly, this ROI is analysed using the POD technique, which ranks modes based on their contribution to the overall variance. This step is followed by the calculation of the power spectral density (PSD) of the POD time coefficient associated with each of these modes. The PSD of the first 20 modes, which contains the majority of the energy (96.86 %), show frequency peaks between 4.3 mHz and 6.5 mHz as shown in the right panel of Figure 1. The PSD of the individual POD modes are then used to determine the dominant frequency, or frequencies if there are a mix of frequencies, associated with a particular POD mode, so that this information could be applied to determine the coherent spatial structure of modes with distinct frequencies using DMD. If there is no exact match between frequencies, the DMD mode closest to the target frequency is selected.
For illustrative purposes we will concentrate on the first branches of the sausage and kink slow body modes, i.e. modes with only one radial node occurring at the umbra/penumbra boundary. The first POD mode shows the clear azimuthal symmetry of a sausage mode as shown in the first column in Figure 2, with a PSD peak at 4.9 mHz as shown in the left panel of Figure 4. The DMD mode that corresponds to the same frequency of 4.8 mHz is shown in the second column in Figure 2. The third column shows the density perturbation of the slow body sausage mode from the cylindrical magnetic flux tube model. This is important for comparison since the MHD wave modes in a cylindrical flux tube are, by definition, spatially orthogonal. Since the observed umbra is approximately circular, POD, which defines modes by spatial orthogonality, should perform well in this particular case study. What is more remarkable is that the DMD technique, which does not have any such criteria, still manages to identify the sausage mode. From both the POD and DMD analysis there is strong oscillatory power in the penumbra at 4.8 mHz and the penumbral filaments can cleary be identified. Obviously, the idealised cylindrical magnetic flux tube model cannot recreate this oscillatory behaviour since it assumes a simple quiescent environment without complex fibril structuring. In addition, it is important to state that even within the umbra, disagreement between observations and the eigenmodes of a magnetic cylinder could simply be due to the fact that the observed oscillations are being continually forced and are not free.
The next POD mode that can be interpreted as a physical MHD wave mode is the 13 th mode which has the azimuthal asymmetry of a kink mode as shown in the first column of Figure 3, with a peak at 6 mHz as shown in the right panel of Figure 4. The DMD mode with frequency of 6 mHz is shown in the second column in the kink mode frequency from POD and DMD is certainly in the same frequency range as the filter applied by J17. Our analysis reveals that the time coefficients of the POD modes are almost sinusoidal. This is remarkable since POD puts no such condition on these coefficients. Hence, Fourier analysis, which has a sinusoidal basis in the temporal domain, in retrospect was a valid approach. The problem with Fourier analysis is the assumption of a sinusoidal basis in the spatial domain, since in the cylinder model, the basis functions in the radial direction are Bessel functions. The strength of POD is that it calculates a spatially orthogonal basis, regardless of the geometry of the observed waveguide. Also, the further advantage of both POD and DMD over Fourier analysis is that these methods cross-correlate individual pixels in the ROI, in the spatial and temporal domain, respectively. This is a distinct advantage in identifying a coherent oscillations across the whole umbra. In agreement with the sausage mode identification in Figure 2, the spatial structure of the POD and DMD modes in the first and second columns of Figure 3 is very similar even though the DMD places no restriction on the mode being orthogonal. This further strengthens the argument that the kink mode interpretation is indeed physical.
Here we would like to investigate the apparent rotational motion of the kink mode detected by J17 who constructed a time-azimuth diagram around the circumference of the umbra and estimated an angular velocity of approximately 2 deg s −1 and a periodicity of about 170 s. Physically, the rotational motion could be explained by having either (i) a kink mode that is standing in the radial direction but propagating in the azimuthal direction or (ii) it could be the superpostion of two approximately perpendicular kink modes (both standing in the radial and azimuthal directions). Before attempting to recover this rotational motion with the POD and DMD techniques it should be emphasised that the filtering process performed by J17 crudely oversimplified the complexity of the swirling "washing machine" motion in the original signal. In particular, the 40 s wide temporal filter could contain at least least 7 DMD modes. Spatially, the filter effectively divided the umbra into quadrants. To recreate the apparent rotational motion (or approximate circular polarisation) with POD we need to superimpose at least two spatially perpendicular kink modes with similar, but not necessarily identical, periods. From our analysis this requires the superposition of POD 10 (shown on the left panel of Figure 5 and POD 13 shown on the first column on Figure 3). The PSD of POD 10 has a peak at 5.4 mHz as shown on the left panel of Figure 6, while PSD of POD 13 has a peak at 6 mHz as shown on the right panel of Figure 4. Both these frequencies lie within the temporal filter chosen by J17. We can also recreate this rotational motion by superimposing at least two DMD modes. Although DMD modes are not defined to be orthogonal in space, we still find two examples of kink modes with DMD that are approximately perpendicular to each other and are also in the same frequency range of J17. These modes correspond to a frequency of 5.4 mHz (see the right panel of Figure 5) and 6 mHz shown on the second column of Figure 3. A similar time-azimuth analysis to J17 was performed on the superposition of these two DMD modes along the solid black circle shown on the right panel of Figure 5 where the signal was strongest. This resulted in an angular velocity of about 2 deg s −1 and periodicity of approximately 170 s (see the right panel of Figure 6), consistent with the result of J17.
To compare the cylinder model MHD modes with the POD modes from the observational data, we also performed a Pearson correlation analysis, calculated on a pixel-by-pixel basis for the sausage (see Figure 2) and kink (see Figure 3) modes using as shown on Figure 7. The result of the correlation is a number between 1 and -1, where 1 means that the pixels have a linear correlation while -1 denotes a linear anti-correlation. Certainly, there is a better correlation for the sausage than the kink, but this is not surprising since it is clearly visible from Figure 3 that signal for the kink mode is weaker than for sausage mode (see Figure 2). However, the kink mode stills shows a good correlation in the regions where amplitude is maximum (see Figure 7).

Conclusions
All the methods used to identify coherent oscillations across sunspots and pores have their particular strengths and weaknesses. We have demonstrated here that a more considered and multi-faceted approach can be more robust in pinpointing modes that are actually physical. For    example, the previous analysis by J17 required fine tuning the Fourier filters in the temporal and spatial domain to reveal the umbral kink mode confirmed by our POD and DMD analysis. In contrast, POD requires no such filtering, and indeed, such filtering would completely skew the results. The inherent problem with POD is identifying which modes are physical since this method produces as many modes as there are time snapshots. This is where further analysis is required as demonstrated in this study and previously by Higham et al. [17]. By calculating the PSD of each POD mode the dominant frequency (or frequencies) of each mode can be identified and these can be paired with the unique frequencies associated each DMD mode allowing for comparison between the spatial structure of the modes produced by both methods. If there is agreement between the spatial structure of both the POD and DMD modes (up to some specified accuracy), then this provides compelling evidence that the mode is indeed physical. To our knowledge, this is the first time the combined approach of using POD and DMD has been used on sunspot data to identify more than one MHD wave mode. We, therefore, suggest that in combination, POD and DMD could prove to be indispensable tools for decomposing the many possible MHD wave modes that could be excited in sunspots and pores, especially with the advent of high resolution  and its public directory can be found here https://www.nso.edu/data/historical-archive/. ROSA is a sixcamera high-cadence solar imaging instrument developed by Queen's University Belfast. Although Queen's University currently do not have the facilities to host such large data sets publicly for download, the ROSA data archive and ROSA reconstructed data are documented and can be requested and transferred (e.g. FTP or on a physical disk drive) from here https://star.pst.qub.ac.uk/wiki/doku.php/public/research-areas/solarphysics/rosa-archive.