Data Analysis in Space Physics

Diffusion Region Detection

Magnetic reconnection occurs in a small diffusion region around the X-line but its consequences are large scale. There are different approaches to determine whether or not we have encountered the diffusion region. The expected kinetic signatures of diffusion region:

• enhanced dissipation
• non-gyrotropic behavior
• electron distribution functions

The first proposed non-gyrotropy measure by Jack Scudder and William Daughton is

which only considers the non-gyrotropy in the plane perpendicular to the magnetic field.

Another measure of non-gyrotropy $D_{ng}$ of a specific population proposed by researchers can be think of a ratio between off-diagonal terms and diagonal terms of the pressure tensor:

Note that this is frame independent and population specific.

However, Marc Swisdak later demonstrated in his paper that a better measure from positive semi-definite matrix analysis should be

where $I_1=P_{xx}+P_{yy}+P_{zz}$ is the trace and $I_2=P_{xx}P_{yy}+P_{xx}P_{zz}+P_{yy}P_{zz}-(P_{xy}P_{yx}+P_{xz}P_{zx}+P_{yz}P_{zy})$ is the principle minor, both of which are invariants under coordinate rotation. (The third invariant is the determinant.) This indeed peak near the center of X-line while is enhanced along the separatrices, which is better than $D_{ng}$ in the sense that the latter may not peak at the diffusion region.

It has been found from PIC simulations that especially in asymmetric reconnection, $E_M^\ast\neq 0$ where $\mathbf{E}^\ast=\mathbf{E}+\mathbf{V}_e\times\mathbf{B}$ is not a good enough indicator. Following the principle ideas of finding a measure:

• dissipation should be related to non-ideal energy conversion;
• scalar indicator;
• insensitive to the relative motion,

a frame independent dissipation measure is constructed,

which is a Lorentz-invariant scalar representing the energy conversion rate in the frame of electron bulk motion. In the nonrelativistic limit, one can simplify the equation by setting $\gamma_e\rightarrow 1$. This can be easily confirmed by multiplying $\mathbf{j}^\prime=qn_i\mathbf{V}_i=\mathbf{j}-\rho_c\mathbf{V}_e$ and $\mathbf{E}^\prime = \mathbf{E}+\mathbf{V}\times\mathbf{B}$.

This is another good indicator of the electron diffusion region, although it is slightly enhanced along the separatrices as well.

An alternative and complementary approach is to identify diffusion region candidates by a combination of representative properties. In practice for reconnection with no guide field for instance, the overlapped region of $B_L=0$ and $E_N$ peak, where some researchers call 'shoulder', is already good enough for the identification.

Spectral Analysis

FFT

Periodogram

This is a group of techniques to determine the periodicity of data. Essentially it is the same as Fourier transform, but with the inverse of frequency as the x axis. Julia has implementations in the DSP package. Here we introduce the usage by looking at practical examples.

A Savitzky-Golay filter can be used to smooth data before applying any further analysis. The idea is to use local polynomial approximations to replace the original data.

Spectrogram

Spectrogram is used a lot in wave analysis. For my purpose, I use it as an approach to visualize time dependent simulation data along a continuous line region.

Minimum Variance Analysis

A nice introduction is given by Bengt U.Ö.Sonnerup and Maureen Scheible. Here is a brief summary of the idea. The implementation of MVA can be found in MVA.jl.

The main purpose of minimum or maximum variance analysis (MVA) is to find, from single-spacecraft data, an estimator for the direction normal to a one-dimensional or approximately one-dimensional current layer, wave front, or other transition layer in a plasma.

For real transition layers observed in space there are usually more or less pronounced deviations from the ideal 1-D model. The layer is likely to have 2-D or 3-D internal structures which evolve in time and to have temporal fluctuations in the orientation of its normal as well.

The minimum variance technique is designed to deal with the situation where some or all of the non-ideal effects mentioned above, except a systematic temporal change in the normal direction, $\widehat{n}$, are present. As the estimate of $\widehat{n}$, the method identifies that direction in space along which the field-component set {$\mathbf{B}^{(m)}\cdot\widehat{n}$} $(m = 1, 2, 3...M)$ has minimum variance. In other words, $\widehat{n}$ is determined by minimisation of

where the average $\langle\mathbf{B}\rangle$ is defined by

and where the minimisation is subject to the normalisation constraint $|\widehat{n}|=1$. Using a Lagrange multiplier $\lambda$ to implement this constraint, one then seeks the solution of the set of three homogeneous linear equations

where $\sigma^2$ is given by the equation above and $\widehat{n}$ is represented in terms of its three components $(n_x, n_y, n_z)$ along the cartesian coordinate system X, Y, Z (e.g., GSE or GSM) in which the field data $\{\mathbf{B}^{(m)} \}$ are given. When the differentiations in equations above have been performed, the resulting set of three equations can be written in matrix form as

where the subscripts $\mu,\nu = 1,2,3$ denote cartesian components along the X, Y, Z system and

is the magnetic variance matrix. It is seen from the equation that the allowed $\lambda$ values are the eigenvalues $\lambda_1,\lambda_2,\lambda_3$ (given here in order of decreasing magnitude) of $M_{\mu\nu}^{B}$. Since $M_{\mu\nu}^{B}$ is symmetric, the eigenvalues are all real and the corresponding eigenvectors, $x_1$, $x_2$, and $x_3$, are orthogonal. The three eigenvectors represent the directions of maximum, intermediate, and minimum variance of the field component along each vector.

A not so good example case is the Galileo G8 Ganymede flyby magnetometer data: where the ratio between the intermediate and minimum eigenvalue is only about 2.

Mean Field-Aligned Coordinates

The magnetic field satellite data are usually referred to geocentric coordinate reference frames such as geocentric solar ecliptic (GSE) frame. Conversely, the MHD waves modes in magnetized plasma depend on the ambient magnetic field, and is then useful to rotate the magnetic field measurements into the mean field-aligned (MFA) coordinate system. This reference frame is useful to study the ultra-low frequency magnetic field variations along and perpendicular to the direction of the mean field.

The MFA coordinate system is a local coordinate system that relates to the geomagnetic field. In the ecliptic plane, the parallel direction is more or less aligned with z in GSE, the first perpendicular direction pointing from Earth to the region of interest, and the second perpendicular direction completing the right hand rule.

In order to identify the mean magnetic field the classical moving average (MAVG) approach is usually adopted, but does not always give reliable performance. In the solarwind, a technique widely used to study parallel and perpendicular magnetic field components of waves is the minimum variance analysis (MVA). However, the minimum variance direction does not necessarily coincide with that of the ambient magnetic field.

If a time scale separation exists within the magnetic field measurements, these time series can be thought as a superposition of a slowly varying (amplitude) ambient field $\mathbf{B}_0(t)$, an higher frequency signal $\mathbf{b}(t)$(or perturbation) and an incoherent noise $\mathbf{n}(t)$:

MFA coordinate system in the magnetosphere: the MFA directions $\mathbf{\mu}$, $\mathbf{\phi}$ and $\mathbf{\nu}$ are showed together with the satellite position (red dot). The geomagnetic field line (blue line) is computed through T96 magnetosphere model [Tsyganenko 1995] during solar quiet conditions. [Regi+, 2016][Regi2016]

The MFA coordinates system, showed in the figure for an assigned position in the inner magnetosphere, is established by means of the unit vectors defined as

where $\mathbf{\mu}$, $\mathbf{\phi}$ and $\mathbf{\nu}$ are usually associated with compressional, toroidal and poloidal ULF waves modes respectively, while $\mathbf{r}(t)$ represents the position vector of the spacecraft.

This definition may also be extended in upstream regions, but in this case the $\mathbf{\phi}, \mathbf{\nu}$ components are simply related with transverse oscillations in the interplanetary region (e.g. foreshock upstream waves).

Then we can define the instantaneous rotation matrix from geocentric to MFA reference frame as

that allows us to project the instantaneous magnetic field vector from the original geocentric reference frame into the MFA one:

It is clear that in order to obtain the time series $\mathbf{B}^{MFA}(t)$ in the MFA reference system, the slowly varying mean field vector $\mathbf{B}_0(t)$ must be found through an appropriate filtering procedure. The moving average (MAVG) is a good choice with respect to any low-pass filter, as it introduces no artifacts. However in applying the MAVG method it is assumed that the characteristic fluctuation time $T_0$ related to $\mathbf{B}_0(t)$ is much greater than the period $T_b$ of the perturbation $\mathbf{b}(t)$. Moreover, $T_0$ depends on both satellite motion and natural magnetic field variation (e.g. high velocity stream, coronal mass ejections, corotating interaction regions), and could be related to nonlinear and non stationary phenomena. Under these conditions the MAVG might be unsuitable in the rotation procedure, while a method such as the empirical mode decomposition (EMD), is useful to identify nonlinear and nonstationary processes.

This is discussed thoroughly in [Regi+, 2016][Regi2016].

ULF Wave Detection

ULF waves are MHD waves: Alfvén wave, fast wave and slow wave. One basic approach to identify waves is to check the correlation of quantity perturbations.

The phase speed of shear Alfvén wave is

where $v_A$ is the Alfvén speed and $\theta$ is the angle between wave vector $\mathbf{k}$ and magnetic field $\mathbf{B}$.

The perturbed quantities of Alfvén waves follow these relations:

where $\delta \mathbf{v}$, $\delta \mathbf{B}$, and $\delta \rho$ are perturbed plasma velocity, magnetic fields, and plasma density, respectively, and $B_0$ is the background magnetic magnitude.

For slow and fast waves, the phase speeds are

The "+" is for fast waves and "−" for slow waves, and $v_S$ is the sound speed. The perturbed quantities for fast and slow waves are

Thus generally the Alfvén wave is identified by the correlations between velocity and magnetic field perturbations, and the fast and slow waves are identified by the negative (for slow waves) or positive (for fast waves) correlations between either density and magnetic field perturbation or thermal pressure and magnetic pressure perturbation.

The second equation above can also be expressed in terms of magnetic and thermal pressure pertubations:

See the lecture notes for more details.

For the magnetosonic waves, consider using $\delta \mathbf{E}$ and $\delta \mathbf{B}$ for identifying speed. The slopes of the curves $\delta E∕\delta B$ correspond to the wave propagation speed in the spacecraft frame.

Transverse and shear Alfvén wave refer to actually the same thing: the descriptions arise from $\mathbf{k}\cdot\mathbf{V}=0$ and $\mathbf{V}\cdot\mathbf{B}_0=0$.

The fast and slow magnetosonic waves are associated with non-zero perturbations in the plasma density and pressure, and also involve plasma motion parallel, as well as perpendicular, to the magnetic field. The latter observation suggests that the dispersion relations $\omega=kV_{\pm}$ are likely to undergo significant modification in collisionless plasmas. In order to better understand the nature of the fast and slow waves, let us consider the cold-plasma limit, which is obtained by letting the sound speed $V_S$ tend to zero. In this limit, the slow wave ceases to exist (in fact, its phase velocity tends to zero) whereas the dispersion relation for the fast wave reduces to

This can be identified as the dispersion relation for the compressional-Alfvén wave. Thus, we can identify the fast wave as the compressional-Alfvén wave modified by a non-zero plasma pressure.

In the limit $V_A\gg V_S$, which is appropriate to low-$\beta$ plasmas, the dispersion relation for the slow wave reduces to

This is actually the dispersion relation of a sound wave propagating along magnetic field-lines. Thus, in low-$\beta$ plasmas the slow wave is a sound wave modified by the presence of the magnetic field.

In reality, waves can be mixed together with mode conversions. Also, notice that the classical wave theory is based on spatially homogeneous plasma assumption, which is rarely the case in nature such as the magnetosphere.

A tricky part in practice is how to get the average through smoothing. Note that a real satellite moves both in time and space. Usually people do moving-box-average to get an average state within a short period.

A more careful analysis is called Walén test.

However, always keep in mind that the most reliable way of identifying waves is to calculate the dispersion relation.

Wave Polarization Analysis

FFT, MVA and singular value decomposition (SVD) techniques are the common tools for analyzing wave signals from EM field observations. Here we present a set of single satellite field analysis following Santolík+ 2003.

SVD of the B-field

Early satellites do not have E field detector onboard, so we start from utilizing B field only. Assuming a single plane wave at frequency f. The linearized Faraday's law gives

where $\omega = 2\pi f$ and $\mathbf{k}$ is the wave vector. Apparently the transverse wave follows

Let $\hat{S}_{ij} = <\hat{B}_i\hat{B}_j^\ast>$, where $\hat{B}(f) = \textrm{fft}(B(t))$ and * denotes the complex conjugate. We have

or

There are three complex equations. For applying real SVD procedure, we turn these into six real equations

where $\mathbf{k}$ is now expanded into $(k_1, k_2, k_3, k_1, k_2, k_3)^T$ and the magnetic field spectral matrix A, which is obtained by performing short-time Fourier transforms (STFT) of the waveforms observed, can be written as^{[complex_property]}

where N denotes the number of spectral matrices averaged over time(???), $B_x^n$, $B_y^n$ and $B_z^n$ denote the three orthogonal components of the n-th magnetic field Fourier spectra in the Cartesian coordinate system for the obtained waveform.

The spectral matrix A can be decomposed by SVD as follows:

where $U$, $S$, and $V^T$ denote a $6\times 3$ matrix, a $3\times 3$ of three nonnegative singular values, and a $3\times 3$ matrix with orthonormal rows, respectively.

Now we have 6 real equations but only 2 unknowns. If this modified set does not degenerate to less than three equations it is not possible to find a vector equation $\mathbf{k}$ which simultaneously solves all six equations.^{[other_methods]} Instead, a "solution" can be defined in the least squares sense, meaning that we search for a column vector $\mathbf{k}$ which gives the minimum modulus of a six-dimensional vector $\hat{A}\cdot\mathbf{k}$. This is essentially what SVD does for us!

The wave normal polar angle $\theta$, polarization ellipse $L_p$ and planarity of polarization $F$ can be calculated as follows:

where $s_1$, $s_2$, and $s_3$ denote the ascending series of singular values $S$, and $v_1$, $v_2$, and $v_3$ are the elements in the row of $V^T$ corresponding to the minimum singular value at the diagonal of W. $|L_p|=0$ and 1 indicate linear polarization and circular polarization, respectively.

In special coordinate systems we may be able to determine the sign of ellipticity, but mathmatically there is an ambiguity.^[MAG] The planarity of polarization $F$ indicates the validity of the single plane-wave assumption. The SVD technique provides reasonable results when the single plane wave approximation is applicable. In that case, $F$ is near unity; otherwise, $F$ approaches zero.

It should be noted that $\theta_k$ has an ambiguity in the sign of polarity because the correlation between electric and magnetic field is not calculated. To determine the polarity of k-vector, we calculate the Poynting vector angle using electric and magnetic field complex spectra. The Poynting vector angle $\theta_p$ can be derived as $\theta_p = \cos^{-1}(P_3/\sqrt{P_1^2 + P_2^2 + P_3^2})$. $\theta_p=0^o$ and $180^o$ indicate parallel and anti-parallel propagation along the field line, respectively. $P_1$, $P_2$, and $P_3$ denote the three components parallel to the Poynting vector $\mathbf{P}$ as follows:

where $\hat{E}_x$, $\hat{E}_y$, and $\hat{E}_z$ denote the three orthogonal components of the electric field Fourier spectra.

SVD of the EM fields

With the vector of magnetic field fluctuations the propagation direction cannot be fully determined –- two antiparallel directions cannot be distinguished. We have to use both the magnetic and electric fields to fully determine the wave vector from the Faraday's law $\mathbf{k}\times\mathbf{E} = \omega \mathbf{B}$. For the sake of simplicity we rewrite it to

where $\mathbf{n} = \mathbf{k}c/\omega$ is a dimensionless vector having the modulus of the wave refractive index and the direction of the wave vector, c being the speed of light.

