I am attempting to numerically solve the ideal MHD equations in normal mode form for a Harris current sheet. The linearized perturbed MHD equations can be written in a normal mode form: $$ -\rho^{}_{m0} \omega_n^2 \boldsymbol{\xi}_n = \mathbf{F}(\boldsymbol{\xi}_n) \tag{1} $$

where $\mathbf{F}(\boldsymbol{\xi}_n)$ is the linear force operator given by: $$ \mathbf{F}(\boldsymbol{\xi}_n) = \frac{1}{\mu_0}(\nabla \times \mathbf{B}_0) \times [\nabla \times (\boldsymbol{\xi}_n \times \mathbf{B}_0)] + \frac{1}{\mu_0} \Big\{\nabla \times \big[\nabla \times (\boldsymbol{\xi}_n \times \mathbf{B}_0)\big] \Big\}\times \mathbf{B}_0+ \nabla \Big[\boldsymbol{\xi}_n \cdot \nabla P_0 + \gamma P_0 \nabla \cdot \boldsymbol{\xi}_n \Big] $$

Here, $\mathbf{B}_0$ is the equilibrium magnetic field given by: $$ \mathbf{B}_0 = B_0 \tanh\left(\frac{y}{a}\right) \hat{\mathbf{x}}+B_T \hat{\mathbf{z}} $$

where $B_0$ and $B_T$ are constants. $\boldsymbol{\xi}_n$ is the plasma displacement vector in two dimensions (the focus of our problem is in 2D), hence: $$ \boldsymbol{\xi}_n = \xi_x \hat{\mathbf{x}} + \xi_y \hat{\mathbf{y}} $$

and finally, $\gamma$ is the polytrope index and $P_0$ is the equilibrium pressure which can be shown to be, using the momentum equation: $$ P_0 = - \frac{B_0^2}{2\mu_0}\tanh^2\left(\frac{y}{a}\right) + C $$

where $C$ is a constant. Furthermore, we assume the plasma is incompressible, that is $\nabla \cdot \boldsymbol{\xi}_n = 0$.

Actual Problem

After a super long and tedious calculation, equation (1) takes on the simple form (I dropped the subscript $n$ for simplicity): $$ \frac{B_0^2}{\mu_0} \tanh^2 \left(\frac{y}{a}\right) (\partial_{xx} + \partial_{yy}) \xi_y + \frac{2B_0^2}{\mu_0a} \text{sech}^2\left(\frac{y}{a}\right)\tanh\left(\frac{y}{a}\right) \partial_y \xi_y = -\rho^{}_{m0} \omega_n^2 \xi_y \tag{2} $$

and $\xi_x = 0$. We further simplify the problem by assuming that $\xi_y = f(y) e^{ikx - \omega t}$ and hence equation (2) becomes: $$ \frac{B_0^2}{\mu_0} \tanh^2 \left(\frac{y}{a}\right) (-k^2 + \partial_{yy}) \xi_y + \frac{2B_0^2}{\mu_0a} \text{sech}^2\left(\frac{y}{a}\right)\tanh\left(\frac{y}{a}\right) \partial_y \xi_y = -\rho^{}_{m0} \omega_n^2 \xi_y \tag{3} $$

Equation (3) looks like a fairly simple eigenvalue problem (I want to solve for the eigenvalues $\omega^2$) that can be solved using the central difference approximation: $$ \frac{B_0^2}{\mu_0} \tanh\left(\frac{y_i}{a}\right) \frac{1}{h} \left[\frac{1}{h}\tanh\left(\frac{y_i}{a} \right)- \frac{1}{a}\text{sech}^2\left(\frac{y_i}{a}\right)\right]f_{i-1} + \frac{B_0^2}{\mu_0} \tanh^2\left(\frac{y_i}{a}\right)\left(-\frac{2}{h^2} - k^2 \right)f_{i} + \frac{B_0^2}{\mu_0} \tanh\left(\frac{y_i}{a}\right) \frac{1}{h} \left[\frac{1}{h}\tanh\left(\frac{y_i}{a} \right)+\frac{1}{a}\text{sech}^2\left(\frac{y_i}{a}\right)\right]f_{i+1} = -\rho^{}_{m0} \omega_n^2 f_i $$

Considering above to be a sparse matrix, I used both Python's scipy.sparse.linalg.eigs and scipy.linalg.eig to solve this eigenvalue equation for the first few smallest eigenvalues and normal modes but got nonsensical results. For example, when I set $B_0 = 1$, $\rho_{m0} = 1$, $k = 0.1$ and $a = 0.1$, I get complex eigenvalues $\omega_n^2$ when using the former and negative values for $\omega_n^2$ when I use the latter. The first is incorrect because it is well known in MHD that $\omega_n^2$ is either real or imaginary. The second is incorrect because I know the eigenmodes should be stable and hence should have $\omega_n^2 > 0$.

I cannot understand where I went wrong in my implementation and I am looking for suggestions if there is a better method to go about it? I am not sure if finite difference is good for this kind of problems.

  • $\begingroup$ It may also help if you include your code for the implementation. Your equations may be correct, but the transcription into Python may have errors. $\endgroup$
    – Tyberius
    Dec 29, 2021 at 17:32


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy