Normal modes of transverse coronal loop oscillations from numerical simulations: I. Method and test case

TL;DR

This paper introduces an iterative numerical method combining time-dependent simulations and CEOF analysis to accurately determine the normal modes of coronal loops, validated against an analytical test case.

Contribution

The paper presents a novel iterative approach that integrates numerical simulations and CEOF analysis to extract coronal loop normal modes, improving accuracy over previous methods.

Findings

Method achieves less than 0.7% error in frequency and eigenfunctions after 6 iterations.

Eigenfunctions with discontinuities show larger errors at discontinuity positions.

Validated approach with an analytical test case before applying to real data.

Abstract

The purpose of this work is to develop a procedure to obtain the normal modes of a coronal loop from time-dependent numerical simulations with the aim of better understanding observed transverse loop oscillations. To achieve this goal, in this paper we present a new method and test its performance with a problem for which the normal modes can be computed analytically. In a follow-up paper, the application to the simulations of \citet{rial2013} is tackled. The method proceeds iteratively and at each step consists of (i) a time-dependent numerical simulation followed by (ii) the Complex Empirical Orthogonal Function (CEOF) analysis of the simulation results. The CEOF analysis provides an approximation to the normal mode eigenfunctions that can be used to set up the initial conditions for the numerical simulation of the following iteration, in which an improved normal mode approximation is obtained. The iterative process is stopped once the global difference between successive approximate eigenfunctions is below a prescribed threshold. The equilibrium used in this paper contains material discontinuities that result in one eigenfunction with a jump across these discontinuities and two eigenfunctions whose normal derivatives are discontinuous there. After 6 iterations, the approximation to the frequency and eigenfunctions are accurate to $\lesssim 0.7$\% except for the eigenfunction with discontinuities, which displays a much larger error at these positions.

Figures

Equations