Suppose the vectors $\mathbf{B}$ and $\mathbf{E}$ are represented by three magnetic components and three electric components in the same Cartesian coordinate system. We now proceed similarly as in the B-only case but instead of using the condition of perpendicularity between $\mathbf{B}$ and $\mathbf{k}$ we can directly use the equation above. Multiplying the three complex equations successively by the three cartesian components of $c\mathbf{B}^\ast$, and by the three cartesian components of $\mathbf{E}^\ast$ we have 18 mutually dependent complex equations

where $\epsilon_{ijk}$ is the Kronecker delta, indexes j,k follows the Einstein summation, and index i = 1...3, l = 1...6. $\zeta$ is a 6-D "electromagnetic" vector defined as

Terms $E_k\zeta_l^\ast$ and $cB_i\zeta_l^\ast$ can be respectively written as selected components $Q_{(k+1)l}$ and $Q_{il}$ of a 6x6 spectral matrix Q

The result can be rewritten as a set of 36 real equations for 3 components of a column vector $\mathbf{n}$,

where $A_E$ is a real $36\times 3$ matrix, and $\mathbf{b}$ is a 36-dimension real column vector. More explicitly,

The system contains 36 equations for only 3 unknowns. However, in the ideal case of analytic signals these equations are obtained from only six real equations from Faraday's law. One of these six equations could be, moreover, eliminated using an appropriate phase factor for the complex amplitudes, and two of them are already contained in the condition $\mathbf{E}\cdot\mathbf{B} = 0$ which internally bounds the complex electric and magnetic field amplitudes. For analytic signals, the system thus should degenerate into only 3 independent real equations for the three components of $\mathbf{n}$.

With real data, the multidimensional spectral analysis at a given frequency results in an estimate $\hat{Q}$ of the spectral matrix which can be used instead of $Q$,

The system then becomes

where $\hat{A}_E$ and $\hat{\mathbf{b}}$ are composed in the same way as equation $A_E$ and $\mathbf{n}$, but this time we use components of $\hat{Q}$ instead of $Q$.

The matrix $\hat{A}_E$ can thus be decomposed,

where $U_E$, $S_E$, and $V_E^T$ denote a $36\times 3$ orthonormal matrix, a $3\times 3$ of three nonnegative singular values, and a $3\times 3$ matrix with orthonormal rows, respectively.

The least squares solution is then directly given by

Similar to the pure B-field case, we have ellipticity

and planarity

where $s_1$, $s_2$, and $s_3$ denote the ascending series of singular values $S_E$. This time we know the direction of the wave vector, so there's no ambiguity about the left or right hand polarizations.

Correlation Test Between Two Variables

This part takes the reference from R.

Correlation test is used to evaluate the association between two or more variables.

Methods for correlation analyses

There are different methods to perform correlation analysis:

• Pearson correlation (r), which measures a linear dependence between two variables (x and y). It’s also known as a parametric correlation test because it depends to the distribution of the data. It can be used only when x and y are from normal distribution. The plot of y = f(x) is named the linear regression curve.
• Kendall tau and Spearman rho, which are rank-based correlation coefficients (non-parametric).

Correlation formula

In the formula below,

• x and y are two vectors of length n
• $\bar{x}$ and $\bar{y}$ corresponds to the means of $x$ and $y$, respectively.

Pearson correlation formula

The p-value (significance level) of the correlation can be determined :

1. by using the correlation coefficient table for the degrees of freedom : $df=n−2$, where $n$ is the number of observation in $x$ and $y$ variables.

2 or by calculating the $t$ value as follows:

where the corresponding p-value is determined using t table distribution for $df=n-2$. If the p-value is $< 5\%$, then the correlation between $x$ and $y$ is significant.

Time-dependent correlation

In real lift, the properties of the signals may change as a function of time: different correlation relation may appear at different time intervals. Thus we usually need to perform the correlation analysis for each subsequent interval to obtain a time-dependent correlation, similar to the concept of short time Fourier transform.

[Regi2016]: https://doi.org/10.4401/ag-7067