# Background

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.

• 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. Dec 29, 2021 at 17:32