$${\rho }_0(x)=\{\begin{array}{ll}{\rho }_i,& |x|\le a,\\ {\rho }_e,& |x|>a.\end{array}$$

$$v_A(x)=\{\begin{array}{ll}{v}_{Ai}\equiv \frac{B_0}{\sqrt{\mu {\rho }_i}},& |x|\le a,\\ v_{Ae}\equiv \frac{B_0}{\sqrt{\mu {\rho }_e}},& |x|>a,\end{array}$$

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot (\rho \mathbf{V}),$$

$$\rho \frac{\partial \mathbf{V}}{\partial t}=-\rho (\mathbf{V}\cdot \nabla )\mathbf{V}+\frac{1}{\mu }\left(\nabla \times \mathbf{B}\right)\times \mathbf{B},$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{V}\times \mathbf{B}\right).$$

$${\rho }_0\frac{\partial \mathbf{v}}{\partial t}=\frac{1}{\mu }\left(\nabla \times \mathbf{b}\right)\times {\mathbf{B}}_0,$$

$$\frac{\partial \mathbf{b}}{\partial t}=\nabla \times \left(\mathbf{v}\times {\mathbf{B}}_0\right),$$

$$\frac{\partial v_x}{\partial t}=\frac{B_0}{\mu {\rho }_0}\left(\frac{\partial b_x}{\partial z}-\frac{\partial b_z}{\partial x}\right),$$

$$\frac{\partial v_y}{\partial t}=\frac{B_0}{\mu {\rho }_0}\left(\frac{\partial b_y}{\partial z}+i{k}_y{b}_z\right),$$

$$\frac{\partial b_x}{\partial t}=B_0\frac{\partial v_x}{\partial z},$$

$$\frac{\partial b_y}{\partial t}=B_0\frac{\partial v_y}{\partial z},$$

$$\frac{\partial b_z}{\partial t}=-B_0\left(\frac{\partial v_x}{\partial x}-i{k}_y{v}_y\right).$$

$$v_x(x,z,t)$$

$$b_x(x,z,t)$$

$$b_z(x,z,t)$$

$$\omega {\hat{v}}_x=-\frac{B_0}{\mu {\rho }_0}\left[k_z(i{\hat{b}}_x)-\frac{d}{dx}(i{\hat{b}}_z)\right],$$

$$\omega (i{\hat{v}}_y)=\frac{B_0}{\mu {\rho }_0}\left[k_z{\hat{b}}_y+k_y(i{\hat{b}}_z)\right],$$

$$\omega (i{\hat{b}}_x)=-B_0{k}_z{\hat{v}}_x,$$

$$\omega {\hat{b}}_y=B_0{k}_z(i{\hat{v}}_y),$$

$$\omega (i{\hat{b}}_z)=-B_0\left[\frac{d{\hat{v}}_x}{dx}-k_y(i{\hat{v}}_y)\right].$$

$$\frac{d^2{\hat{v}}_x}{d{x}^2}=m^2{\hat{v}}_x,$$

$$m^2=k_y^2+k_z^2-\frac{{\omega }^2}{v_A^2}.$$

$$i{\hat{v}}_y=\frac{k_y}{m^2}\frac{d{\hat{v}}_x}{dx},$$

$$\frac{i{\hat{b}}_z}{B_0}=-\frac{1}{\omega }\frac{{\kappa }^2}{m^2}\frac{d{\hat{v}}_x}{dx},$$

$${\kappa }^2=k_z^2-\frac{{\omega }^2}{v_A^2}\equiv m^2-k_y^2.$$

$${\hat{v}}_x(x)=\{\begin{array}{ll}C\exp (m_{e}x),& \text{for }x<-a,\\ A\cosh (m_{i}x),& \text{for }-a\le x\le a,\\ C\exp (-m_{e}x),& \text{for }x>a,\end{array}$$

$$C=A\exp (m_{e}a)\cosh (m_{i}a).$$

$$\tanh (m_{i}a)=-\frac{{\kappa }_e^2}{{\kappa }_i^2}\frac{m_i}{m_e}.$$

$$v_x(x,z,t=0)=v_0\exp \left(-\frac{x^2}{a^2}\right)\exp \left(-\frac{z^2}{a^2}\right),$$

$$Re\left\{R(t_l)\exp [i\varphi (t_l)]S(x_i,z_k)\exp [-i\theta (x_i,z_k)]\right\},$$

$${\tilde{v}}_x(x_i,z_k)=S_{v_x}(x_i,z_k)\cos {\theta }_{v_x}(x_i,z_k),$$

Full text

Sun: oscillations — methods: numerical — techniques: miscellaneous

1 INTRODUCTION

The solar atmosphere is the site of a diversity of magnetohydrodynamic waves and oscillations. Transverse coronal loop oscillations are a prominent example of such events. They take place when a large energy deposition, usually caused by a flare, perturbs an active region magnetic structure, which sets some loops into oscillation (see, e.g., Aschwanden et al., 1999; Nakariakov et al., 1999, for some early observations). These events have been modeled with the help of slab and straight cylindrical loop models, whose normal modes can often be obtained by either analytical or numerical means (see Ruderman & Erdélyi, 2009, for a review). Starting with the simplest model that considers the fundamental transverse oscillation of a magnetic flux tube (Roberts, 1981; Edwin & Roberts, 1983; Nakariakov & Verwichte, 2005) several model improvements have included other effects, such as the curvature of coronal loops (van Doorsselaere et al., 2004, 2009; Terradas et al., 2006), longitudinal density stratification (Andries et al., 2005b, a), magnetic field expansion (Ruderman et al., 2008), departure from circular cross section of the tubes (Ruderman, 2003), or coronal loop cooling (Aschwanden & Terradas, 2008; Morton & Erdélyi, 2009). These ingredients have been seen to produce effects on the main wave properties, such as shifts on the frequency and position of the antinodes of the eigenfunctions. Also, the presence of internal fine structuring (Terradas et al., 2008) and/or a continuous cross-field inhomogenity in density is known to produce important effects, making possible physical processes such as phase-mixing (Heyvaerts & Priest, 1983) and resonant damping (Hollweg & Yang, 1988; Ruderman & Roberts, 2002; Goossens et al., 2002). The more general the model, the more difficult it is to calculate the eigenmodes of the structure and one has to resort to the study of time-dependent numerical simulations to study these transverse oscillations (Selwa et al., 2006, 2007, 2011a, 2011b; Rial et al., 2013). However, the comparison between the obtained numerical results to observed properties is not as straightforward as with the use of simple models. In these simulations, the initial disturbance excites different oscillatory harmonics, whose presence in the results is easily detected by a Fourier analysis of the variables collected at different points in the numerical domain, but this does not give information about the spatial structure of the eigenmodes. Hence, direct comparison between observed wave properties and the possibly present normal modes becomes difficult. For this reason, we have decided to devise the algorithm described in this paper, which allows us to isolate the eigenmodes present in a numerical simulation. Given the space required to present the algorithm, its application to the time-dependent numerical simulations of Rial et al. (2013), who use a model that takes into account effects such as density stratification, curvature, etc., is left for the second part of this work (Rial et al., 2019).

Normal modes provide a physical basis to understand the dynamics of a system. When the equilibrium configuration does not allow a simple solution of the normal mode problem, numerical techniques must be used to determine the normal modes’ eigenfunctions and eigenfrequencies. However, general purpose (i.e., for arbitrary equilibria) numerical codes that provide this information cannot be readily found. On the other hand, general purpose numerical codes to solve time-dependent equations are much more abundant. For this reason, being able to determine the normal modes of a system from time-dependent numerical simulations is a practical effort. A spectral analysis of the variables at different points in the spatial domain do give a good indication of the frequencies present in the results, but the very relevant spatial structure of the associated eigenmodes cannot be easily achieved with such analysis. Hence, a means of extracting the spatial profile of eigenfunctions together with their associated oscillatory frequencies from time-dependent simulations is desirable. In this way, the results can be compared to observations to ascertain the presence of a given normal mode in the coronal structure under study. The Complex Empirical Orthogonal Function (CEOF; see von Storch & Zwiers, 1999; Hannachi et al., 2007) analysis is a tool that satisfies these requirements: it takes as its input the numerical values of one or more variables over a spatial domain and for a given time span, and returns the spatial and temporal information about the main modes of variability contained in the data, which in our case will be the main eigenmodes present in the time-dependent numerical simulations. Thus, the aim of retrieving the normal mode features is feasible with this procedure.

The main advance of this paper is the repeated application of the described combination of time-dependent numerical simulations and CEOF analysis. The later results allow to determine initial conditions (for the numerical simulations) that more accurately resemble those of the normal mode, resulting in a numerical simulation in which the amplitude of all other normal modes is reduced with respect to the previous iteration. Therefore, a repetition of this process leads to successively better approximations to a normal mode and convergence to a prescribed accuracy can be achieved. Since our aim is to test the feasibility of the new method, we keep our model as simple as possible, considering a slab loop model and neglecting the model improvements mentioned above (coronal loop curvature, longitudinal density stratification, magnetic field expansion, …) In the presentation of the iterative method we follow a textbook approach: a simple test case with known solution is used, approximate solutions are found, the evolution of the error with the iterations is studied, and a proxy for this error that can be used in the stopping criterion is defined in terms of two successive approximations to the solution.

The outline of this paper is as follows: the equilibrium configuration and the equations for small amplitude perturbations are presented in Section 2. Analytical expressions for the normal modes of this system are introduced in Section 3. The time-dependent equations are solved in Section 4 for a prescribed initial condition and the CEOF analysis is applied to the results of this simulation; hence, the first iteration is complete, which allows us to give an approximation to the normal mode eigenfunctions and eigenfrequency. We next apply repeatedly the last two steps in an iterative process that improves the accuracy of the normal mode approximation (Section 6). Our conclusions are finally discussed in Section 7.

2 EQUILIBRIUM AND ZERO-$\beta$ GOVERNING EQUATIONS

We here use the Cartesian coordinate system shown in Figure 1. The equilibrium is invariant in the $y$-direction and consists of a dense plasma slab of width $2a$ that extends between $x=-a$ and $x=a$ and is embedded in a rarer environment that fills the space $|x|>a$. The whole system is bounded by the two planes $z=\pm L/2$, with $L$ the slab length. In the equilibrium the magnetic field is uniform and points in the direction of the slab axis: ${\boldsymbol{B}}_{0}=B_{0}\hat{\boldsymbol{e}}_{z}$; in addition, the plasma is at rest. This configuration has been often used to study the oscillations of a coronal loop. The $x$- and $z$-coordinates represent the directions transverse and longitudinal to the loop, respectively.

The equilibrium density is expressed as

[TABLE]

The internal (i.e., inside the slab) and external Alfvén velocities are

[TABLE]

with $\mu$ the permeability of free space.

We next introduce perturbations whose evolution is described by the ideal MHD equations, that in the zero-$\beta$ limit (i.e., zero plasma pressure) and in the absence of gravity read (Priest, 2014)

[TABLE]

[TABLE]

[TABLE]

Here $\rho$, ${\boldsymbol{V}}$, and ${\boldsymbol{B}}$ are the total (equilibrium plus perturbed) density, velocity, and magnetic field. Assuming small amplitude perturbations, Equations (3)–(5) can be linearized. The density perturbation is only present in the first of these equations, so that it is a secondary quantity that can be obtained once the velocity ($\boldsymbol{v}$) and magnetic field ($\boldsymbol{b}$) perturbations are known. The linearized momentum and induction equations can be expressed as follows:

[TABLE]

[TABLE]

where $\boldsymbol{v}$ and $\boldsymbol{b}$ are both functions of position and time.

Now, perturbations are assumed to propagate in the $y$-direction with wavenumber $k_{y}$ and so the $y$-dependence of $\boldsymbol{v}$ and $\boldsymbol{b}$ is of the form $\exp(-ik_{y}y)$. The Cartesian components of Equations (6) and (7) then reduce to111The right-hand side of the $z$-component of Equation (6) is equal to zero and so it leads to $v_{z}=0$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The velocity and magnetic field perturbations in these expressions are $\boldsymbol{v}(x,z,t)=v_{x}(x,z,t)\hat{\boldsymbol{e}}_{x}+v_{y}(x,z,t)\hat{\boldsymbol{e}}_{y}$ and $\boldsymbol{b}(x,z,t)=b_{x}(x,z,t)\hat{\boldsymbol{e}}_{x}+b_{y}(x,z,t)\hat{\boldsymbol{e}}_{y}+b_{z}(x,z,t)\hat{\boldsymbol{e}}_{z}$.

In this paper we impose that the slab has a finite length, $L$, in the $z$-direction (Figure 1) and that its ends are line-tied, that is, that the velocity perturbations are zero there. Moreover, in what follows we use the parameter values $\rho_{i}/\rho_{e}=10$, $L=50a$, and $k_{y}a=0.5$. Dimensionless values are obtained with the help of the length $a$, the velocity $v_{Ai}$, and the time $\tau_{Ai}=a/v_{Ai}$.

3 NORMAL MODES

Given that the plasma properties are uniform along the slab, the $z$-dependence of $v_{x}$ and $v_{y}$ is $\cos(k_{z}z)$. Equations (8)–(12) then reveal that the $z$-dependence of $b_{z}$ is $\cos(k_{z}z)$, while that of $b_{x}$ and $b_{y}$ is $\sin(k_{z}z)$. To satisfy the boundary conditions at the slab ends, $k_{z}$ must be equal to $(n+1)\pi/L$, with $n=0$ for the longitudinally fundamental mode, $n=1$ for its first longitudinal overtone, etc. To study normal modes a temporal dependence of the form $\exp(i\omega t)$ is also imposed and so the perturbed velocity and magnetic field components are (the $y$-dependence is omitted)

[TABLE]

Equations (8)–(12) now reduce to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now, the problem is to compute the $x$-dependence of the eigenfunctions $\hat{v}_{x}$, $i\hat{v}_{y}$, $i\hat{b}_{x}$, $\hat{b}_{y}$, and $i\hat{b}_{z}$, which are all real, and the eigenvalue $\omega$.

It is straightforward to eliminate all variables in favor of $\hat{v}_{x}$, which leads to the following ordinary differential equation,

[TABLE]

with

[TABLE]

The parameter $m$ takes the value $m_{i,e}$ when the Alfvén speed is substituted by its value $v_{Ai,e}$ inside and outside the slab, respectively. After determining $\hat{v}_{x}$ one can obtain $i\hat{v}_{y}$ and $i\hat{b}_{z}$ from

[TABLE]

[TABLE]

where

[TABLE]

Again, $\kappa$ takes the values $\kappa_{i,e}$ inside and outside the slab, respectively. The eigenfunctions $i\hat{b}_{x}$ and $\hat{b}_{y}$ follow from Equations (18) and (19), and are just proportional to $\hat{v}_{x}$ and $i\hat{v}_{y}$, respectively.

To solve Equation (21) one must impose boundary conditions in the $x$-direction, together with the proper jump conditions at the $x=\pm a$ interfaces, which according to, e.g., Goedbloed & Poedts (2004) are the continuity of the normal velocity component ($\hat{v}_{x}$) and of the total pressure, which in turn leads to the continuity of $i\hat{b}_{z}$.

Because of the symmetry222The imposed boundary conditions in the $x$-direction are also symmetric: see Sections 3.1 and 3.2. of the equilibrium and of Equations (16)–(20) with respect to $x=0$, eigenfunctions are either even or odd: for kink modes $\hat{v}_{x}$ and $i\hat{b}_{x}$ are even about the slab axis, while $i\hat{v}_{y}$, $\hat{b}_{y}$, and $i\hat{b}_{z}$ are odd; for sausage modes, the parity of the 5 eigenfunctions is the opposite. In our simulations only kink solutions are excited and so we restrict our analysis to these normal modes.

3.1 Laterally evanescent normal modes

Arregui et al. (2007) solved the eigenproblem of Equations (16)–(20) for solutions that are laterally evanescent, that is, for which the perturbations vanish as $x\rightarrow\pm\infty$; see their Section 3 and also Roberts (1981) for the treatment of the $k_{y}=0$ case. The kink solution that satisfies these constraints has the following $x$-velocity component:

[TABLE]

where the positive value of $m_{e}$ is taken and

[TABLE]

The constant $A$ can be arbitrarily chosen and so we set $A=1$. The eigenfrequency is the solution to the dispersion relation

[TABLE]

Figure 2 displays the eigenfunctions $\hat{v}_{x}$, $i\hat{v}_{y}$, and $i\hat{b}_{z}$ for the fundamental kink mode. They possess the parity and continuity properties described above: $\hat{v}_{x}$ is even, $i\hat{v}_{y}$ and $i\hat{b}_{z}$ are odd, and $\hat{v}_{x}$ and $i\hat{b}_{z}$ are continuous at the interfaces $x=\pm a$. On the other hand, $i\hat{v}_{y}$ jumps abruptly at these boundaries. In addition, these functions decay exponentially with $x$, as described by Equation (26). The longitudinal harmonics have a similar spatial structure of eigenfunctions. The frequencies for the $n=0,2,4$ longitudinal harmonics are $\omega/\tau_{Ai}=0.1011,0.2989,0.4852$.

3.2 Laterally confined normal modes

In Section 4 we solve numerically the initial value problem made of Equations (8)–(12) with suitable initial and boundary conditions. We consider the spatial domain $-20a\leq x\leq 20a$. Given that the boundaries are sufficiently far from the slab, the evanescent eigensolution of Figure 2 is almost zero at the edges of the numerical domain and so it is, in practice, a solution to the initial value problem. Placing the boundaries at a finite distance from the slab, however, adds new, non-evanescent eigensolutions that can be excited by the initial perturbation. By replicating the analysis of Section 3.1 with the boundary condition $\hat{v}_{x}=0$ at $x=\pm L_{x}$ the laterally confined eigenfunctions can be obtained. It will suffice to say that the fundamental confined mode has $\omega=1.676/\tau_{Ai}$.

4 TIME-DEPENDENT NUMERICAL SIMULATIONS

4.1 Simulation setup and numerical method

We now solve numerically Equations (8)–(12) in the region $-L_{x}\leq x\leq L_{x}$, $-L/2\leq z\leq L/2$; see Figure 1. The coefficients of the system of partial differential equations can be made real by using the independent variables $v_{x}$, $iv_{y}$, $b_{x}$, $ib_{y}$, and $b_{z}$. Our initial disturbance is such that the full slab is subject to an initial transverse forcing given by

[TABLE]

while all other variables are initially zero. This initial perturbation represents a sudden deposition of energy at the slab center. The $v_{x}$ perturbation is even about $x=0$ and so can only excite kink modes. Since the transverse profile of $v_{x}$ resembles that of the laterally evanescent kink modes (top panel of Figure 2), a large portion of the energy in the initial disturbance will go to these modes. But one can expect that the laterally confined mode of Section 3.2 will also be excited.

The simulation box is determined by the lengths $L_{x}=20a$ and $L=50a$, and a uniform grid of 4001$\times$51 points in the $x$- and $z$-directions is used. The grid is coarser along the slab because it is sufficient to capture well the smooth sinusoidal dependence of normal modes in the $z$-direction; on the other hand, the grid is much finer across the slab because normal modes have much more structure in this direction. The numerical simulation is stopped at $t\simeq 280\tau_{Ai}$, which is $\sim 4.5$ and $\sim 75$ periods of the fundamental laterally evanescent and confined modes, respectively. The time step is $\Delta t=0.704\tau_{Ai}$.

The numerical method used to solve the linearized wave equations is based on the method of lines (MOL). Time and space are treated independently, using a third order Runge-Kutta method and a six order finite difference method, respectively. Artificial dissipation is included to avoid oscillations on the grid scale. This method has been used successfully in the past (e.g., Bona et al., 2009) and has a weak effect on the attenuation of the physical oscillations reported in the simulations. Since the linear hyperbolic MHD equations are solved explicitly, the time step is subject to the CFL condition. Note that in the linearized MHD equations there are terms proportional to $k_{y}$: these terms are incorporated to the code as simple source terms.

Although we solve the linearized MHD equations, there are jumps in the perturbed variables (in $iv_{y}$ and $b_{y}$) due to the discontinuities in the equilibrium variables. We have decided to use a simple numerical scheme that is not shock-capturing (better suited for discontinuities) since the effect of the jump is rather small in the temporal evolution of the different quantities.

Line-tying conditions are applied at $z=\pm L/2$, meaning that the velocities are zero, while for the rest of variables the derivatives with respect to $z$ are zero. At $x=\pm L_{x}$ we impose that the derivatives with respect to $x$ of all the variables are zero. This condition does not allow a perfect outward transmission of the waves and some reflections are produced. A direct consequence of these reflections is the presence of the laterally confined normal mode in our simulations.

4.2 Results

The initial condition excites a large number of longitudinal harmonics, both evanescent and confined in the $x$-direction, together with leaky waves that travel away from the slab. The emission of these leaky waves is clear until $t\simeq 50\tau_{Ai}$, after which only the kink normal modes of Sections 3.1 and 3.2 remain. Evidence of the presence of these normal modes comes from the spectral analysis of $v_{x}$, $iv_{y}$, and $b_{z}$ at a given location, which is selected so that normal modes have a non-negligible amplitude. Thus, for the transverse velocity component, $v_{x}$, we choose the point $x=0$, $z=0$, while for the other two variables the position $x=a$, $z=0$ is preferred. The Lomb-Scargle periodograms (Lomb, 1976; Scargle, 1982) at these points are shown in Figure 3. The three panels display the largest power peaks at $\nu=0.01601/\tau_{Ai}$ (i.e., $\omega=0.1006/\tau_{Ai}$), $\nu=0.04721/\tau_{Ai}$ ($\omega=0.2966/\tau_{Ai}$), and $\nu=0.07681/\tau_{Ai}$ ($\omega=0.4826/\tau_{Ai}$); these values are in excellent agreement with those of the lowest three laterally evanescent harmonics. The periodograms also show large power above $\nu=0.2/\tau_{Ai}$ caused by the excitation of laterally confined normal modes. Indeed, the largest peak in this frequency range is at $\nu=0.2544/\tau_{Ai}$ (i.e., $\omega=1.598/\tau_{Ai}$), again in good agreement with the value quoted in Section 3.2. It is worth noting that the power at $\omega=0.1011/\tau_{Ai}$ is 2–3 orders of magnitude higher than that at $\omega=1.676/\tau_{Ai}$ for $v_{x}$ and $iv_{y}$, although in the case of $b_{z}$ the two peaks are of similar magnitude. The reason for this is that the height of a power peak comes from the combination of the energy deposited by the initial disturbance in each normal mode (which is much larger for the evanescent one) and the amplitude of each variable (which in the case of $b_{z}$ is much smaller for the evanescent normal mode than for the confined one). The combination of these two factors results in the function $b_{z}$ containing similar power in the evanescent and confined normal modes in this numerical simulation.

5 COMPLEX EMPIRICAL ORTHOGONAL FUNCTION ANALYSIS

In this section we go beyond the normal mode frequency we just obtained and attempt to determine the normal mode eigenfunction structure using the Complex Empirical Orthogonal Function (CEOF) analysis. We will here give a very brief summary of this method; a detailed description can be found in Horel (1984); Wallace & Dickinson (1972); von Storch & Zwiers (1999); Hannachi et al. (2007, where it is called Complex Hilbert EOF) and an application to the study of coronal oscillations in Terradas et al. (2004).

The CEOF analysis is a numerical method that takes as its input a field $U(\vec{r},t_{l})$ discretized over a spatial mesh of points $\vec{r}=(x_{i},y_{j},z_{k})$ and evaluated at the discrete times $t_{l}$. Its output is a set of CEOF modes, which are not necessarily associated to physical modes of the system under study, each of them described by four measures called the temporal amplitude and phase and the spatial amplitude and phase. Together with these measures, the CEOF analysis associates to each mode a fraction of the total field variance. Once the CEOF code is fed with the input field, the “highest contributing” CEOF mode, that is, the one associated to the largest fraction of the total field variance, is retrieved first and other CEOF modes are obtained next in decreasing order of their fraction of the total field variance. The sum of the modes’ fraction of the total field variance tends to 1 as the number of CEOF modes is increased. The execution is stopped when the percentage of the total field variance accounted for by all the retrieved CEOF modes exceeds a pre-established value, here taken as 99.9%.

In our case the field $U$ can be, for example, the velocity component $v_{x}(x_{i},z_{k},t_{l})$ obtained in the numerical simulation of Section 4. This means that the input field is a three-dimensional data cube. Such as mentioned above, each of the obtained CEOF modes has empirically computed temporal and spatial measures, called the temporal amplitude, $R(t_{l})$, the temporal phase, $\phi(t_{l})$, the spatial amplitude, $S(x_{i},z_{k})$, and the spatial phase, $\theta(x_{i},z_{k})$. The spatial and temporal variability of the field described by this CEOF mode is

[TABLE]

where $Re$ denotes the real part. A CEOF mode that, for example, represents a propagating wave has a temporal phase that varies linearly with $t_{l}$ and a spatial phase that varies linearly with $x_{i}$ and $z_{k}$. The CEOF modes that we expect to find when analyzing the results of the numerical simulation, however, are standing waves. In this case, the temporal phase also varies linearly with $t_{l}$, but the spatial phase is such that two regions in which the difference of $\theta(x_{i},z_{k})$ is an integer multiple of $2\pi$ correspond to in-phase oscillations, while oscillations that are in anti-phase display a phase difference that is an odd multiple of $\pi$. In our results we will also find that a standing wave can have a phase that slowly varies in space, which is nothing but a modulation of $S(x_{i},z_{k})$ by the factor $\exp[-i\theta(x_{i},z_{k})]$. Section 3 of Terradas et al. (2004) gives simple two-dimensional examples of the outcome of the CEOF analysis when applied to a synthetic signal made of the sum of a propagating and a standing wave.

Our hypothesis is that the CEOF analysis applied to the results of the numerical simulation of Section 4 will provide an approximation, by means of Equation (30), to the evanescent normal mode eigenfunctions. Given that the eigenfunctions do not depend on time, we will ignore the temporal variation given by the measures $R(t)$ and $\phi(t)$ in Equation (30) and will only retain the real part of the spatial measures. Let us take, for example, the variable $v_{x}$, which for a normal mode has the eigenfunction $\hat{v}_{x}(x)\cos(k_{z}z)$. The CEOF approximation to this eigenfunction is

[TABLE]

where $S_{v_{x}}$ and $\theta_{v_{x}}$ are the spatial amplitude and phase of $v_{x}$ and the tilde in $\tilde{v}_{x}$ indicates that this is an approximation to the normal mode $v_{x}$. A numerical comparison between the normal mode eigenfunction and its approximation from the CEOF analysis is obtained with the help of the difference

[TABLE]

An analogous expression can be built for all other eigenfunctions.

Regarding $iv_{y}$, its eigenfunction for the confined mode is $i\hat{v}_{y}(x)\cos(k_{z}z)$ with $i\hat{v}_{y}(x)$ given by Equations (23) and (26). The CEOF approximation to this eigenfunction is

[TABLE]

with $S_{iv_{y}}$ and $\theta_{iv_{y}}$ the spatial amplitude and phase of $iv_{y}$. The case of $b_{z}$ requires special attention. Its eigenfunction is $i\hat{b}_{z}(x)\cos(k_{z}z)$, where $i\hat{b}_{z}(x)$ can be obtained from Equations (24) and (26). In the numerical simulation, however, we have not used the variable $ib_{z}$ but $b_{z}$. For this reason, the CEOF approximation to $ib_{z}$ requires inserting a factor $i$ inside $Re\left\{\ldots\right\}$ of Equation (30). We then have that the CEOF approximation to $ib_{z}$ is

[TABLE]

where $S_{b_{z}}$ and $\theta_{b_{z}}$ are the spatial amplitude and phase of $b_{z}$. Concerning $b_{x}$ and $b_{y}$, Equations (18) and (19) tell us that their respective CEOF approximations can be obtained from those of $v_{x}$ and $iv_{y}$. Finally, approximations to $\hat{v}_{x}(x)$, $i\hat{v}_{y}(x)$, and $i\hat{b}_{z}(x)$ can be derived by taking a cut along $z=z_{k}$ of Equations (31), (33), and (34).

Before applying the CEOF method to the results of our simulation, two more comments are in order. First, the transient phase is excluded from the analysis by considering $t\geq 50\tau_{Ai}$ only. Second, to reduce the memory requirements and speed up the computation of the CEOF modes, the values of $v_{x}$, $iv_{y}$ and $b_{z}$ are interpolated from the 4001$\times$51 numerical grid to a grid of $N_{x}\times N_{z}$ points (here we use $N_{x}=201$, $N_{z}=25$). To do so, in the $x$- and $z$-directions only 1 every 20 points and 1 every 2 points, respectively, from the numerical simulation are kept for the CEOF analysis.

5.1 Results

The CEOF method has the possibility of analyzing several fields simultaneously, which allows a better characterization of the physical modes because more restrictions are imposed by the higher complexity of the combined fields. Thus, we run the CEOF code on the fields $v_{x}(x_{i},z_{k},t_{l})$, $iv_{y}(x_{i},z_{k},t_{l})$, and $b_{z}(x_{i},z_{k},t_{l})$ together. To do this, the three data cubes are put next to each other and a larger data cube is created. We choose to join the three 201$\times$25$\times$330 cubes by attaching their $xt$-faces, so that the CEOF input is a cube of 201$\times$75$\times$330 data values. After the CEOF analysis is complete, we obtain a collection of CEOF modes, each of them characterized by its temporal amplitude and phase, $R(t)$ and $\phi(t)$, and its spatial amplitude and phase, $S(x,z)$ and $\theta(x,z)$, that can be split into the spatial measures $S_{v_{x}}$ and $\theta_{v_{x}}$ of the field $v_{x}(x_{i},z_{k},t_{l})$, $S_{iv_{y}}$ and $\theta_{iv_{y}}$ of the field $iv_{y}(x_{i},z_{k},t_{l})$, and $S_{b_{z}}$ and $\theta_{b_{z}}$ of the field $b_{z}(x_{i},z_{k},t_{l})$. These measures can in turn be inserted into Equations (31), (33), (34) to obtain the approximate CEOF eigenfunctions.

The first CEOF mode accounts for 64.8% of the total field variance and corresponds to the fundamental evanescent normal mode. Its frequency is determined by fitting the straight line $\phi(t)=\omega t+\phi_{0}$ to the temporal phase, which yields $\omega=0.1016/\tau_{Ai}$. This value is in excellent agreement with the normal mode frequency $\omega=0.1011/\tau_{Ai}$. The goodness of the CEOF approximation to the normal mode can also be judged with the help of the differences $\Delta v_{x}$, $\Delta v_{y}$, $\Delta b_{z}$, which are presented in Figure 4. To make this figure, the $\hat{v}_{x}$ eigenfunction is normalized to a maximum value of 1 and the CEOF approximation $\tilde{v}_{x}$ is also normalized to 1 at the position where the eigenfunction is maximum. The conclusion from this figure is that the CEOF analysis of the numerical simulation results allows us to recover the normal mode eigenfunction $\hat{v}_{x}$ with an error below 4%. We next turn our attention to the error of $iv_{y}$ and $ib_{z}$. The middle row of Figure 4 gives the difference $\Delta v_{y}$. Except for the points on the boundaries, $x=\pm a$, the error is smaller than 10% inside the slab ($|x|<a$) and practically zero outside the slab ($|x|>a$). The bottom row of Figure 4 gives the difference $\Delta b_{z}$, which also attains its largest value, of the order of 15% the eigenfunction value, at the slab boundary.

The second CEOF mode accounts for 15.9% of the total field variance and corresponds to the second longitudinal evanescent overtone. A linear fit to the temporal phase results in the frequency $\omega=0.2969/\tau_{Ai}$, which is very close to the analytical value $\omega=0.2989/\tau_{Ai}$. Figure 5 shows the difference between the normal mode and CEOF approximation to the eigenfunctions. We see that the error in the second longitudinal overtone is almost a factor of 2 better than that of the fundamental mode. When comparing the exactness of the CEOF approximation for the fundamental and the second overtone, we find a better agreement in the second case because there are more periods of this normal mode in the numerical simulation.

The conclusion of this section is that, while the CEOF approximation to $\hat{v}_{x}$ is acceptable, those to $i\hat{v}_{y}$ and $i\hat{b}_{z}$ are not too good. This situation will be improved by the application of the iterative method presented in the following section.

6 ITERATIVE METHOD

6.1 Description of the method

The scheme we have used so far consists of two steps: (i) a time-dependent numerical simulation of Equations (8)–(12) followed by (ii) the CEOF analysis of the obtained results. If we imagine that an eigenmode is perfectly described by a CEOF mode, then one could run a numerical simulation with initial conditions given by the eigenfunctions and so the obtained temporal evolution would be that of the eigenmode. At this point, this is not the case, but we have seen that the CEOF analysis produces an approximation to a normal mode eigenfunctions. We thus devise an iterative method that is made of the repeated application of steps (i) and (ii), in which the initial conditions of the numerical simulation of a given iteration are taken from the CEOF method of the previous iteration. The iterations will be stopped once a given measure of goodness is reached. The iterative method is carried out separately for each normal mode.

We first need to determine which information is required from the CEOF analysis to fix the initial conditions. Rather than using the time dependence $\exp(i\omega t)$ of Equations (13)–(15) we assume that $v_{x}(x,z,t)$ is maximum at $t=0$ and so it has the form $v_{x}(x,z,t)=\hat{v}_{x}(x)\cos(k_{z}z)\cos(\omega t)$. Now, Equations (8)–(12) tell us that $v_{x}$ and $iv_{y}$ are in phase (in time) and that they are a quarter of a period out of phase with respect to $b_{x}$, $ib_{y}$, and $b_{z}$. This implies that $iv_{y}(x,z,t)$ is also maximum at $t=0$ and that the perturbed magnetic field components vanish at the start of the numerical simulation. Hence, the information that the CEOF analysis must provide to repeat step (i) is the approximation to $\hat{v}_{x}$ and $i\hat{v}_{y}$ provided by Equations (31) and (33).

6.2 Results

We are then ready to carry out the iterative process. Iteration #1 consists of the numerical simulation of Section 4 and the CEOF analysis of Section 5.1. The results we present now are a summary of the performance of 6 iterations, which are carried out independently for the fundamental evanescent mode and its second longitudinal overtone. An excellent assessment of the performance of the iterative method can be gained from the power spectra of $v_{x}$, $iv_{y}$, and $b_{z}$ for the numerical simulation of the last iteration. These power spectra are shown in Figure 6 for the two normal modes. Whereas the numerical simulation of iteration #1 displays power peaks at the frequencies of many harmonics, both simulations of iteration #6 only show a power peak for a single normal mode. The left (right) panels of Figure 6 have non-negligible power around the maximum $\nu=0.01601/\tau_{Ai}$ ($\nu=0.04721/\tau_{Ai}$), that are identical to those obtained from the power spectra in the first iteration. All other normal modes are virtually absent in the numerical simulations of the last iteration.

The iterative method yields another approximation to the frequency that comes from the CEOF analysis of the numerical simulation of iteration #6. A linear least-squares fit to the temporal phase gives the frequency $\omega=0.1016/\tau_{Ai}$ for the fundamental evanescent mode and $\omega=0.2969/\tau_{Ai}$ for the second longitudinal evanescent overtone. These values are identical to those of the first iteration. The error associated to these approximate values is 0.5% and 0.7%, so that the obtained accuracy is excellent.

6.3 Error and stopping criterion

We next examine in detail the error of the CEOF approximation to the eigenfunctions. Again we consider the results of iteration #6 and show these errors in Figures 7 and 8. Although their maximum values are reduced by a factor of 2, the errors display some of the patterns of iteration #1: the error of $v_{x}$ is maximum in the slab, the error of $iv_{y}$ has a dominant component at the slab boundaries, but has been strongly reduced inside the slab during the iterative process, and the error of $b_{z}$ is maximum at the slab boundaries, but has become much more confined to the slab neighborhood.

We finally analyze the evolution of the error with the iterations333Before computing the errors described here we normalyze the normal mode eigenfunctions and their approximation from the CEOF analysis so that $v_{x}$ equals zero at $x=0$, $z=0$.. At the end of iteration #$n$, with $n=1,2,\ldots$, Equations (31), (33), and (34) provide us with approximations for the three main eigenfunctions; we denote these approximations with the superscript $n$, i.e., $\tilde{v}_{x}^{n}$, $i\tilde{v}_{y}^{n}$, $i\tilde{b}_{z}^{n}$. For each eigenfunction, we define a global measure of the error, $\varepsilon$, by summing over the spatial domain the squares of the difference between the normal mode eigenfunction and the CEOF approximation. For example, for the iteration #$n$ and the variable $v_{x}$ this global error is:

[TABLE]

where the factor $\max_{i,k}{|\hat{v}_{x}(x_{i})\cos k_{z}z_{k}|}$ in the denominator provides the right normalization that enables us to compare the error of different eigenfunctions. The additional factors $N_{x}$ and $N_{z}$ give an additional normalization that removes the dependence of $\varepsilon_{v_{x}}^{n}$ on the number of points in the CEOF analysis. The definitions of $\varepsilon_{iv_{y}}^{n}$ and $\varepsilon_{ib_{z}}^{n}$ are done in a similar way.

The top panels of Figure 9 present the global errors for the first 6 iterations for the fundamental evanescent mode (left column) and its second longitudinal overtone (right column). In each iteration, $v_{x}$ has the smallest error (possibly because it is the eigenfunction with less “contamination” from the confined normal mode) and $iv_{y}$ displays the largest global error (because of the large contributions at the slab boundaries, that do not disappear with the iterative process). We see that the biggest improvement in the global error is obtained in iteration #2, for which a remarkable reduction in $\varepsilon_{v_{x}}$ and $\varepsilon_{ib_{z}}$ is found. The subsequent variation of the three global errors is much more moderate and so for this problem two iterations give a good compromise between the error associated to the CEOF approximations and the computer time spent.

The case studied in this paper allows us to compute the global error because of our knowledge of the exact eigenfunctions. In a general case, in which the eigenfunctions are unknown and our aim is just to obtain them, a proxy for the global error can be computed by comparing the approximate eigenfunctions of two successive iterations. To define this new global uncertainty measure, in Equation (35) we replace the normal mode eigenfunction $\hat{v}_{x}(x_{i})\cos k_{z}z_{k}$ by its CEOF approximation in the iteration $n+1$. We also rewrite the iteration indices and substitute $n+1$ by $n$. This gives the following definition for the global difference between the $v_{x}$ eigenfunction of iterations #$n$ and #$n-1$:

[TABLE]

Analogous expressions can be written for $\delta_{iv_{y}}$ and $\delta_{ib_{z}}$.

The variation of the global difference with the iterations is displayed in the bottom panels of Figure 9. The convergence is quite fast, with $\delta_{ib_{z}}$ showing an improvement of roughly two orders of magnitude per iteration during the first iterations, the highest convergence rate of all variables. $\delta_{v_{x}}$ and $\delta_{iv_{y}}$ are reduced at a slower pace, namely, almost an order of magnitude per iteration. We see that all variables attain a global difference smaller than $10^{-5}$ in iteration #5 (fundamental evanescent mode) and in iteration #6 (second longitudinal evanescent overtone). Therefore, we adopt the iteration stopping criterion that the global difference must be smaller than $10^{-5}$.

7 CONCLUSIONS

In this paper we have devised a method to determine a physical system normal modes, i.e., their eigenfunctions and eigenfrequencies, by the iterative application of time-dependent numerical simulations of the equations that govern the system dynamics and the CEOF analysis of the simulation results. We have illustrated how the CEOF method can be applied to all the non-redundant variables: in our case, in particular, this means that we can avoid including $b_{x}$ and $ib_{y}$ in the CEOF computation because their eigenfunctions can be readily computed from the other three ($v_{x}$, $iv_{y}$, and $b_{z}$). At the end of each iteration, the CEOF approximations to the eigenfunctions are used as the initial conditions for the time-dependent numerical simulation of the next iteration. Finally, we have examined the global error of the approximate eigenfunctions as a function of the iteration number and have established a convergence criterion based on the global difference between the approximate eigenfunctions of consecutive time steps.

The main disadvantage of our test case is the presence of sharp boundaries in the equilibrium structure, which leads to abrupt jumps of the eigenfunction $i\hat{v}_{y}$ and non-derivable eigenfunctions $\hat{v}_{x}$ and $i\hat{b}_{z}$ at these positions (see Figure 2). We have found that these normal mode features result in the presence of large errors at the slab boundary, which are quite substantial for the approximation to $i\hat{v}_{y}$; see Figures 7 and 8.

We have obtained an approximation to the two normal modes of interest (the fundamental evanescent mode and its second longitudinal overtone) with great accuracy: after 6 iterations, the frequency is wrong by only 0.5%–0.7% and the eigenfunctions $\hat{v}_{x}$ and $i\hat{b}_{z}$ have maximum errors of the order of 0.6% and 0.7%, respectively. The case of $i\hat{v}_{y}$ is worse because of the difficulties of recovering a function that jumps at $x=\pm a$. If these two lines are ignored, the maximum error of $i\hat{v}_{y}$ is also of the order of 0.6%.

In the second paper of this work (Rial et al., 2019) we will apply the technique presented here to the time-dependent numerical simulations of a loop embedded in a coronal arcade carried out by Rial et al. (2013). The equilibrium structure is similar to the one used in the present paper but includes a curved slab in which both the magnetic field strength and plasma density vary along the magnetic field. The initial condition used by Rial et al. (2013) is analogous to Equation (29) and so various longitudinal harmonics are excited. The present paper shows that our technique is well suited for this task because it allows to obtain the features of different normal modes. Its application to the more realistic numerical simulations by Rial et al. (2013) should produce normal mode characteristics comparable to observed loop oscillation events.

Acknowledgements

R.O. and J.T. acknowledge support from the Spanish Ministry of Economy and Competitiveness (MINECO) and FEDER funds through project AYA2017-85465-P. I.A. acknowledges support by the Spanish Ministry of Economy and Competitiveness (MINECO) through projects AYA2014-55456-P (Bayesian Analysis of the Solar Corona) and AYA2014-60476-P (Solar Magnetometry in the Era of Large Solar Telescopes) and FEDER funds. R.O. is indebted to D. W. Fanning for making available the Coyote Library of IDL programs (http://www.idlcoyote.com/).