2 The Okubo–Weiss-type criterion for MHD flows

The equation governing the transport of the magnetic field ${\boldsymbol {B}} = \langle B_1, B_2 \rangle$ in 2-D (or 3-D) incompressible MHD flow is (Goedbloed & Poedts Reference Goedbloed and Poedts2004), in usual notation,

(2.1a)\begin{align} \frac{\partial {\boldsymbol{B}}}{\partial t} = \boldsymbol{\nabla} \times ( {\boldsymbol{v}} \times {\boldsymbol{B}} ) \end{align}

or

(2.1b)\begin{align} \frac{{\rm D} {\boldsymbol{B}}}{{\rm D} t} = \frac{\partial {\boldsymbol{B}}}{\partial t} + \left( {\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) {\boldsymbol{B}} = ({\boldsymbol{B}} \boldsymbol{\cdot} \boldsymbol{\nabla} ){\boldsymbol{v}} \equiv \boldsymbol{\mathcal{A}} \boldsymbol{\cdot}{\boldsymbol{B}}, \end{align}

where $\boldsymbol {\mathcal {A}}$ is the flow-velocity gradient matrix and ${\boldsymbol {v}}$ is defined by ${\boldsymbol {v}} = \langle u, v \rangle$.

It may be of interest to note that (2.1a) for the 2-D transport of a magnetic field in MHD is the same as that for the 2-D transport of the fluid divorticity (which is the curl of the fluid vorticity). This may be understood by noting the well-known analogy between vorticity in 3-D hydrodynamics and a magnetic field in 3-D MHD, and the analogy between vorticity in 3-D hydrodynamics and divorticity in 2-D hydrodynamics (Kuznetsov et al. Reference Kuznetsov, Naulin, Nielsen and Rasmussen2007).

If we assume, following Okubo (Reference Okubo1970) and Weiss (Reference Weiss1991), the straining flow-velocity gradient tensor $\boldsymbol {\nabla }{\boldsymbol {v}}$ to temporally evolve slowly, so the magnetic field can be taken to evolve adiabatically with respect to the straining flow-velocity gradient field,Footnote ¹ (2.1b) may be locally approximated by an eigenvalue problem, with eigenvalues given by

(2.2)\begin{align} \lambda^2 = u_y v_x + {v_y}^2 \equiv Q. \end{align}

An interesting interpretation of the Okubo–Weiss-type parameter $Q$ for the MHD case follows upon noting that the slow-flow-variation ansatz used above may be approximated by the MHD stationary Beltrami flow state (Shivamoggi Reference Shivamoggi2011). From (2.1a), the latter corresponds to the generalized Alfvénic state

(2.3)\begin{equation} {\boldsymbol{B}} = a {\boldsymbol{v}}, \end{equation}

with $a$ being an arbitrary constant.Footnote ² Equation (2.3) implies that the vorticity $\boldsymbol {\omega } = \boldsymbol {\nabla } \times {\boldsymbol {v}}$ is proportional to the electric current density ${\boldsymbol {J}} = \boldsymbol {\nabla } \times {\boldsymbol {B}}$ for a generalized Alfvénic state in 2-D MHD. Using (2.3), (2.2) becomes

(2.4)\begin{equation} \lambda^2\equiv Q = \frac{1}{a^2} \left(\displaystyle B_{1y} {B_2}_x + {B^2_{2y}} \right)=\frac{1}{4a^2}\left( b_1^2+ b_2^2- J^2/c^2\right), \end{equation}

where

(2.5a–c)\begin{equation} b_1\equiv-2B_{2y}, \quad b_2\equiv{B_2}_x+B_{1y},\quad\frac{J}{c}\equiv B_{2x}-B_{1y}.\end{equation}

Equation (2.4) shows that the Okubo–Weiss-type parameter $Q$ is a measure of the relative importance of magnetic shear ($Q>0$, hyperbolic) and electric current ($Q < 0$, elliptic).Footnote ³ This result, in conjunction with (2.2) and (2.3), is in accord with the eigenvalues of the magnetic field gradient tensor $\boldsymbol {\nabla }{\boldsymbol {B}}$ becoming purely imaginary near the $O$-points or real near the $X$-points of the magnetic field lines (Greene Reference Greene1993).

The hyperbolic magnetic field topology is of great interest in plasmas in space and fusion devices because it becomes a seat for magnetic reconnection processes (Giovanelli Reference Giovanelli1949; Shivamoggi Reference Shivamoggi1986, Reference Shivamoggi1999; Rollins & Shivamoggi Reference Rollins and Shivamoggi2007).

In terms of the magnetic vector potential ${\boldsymbol {A}}$ given by

(2.6a,b)\begin{equation} {\boldsymbol{B}} \equiv \boldsymbol{\nabla} \times {\boldsymbol{A}}, \quad {\boldsymbol{A}} = A \,{\hat{\boldsymbol{i}}}_z, \end{equation}

the Okubo–Weiss-type parameter (2.4) becomes

(2.7)\begin{equation} Q = \frac{1}{a^2} \left[ \left( \frac{\partial^2 A}{\partial x \partial y}\right)^2 - \frac{\partial^2 A}{\partial x^2} \frac{\partial^2 A}{\partial y^2} \right]\!. \end{equation}

Equation (2.7) shows that the Okubo–Weiss-type parameter $Q$ for the MHD case is the negative of the Gaussian curvature of the magnetic flux surface, to within a positive multiplicative factor. Since the 2-D MHD magnetic field topology is characterized by the local Gaussian curvature of the underlying magnetic flux manifold, (2.2) has the potential to serve as a useful diagnostic tool to parameterize the magnetic field topology in 2-D MHD flows in terms of magnetic-shear-dominated and electric-current-dominated regions. This is numerically investigated in § 5.

This development may be formally extended to the 2-D EMHD case via the generalized magnetic flux framework to incorporate the effects of electron inertia. The Okubo–Weiss-type parameter for the 2-D EMHD case turns out to be related to the intrinsic metric properties of the out-of-plane magnetic field manifold (see Appendix A).

Example 2.1 As an example of the above formulation, consider a MHD flow with the magnetic flux function given by

(2.8a)\begin{equation} A( x,y)=\frac{k}{2}\left( \alpha x^2-y^2\right)\!. \end{equation}

So, the magnetic field is given by the hyperbolic/elliptic configuration near an $X(O)$-type $(\alpha \gtrless 0)$ magnetic neutral point

(2.8b)\begin{equation} B_1=-ky,\quad B_2=-k\alpha x.\end{equation}

Thus, we obtain

(2.9a–c)\begin{equation} b_1=0,\quad b_2=-k( 1+\alpha),\quad \frac{J}{c}=k( 1-\alpha),\end{equation}

which shows that the current density $J\not =0$, unless $\alpha = 1$. On the other hand, the magnetic field topology is determined only by whether $\alpha \lessgtr 0$, as shown below.

Using (2.9a–c), (2.4) gives

(2.10)\begin{equation} Q=\left( \frac{k^2}{a^2}\right) \alpha \gtrless 0, \quad\alpha \gtrless 0, \end{equation}

as to be expected.

3 Reformulation in polar coordinates

In plane polar coordinates $(r,\theta )$, (2.1b) written in the component form is

(3.1a)\begin{gather} \frac{\partial B_r}{\partial t}+v_r\frac{\partial B_r}{\partial r}+v_\theta\frac{\partial B_r}{r\partial\theta}-\frac{v_\theta B_\theta}{r}=B_r\frac{\partial v_r}{\partial r}+B_\theta\frac{\partial v_r}{r\partial \theta}-\frac{B_\theta v_\theta}{r} , \end{gather}

(3.1b)\begin{gather}\frac{\partial B_\theta}{\partial t}+v_r\frac{\partial B_\theta}{\partial r}+v_\theta\frac{\partial B_\theta}{r\partial\theta}+\frac{v_\theta B_r}{r}=B_r\frac{\partial v_\theta}{\partial r}+B_\theta\frac{\partial v_\theta}{r\partial \theta}+\frac{B_\theta v_r}{r}. \end{gather}

The velocity gradient matrix $\boldsymbol {\mathscr {A}}$ in (2.1b) then becomes

(3.2)\begin{equation} \boldsymbol{\mathscr{A}}\equiv\left[ \begin{array}{@{}cc@{}} \dfrac{\partial v_r}{\partial r} & \dfrac{\partial v_r}{r\partial\theta}-\dfrac{v_\theta}{r}\\[10pt] \dfrac{\partial v_\theta}{\partial r} & \dfrac{\partial v_\theta}{r\partial\theta}+\dfrac{v_r}{r} \end{array} \right]\!. \end{equation}

On using the mass-conservation equation

(3.3)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r}\left( rv_r\right)+\frac{\partial v_\theta}{r\partial\theta}=0, \end{equation}

and assuming again that the magnetic field evolves adiabatically with respect to the straining flow-velocity gradient field, (3.1a,b) may again be locally approximated by an eigenvalue problem, with eigenvalues given by

(3.4)\begin{equation} \lambda^2=\left(\frac{\partial v_r}{\partial r}\right)^2+ \displaystyle \frac{1}{r}\frac{\partial v_r}{\partial\theta}\frac{\partial v_\theta}{\partial r}-\frac{1}{r}v_\theta\frac{\partial v_\theta}{\partial r}\equiv Q. \end{equation}

Equation (3.4) is the same as that for the divorticity field in 2-D hydrodynamics (Shivamoggi et al. Reference Shivamoggi, Van Heijst and Kamp2022).

On using the stationary generalized Alfvénic state condition (2.3), (3.4) becomes

(3.5a)\begin{equation} \lambda^2=\frac{1}{a^2}\left[\left(\frac{\partial B_r}{\partial r}\right)^2+\frac{1}{r}\frac{\partial B_r}{\partial\theta}\frac{\partial B_\theta}{\partial r}-\frac{1}{r}B_\theta\frac{\partial B_\theta}{\partial r}\right]\equiv Q.\end{equation}

On using the Gauss law

(3.6)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r}\left( r B_r\right) +\frac{\partial B_\theta}{r\partial\theta}=0, \end{equation}

equation (3.5a) may be alternatively expressed as

(3.5b)\begin{equation} \lambda^2=\frac{1}{a^2}\left[\left(\frac{\partial B_\theta}{\partial r}\right)\left(\frac{\partial B_r}{r\partial\theta}-\frac{B_\theta}{r}\right)+\left(\frac{B_r}{r}+\frac{\partial B_\theta}{r\partial\theta}\right)^2\right]\equiv Q.\end{equation}

In terms of the magnetic vector potential ${\boldsymbol {A}} = A \hat {{\boldsymbol {i}}}_z$, given by (2.6a–c), (3.5b) becomes

(3.7)\begin{equation} \lambda^2=\frac{1}{a^2r^4}\left[-r^2\frac{\partial^2 A}{\partial r^2}\left(\frac{\partial^2A}{\partial\theta^2}+r\frac{\partial A}{\partial r}\right)+\left(\frac{\partial A}{\partial\theta}-r\frac{\partial^2 A}{\partial r\partial\theta}\right)^2\right] \equiv Q. \end{equation}

Equation (3.7) shows that $\lambda ^2$ is the negative of the Gaussian curvature of the magnetic flux manifold, to within a positive multiplicative factor.

Example 3.1 As an example, consider an axisymmetric MHD flow with magnetic flux function given by

(3.8a)\begin{equation} A=A( r). \end{equation}

Equation (3.8a) leads to the magnetic field

(3.8b)\begin{equation} B_r=0,\quad B_\theta=-\frac{{\rm d}A}{{\rm d}r}. \end{equation}

Using (3.8), (3.7) gives for the eigenvalues

(3.9)\begin{equation} \lambda^2=-\frac{1}{r}\frac{{\rm d}A}{{\rm d}r}\frac{{\rm d}^2A}{{\rm d}r^2}. \end{equation}

Consider the $Z$-pinch problem, where the magnetic field is generated by a uniform, unidirectional current ${\boldsymbol {J}}=J{\hat {\boldsymbol {i}}}_z$,Footnote ⁴

(3.10)\begin{equation} A=-\frac{1}{4c}Jr^2 . \end{equation}

Using (3.10), (3.9) becomes

(3.11)\begin{equation} \lambda^2=-\frac{J^2}{4c^2}<0, \end{equation}

as to be expected.

Consider next a magnetic field generated by a current-carrying cylinder $(0< r< r_0)$

(3.12)\begin{equation} J=J_0h(r_0-r), \end{equation}

where $h\left ( x\right )$ is the unit step function, and $J_0$ is the uniform current density in this cylinder. The magnetic flux function associated with (3.12) is given by

(3.13)\begin{equation} A=\begin{cases} \displaystyle-\frac{1}{2}\frac{J_0r^2}{c}, & r< r_0\\[10pt] \displaystyle-\frac{J_0r_0^2}{c}\ln(r/r_0), & r>r_0, \end{cases} \end{equation}

and the magnetic field is given by

(3.14)\begin{equation} B_\theta= \begin{cases} \displaystyle\frac{J_0r}{c}, & r< r_0\\[10pt] \displaystyle\frac{J_0r_0^2}{rc}, & r>r_0, \end{cases} \end{equation}

which represents an MHD Rankine-type flux tube.

Using (3.14), (3.5a) gives

(3.15)\begin{equation} \lambda^2= \begin{cases} \displaystyle-\frac{J_0^2}{a^2c^2}, & r< r_0\\[10pt] \displaystyle\frac{J_0^2r_0^4}{c^2a^2r^4}, & r>r_0, \end{cases} \end{equation}

signifying an elliptic region inside the current-carrying cylinder and a hyperbolic region outside it, as to be expected.

4 Numerical set-up

The Okubo–Weiss-type criteria for 2-D MHD flows are now illustrated by performing numerical simulations of the relevant equations in dimensionless form. The fluid flow occurs in a planar 2-D domain, with the velocity field given by ${\boldsymbol {v}}=u(x,y,t){\hat {\boldsymbol {i}}}_x+v(x,y,t){\hat {\boldsymbol {i}}}_y$, where ${\hat {\boldsymbol {i}}}_x$ and ${\hat {\boldsymbol {i}}}_y$ are the unit vectors in the $x$- and $y$-directions, respectively. Additionally, there is a (dimensionless) magnetic field, which is normalized using some characteristic magnetic field strength $B_0$, ${\boldsymbol {B}}=B_1(x,y,t){\hat {\boldsymbol {i}}}_x+B_2(x,y,t){\hat {\boldsymbol {i}}}_y$.

For MHD flows, the governing (dimensionless) equations are the equation of motion

(4.1)\begin{equation} \frac{\partial{\boldsymbol{v}}}{\partial t}+\left({\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right){\boldsymbol{v}}=- \frac{1}{\rho} \boldsymbol{\nabla} p+ \frac{1}{R_e}\nabla^2{\boldsymbol{v}}+J{\hat{\boldsymbol{i}}}_z\times{\boldsymbol{B}}, \end{equation}

where ${\hat {\boldsymbol {i}}}_z$ is the unit vector in the $z$-direction and the magnetic field transport equation

(4.2)\begin{equation} \frac{\partial{\boldsymbol{B}}}{\partial t}=\boldsymbol{\nabla}\times\left({\boldsymbol{v}}\times {\boldsymbol{B}}\right)+\frac{1}{R_m}\nabla^2{\boldsymbol{B}}. \end{equation}

In (4.1), velocity has been non-dimensionalized using a reference speed $U$, and length has been normalized using a reference length $L$; $R_e$ is the Reynolds number given by $R_e\equiv UL/\nu$, where $\nu$ is the kinematic viscosity of the fluid; $R_m$ is the magnetic Reynolds number given by $R_m\equiv \mu _0UL/\eta$, $\mu _0$ being the free space permeability; and $\eta$ is the resistivity. The electric current density $J$ follows from Ampère's law $J=(\partial B_2/\partial x-\partial B_1/\partial y)$.

Equations (4.1) and (4.2) have been numerically solved with a finite-element code Comsol,Footnote ⁵ using double-periodic boundary conditions on a computational domain defined by $-5\leq x\leq 5$ and $-5\leq y\leq 5$. The velocity field has been initialized with a statistically steady turbulent flow is field given by $u(x,y,t=0)=u_{0}(x,y)$ and $v(x,y,t=0)=v_{0}(x,y)$, with root-mean-square value equal to 1. In consistency with the generalized Alfvénic state condition (2.3), the magnetic field has also been initialized with this turbulent field, i.e. $B_1(x,y,t=0)=u_0(x,y)$ and $B_2(x,y,t=0)=v_0(x,y)$. This corresponds to putting $a=1$ (on taking $\rho =1$, plasma being incompressible) in (2.3), which as mentioned in footnote 2, reflects the nonlinear Alfvénic state associated with perturbations in a uniformly magnetized plasma (Hasegawa Reference Hasegawa1985). All numerical simulations were performed with $R_e=1000=R_m$. Integrations have been carried out using a second-order backward-difference scheme and an implicit time stepping.

This statistically steady turbulent state is set up by starting with a quiescent fluid which is mechanically forced using an external body force acting on the fluid. The details of this forcing become irrelevant once a turbulent flow develops, and then the forcing is turned off. The kinetic energy of the flow increases with time, eventually reaching a quasi-steady value, which indicates the attainment of a statistically steady state. The root-mean-square value of the flow velocity then becomes quasi-steady, even though the flow is turbulent and time dependent. As time progresses, this turbulent flow state, which is not forced anymore, starts to decay concomitantly creating large-scale coherent structures. Figure 1 gives a transient snapshot of this time-dependent decaying turbulent flow showing these large-scale structures, which may be identified with an inverse cascade.

Figure 1. Snapshot of the vorticity field (a) and the associated Okubo–Weiss function (b) taken during the decaying phase of MHD turbulence. The colour bar denotes the value of the Okubo–Weiss function, which has been rescaled to range from $-1$ to $+1$.

Our numerical simulations are carried out using a finite-element method for which the computational domain is discretized using an ultra-fine triangular mesh – a typical mesh element has a size of 0.01 (in normalized units), resulting in $10^6$ mesh elements. The intended accuracy of the method is, among other things, determined by the number of mesh elements in the computational domain rather than the number of grid points (which are dealt with better by a finite-difference method rather than a finite-element method).

5 Numerical results

For a generalized Alfénic state in 2-D MHD, as discussed in § 2, the spatial pattern of the electrical current density ${\boldsymbol {J}}$ is similar to that of the vorticity field $\boldsymbol {\omega }$. Therefore, the regions of concentrated electrical current density can be visualized via the Okubo–Weiss function $Q = u_y v_x+v_y^2$. Figure 1 shows a snapshot of the numerically calculated vorticity field (a) and the associated Okubo–Weiss function (b) taken during the decaying phase of a MHD turbulent field that was initialized with the generalized Alfvénic state (2.3) and evolved according to (4.1) and (4.2). Note that $Q$ has been rescaled each time so as to have it range from $-1$ to 1.

For the MHD case, the pertinent quantity for consideration, as per (2.4), would be the electric current density ${\boldsymbol {J}}$. However, thanks to the initialization with the generalized Alfvénic state (2.3), the snapshot of the vorticity would be similar to that of the electric current density (as confirmed also by the direct numerical simulation (DNS) of decaying 2-D MHD by Kinney, McWilliams & Tajima Reference Kinney, McWilliams and Tajima1995). Figure 1 shows that coherent vortices are indeed located in elliptic regions ($Q<0$) while divorticity sheets are located in hyperbolic regions ($Q>0$). The Okubo–Weiss-type parameter given in § 2 is therefore found to track the calculated location of the magnetic-shear-dominated and electric-current-dominated regions for the MHD case very well.

It may be mentioned that similar numerical results were obtained by Banerjee & Pandit (Reference Banerjee and Pandit2014, Reference Banerjee and Pandit2019), verifying (2.4) given in our earlier preliminary arXiv preprint (Shivamoggi, van Heijst & Kamp Reference Shivamoggi, van Heijst and Kamp2016).