1. Introductory remarks and description of the model

The magnetohydrodynamic rotating shallow water model (MRSW) was introduced by Gilman (Reference Gilman2000) as a simple model of the solar tachocline. As is the case of usual shallow water models, it can be obtained by applying vertical averaging and the mean field hypothesis to the full primitive equations of magnetohydrodynamics in the (magneto-)hydrostatic approximation (Zeitlin Reference Zeitlin2013). Again, as in the case of usual shallow water models, taking an asymptotic limit of small Rossby and magnetic Rossby numbers, i.e. of strong rotation, results in a self-consistent system of (magneto-)quasi-geostrophic equations for slow motions, where fast waves are filtered out (Zeitlin Reference Zeitlin2013). The magnetohydrodynamic quasi-geostrophic (MQG) model is of use in studies of (quasi-) two-dimensional magnetohydrodynamic turbulence, (e.g. Tobias, Diamond & Hughes Reference Tobias, Diamond and Hughes2007) still in the context of the solar tachocline. In the present paper we construct exact solutions of MQG equations, the magnetic modons, which are steady-moving dipoles of both magnetic field and vorticity. With the help of direct numerical simulations with the Dedalus code (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), which we adapted to the MQG simulations, we show that these solutions can be stable or unstable, depending on the intensity of the magnetic field inside and outside the dipole. The construction of the modon solutions follows the lines of that used recently (Lahaye, Zeitlin & Dubos Reference Lahaye, Zeitlin and Dubos2020) to obtain similar solutions in the quasi-geostrophic limit of the rotating shallow water model with horizontal density and/or temperature gradients, the so-called thermal rotating shallow water (TRSW, which becomes TQG in the quasi-geostrophic limit). It is a generalization of the classical procedure of building the modon solutions of the ‘ordinary’ quasi-geostrophic (QG) equations, first obtained by Larichev & Reznik (Reference Larichev and Reznik1976) following the idea of construction of the Lamb's dipole in two-dimensional (2-D) hydrodynamics. Existence of such solutions emphasizes a known (Dellar Reference Dellar2003) structural similarity between MRSW and TRSW.

Let us recall the equations of the MRSW model in the rotating $(x,\,y)$ plane, i.e.

(1.1)\begin{gather} \partial_{t}\boldsymbol{v} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v} + f \hat{\boldsymbol{z}} \wedge \boldsymbol{v}+ g \boldsymbol{\nabla} h = \frac{1}{h} \boldsymbol{\nabla} \left( h \, \boldsymbol{b} \otimes \boldsymbol{b} \right), \end{gather}

(1.2)\begin{gather}\partial_{t}h + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( h \boldsymbol{v} \right) = 0, \end{gather}

(1.3)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \left( h \boldsymbol{b} \right) = 0, \end{gather}

(1.4)\begin{gather}\partial_{t}\boldsymbol{b} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{b} = \frac{1}{h} \boldsymbol{\nabla} \left( h \, \boldsymbol{v} \otimes \boldsymbol{b} \right), \end{gather}

where $f$ is the Coriolis parameter and $g$ is the gravity acceleration. The dynamical variables of the model are the thickness of the layer $h$, the velocity $\boldsymbol {v} = u\hat {\boldsymbol {x}} + v \hat {\boldsymbol {y}}$ and the magnetic field $\boldsymbol {b}= a\hat {\boldsymbol {x}} + b \hat {\boldsymbol {y}}$ in the plane, where $(\hat {\boldsymbol {x}}, \hat {\boldsymbol {y}}, \hat {\boldsymbol {z}})$ are unit vectors along the respective axes, $\boldsymbol {\nabla } = \hat {\boldsymbol {x}} \partial _{x}+ \hat {\boldsymbol {y}} \partial _{y}$. Note that the magnetic field is rescaled with magnetic permeability of the vacuum and density of the fluid in order to have the dimension of velocity (Alfvèn velocity corresponding to a given value of the magnetic field). The notation $\boldsymbol {\nabla } \boldsymbol {\cdot }\left ( \boldsymbol {A} \otimes \boldsymbol {B} \right )$ is a shorthand for tensor notation: the $i$th component of such an expression is given by $\partial _{j} A_i B_j$, with summation over repeated indices from 1 to 3. The Coriolis parameter $f$ is constant in the $f$-plane configuration, and is a linear function of $y$ in the beta-plane configuration. In the following we work with the former, extrapolation of our results to the latter being rather straightforward. The constraint (1.3) can be resolved by introducing the (scalar) magnetic potential $A$,

(1.5)\begin{equation} h \boldsymbol{b} = \hat{\boldsymbol{z}} \wedge \boldsymbol{\nabla} A. \end{equation}

Up to a sign $A$ is the vertical component of the standard magnetic potential and is, in fact, a magnetic streamfunction for the horizontal magnetic field integrated over the vertical extent of the shallow water layer. In continuity with the previous work (Zeitlin Reference Zeitlin2013) we keep calling it magnetic potential. If $L$ is a characteristic scale, $U$ is a characteristic velocity and $B$ is a characteristic value of the magnetic field, the ‘ordinary’ and magnetic Rossby numbers, $Ro = {U}/{f L}$ and $Ro_{m} = {B}/{f L}$, respectively, measure the relative strength of rotation. For the solar tachocline, $Ro \sim 0.1$ is a reasonable assumption, and observations of magnetic Rossby waves indicate that investigating the weak magnetic Rossby number limit is relevant (Tobias et al. Reference Tobias, Diamond and Hughes2007; Dikpati et al. Reference Dikpati, Belucz, Gilman and McIntosh2018; Zaqarashvili et al. Reference Zaqarashvili2021). If both Rossby numbers are small, and of the same order, as well as the deviations of the free surface from its mean (rest) value $H$, and the characteristic scale $L$ is of the order of the deformation radius $R_{d} = {\sqrt {g H}}/{f}$ (i.e. the Burger number $Bu = R^{2}_{d}/L^{2} = {O}(1)$), the velocity in the leading order is divergenceless, and can be expressed in terms of a streamfunction $\psi$, which is the leading-order thickness anomaly $h - H$. In the leading order of the asymptotic expansion in $Ro$, $Ro_{m}$, the MQG equations on the $f$-plane for $\psi$ and $A$ thus result (e.g. Zeitlin Reference Zeitlin2013),

(1.6a)\begin{gather} \partial_{t}\left( \nabla^2 \psi - \frac{1}{Bu}\psi\right)+ \mathcal{J} (\psi, \nabla^2 \psi)- \mathcal{J} (A, \nabla^2 A) = 0, \end{gather}

(1.6b)\begin{gather}\partial_{t} A + \mathcal{J} (\psi, A) =0. \end{gather}

Here $\mathcal {J}(a,b)=\partial _x a\, \partial _y b - \partial _y a\, \partial _x b$ denotes the Jacobian, and the equations are written in non-dimensional form. Note that non-dimensional potential vorticity (PV) $q$ in this approximation is given by the formula

(1.7)\begin{equation} q = \nabla^2 \psi - \frac{1}{Bu}\psi. \end{equation}

Below we impose $Bu = 1$ for compactness, i.e. we use $R_{d}$ as spatial scale. Potential vorticity is not conserved in the presence of the magnetic field, and the equation (1.6a) describes its evolution. The energy conservation for the system (1.6) can be written in the standard form

(1.8)\begin{equation} E = \tfrac{1}{2} \int {{\rm d}\kern0.7pt x} \,{{\rm d} y} [(\boldsymbol{\nabla} \psi)^{2} + \psi^{2} + (\boldsymbol{\nabla} A)^{2}] = \mbox{const}. \end{equation}

Note, however, that in this form the energy is conserved only for $\psi$ and $A$ decaying at infinity, which is not the case for $A$ if the system evolves in an external magnetic field $\boldsymbol {B_{0}}$. In the latter case $A \rightarrow A + \hat {\boldsymbol {z}} \boldsymbol {\cdot } (\boldsymbol {x} \wedge \boldsymbol {B_0})$ in the formula (1.8).

Let us say again that the MQG system (1.6) is a QG limit of the system (1.1)–(1.4) at $(Ro, Ro_{m}) \rightarrow 0$. As a consequence, it does not contain fast magneto-inertia-gravity waves, although it retains long Alfvèn waves, cf. Zeitlin, Lusso & Bouchut (Reference Zeitlin, Lusso and Bouchut2015) and below. The MQG is mathematically equivalent to 2-D incompressible magnetohydrodynamics (2DMHD) with replacement of the vorticity, $\nabla ^2\psi +f$, in the latter by the PV $q$ in the former. Like the standard hydrodynamic QG model, which in the limit of infinite deformation radius, i.e. ${1}/{Bu} \rightarrow 0$, becomes equivalent to 2-D Euler equations for an incompressible fluid, the MQG model, in the same limit, becomes equivalent to 2DMHD.

Below we give an analytical derivation of the modon solutions in § 2, and then proceed in § 3 with numerical simulations initialized with these solutions, which we take either as are (§ 3.1) or perturbed by a weak noise (§ 3.2), and at different values of parameters. We then discuss the results in § 4.

2. Derivation of the magnetic modon solutions

We are looking for a localized solution that is stationary in a frame moving with constant velocity $U\hat {\boldsymbol {x}}$ (co-moving frame), so $\psi = \psi (x-Ut,y)$, $A = A(x-Ut,y)$. The MQG equations (1.6) thus become:

(2.1a)\begin{gather} -U\partial_x q + \mathcal{J}(\psi,q) - \mathcal{J}(A,\nabla^2 A) = 0, \end{gather}

(2.1b)\begin{gather}-U\partial_x A + \mathcal{J}(\psi,A) = 0. \end{gather}

Equation (2.1b) implies $\mathcal {J}(\psi +Uy,A)=0$; hence, $A=F(\psi +Uy)$, where $F$ is an arbitrary function. Following Lahaye et al. (Reference Lahaye, Zeitlin and Dubos2020) we will take it to be linear,

(2.2)\begin{equation} A = \kappa (\psi+Uy), \end{equation}

where $\kappa$ is a constant measuring the strength of magnetic potential relative to the streamfunction $\varPsi$ in the co-moving frame. Likewise, (2.1a) can be rewritten as

(2.3)\begin{equation} \mathcal{J}(\psi+Uy,q) - \mathcal{J}(A,\nabla^2 A) = 0. \end{equation}

Using the expression of $A$ in terms of $\psi$ (2.2) we obtain

(2.4)\begin{equation} \mathcal{J}(A,\nabla^2 A) = \kappa^2\mathcal{J}(\psi+Uy,\nabla^2\psi), \end{equation}

and, therefore, the PV evolution equation gives

(2.5)\begin{equation} \mathcal{J}\left(\psi+Uy, (1-\kappa^2)\nabla^2\psi - \psi\right) = 0. \end{equation}

This means that

(2.6)\begin{equation} (1-\kappa^2)\nabla^2\psi - \psi = G(\psi+Uy), \end{equation}

where $G$ is another arbitrary function, which we will again choose to be linear: $G(x) =\alpha x$, where $\alpha$ is another constant, which measures, like in the classical modon solution of Larichev & Reznik (Reference Larichev and Reznik1976), the strength of the PV anomaly in terms of the streamfunction in the co-moving frame.

Following the standard derivation (Larichev & Reznik Reference Larichev and Reznik1976), we divide the whole plane into inner ($-$) and outer ($+$) regions, with a circular separatrix at some radius $r = \sqrt {x^{2} + y^{2}}=r_0$. The parameters $\alpha$ and $\kappa$ are allowed to take different values in the outer and inner domains: $\alpha _\pm$ and $\kappa _\pm$. We thus have the following equations for the streamfunctions $\psi _{\pm }$ in the outer ($+$) and inner ($-$) domains, respectively:

(2.7)\begin{equation} (1-\kappa_{{\pm}}^2)\nabla^2\psi_{{\pm}} - \psi_{{\pm}} = \alpha_{{\pm}}(\psi_{{\pm}}+Uy). \end{equation}

Their solutions should be matched using the standard boundary conditions of continuity of velocity and pressure at the separatrix.

Solution in the outer region $r\geqslant r_0$: we are looking for a localized, finite-energy, where the energy is defined in (1.8), solution which should vanish far away from the centre. Hence, $\alpha _+ = 0$ and the following linear equation results:

(2.8)\begin{equation} (1-\kappa_+^2)\nabla^2\psi - \psi = 0. \end{equation}

A decaying at $r \rightarrow \infty$ solution is sought by separation of variables in polar coordinates $(r, \theta )$. Keeping in mind that matching with the inner solution is performed in the co-moving frame, and that a constant zonal velocity corresponds to a linear in $y = r \sin \theta$ streamfunction, we get

(2.9)\begin{equation} \psi_+= C_+ {\rm K}_1\left( r/\sqrt{1-\kappa_+^2}\right) \sin \theta. \end{equation}

Here ${\rm K}_{1}$ is the modified Bessel function of the second kind. Note that the decaying solution is possible only if $1-\kappa _+^2 > 0$, and that the particular case $\kappa _+ =0$ corresponds to the magnetic field being confined in the inner region, while if $\kappa _+ \neq 0$ there is a constant background zonal magnetic field $B = -\partial _{y} A = -\kappa _{+} U$ everywhere in the plane, cf. (2.2).

Solution in the inner region, $r\leqslant r_0$: the linear equation to solve is

(2.10)\begin{equation} (1-\kappa_-^2)\nabla^2\psi - \psi = \alpha_-(\psi+Uy). \end{equation}

The solution is represented as a superposition of a solution to the homogeneous problem, $\psi ^h$ and a particular solution $\psi ^p$ of the inhomogeneous problem,

(2.11a)\begin{gather} (1-\kappa_-^2)\nabla^2\psi^{h} - (1+\alpha_-)\psi^{h} = 0, \end{gather}

(2.11b)\begin{gather}\psi^p ={-}\frac{\alpha_- Uy}{1+\alpha_-}. \end{gather}

Anticipating matching with the outer solution, we impose the same angular structure and, hence, obtain

(2.12)\begin{equation} \psi_-= \left[C_- {\rm J}_1(\lambda r) - \frac{\alpha_- Ur}{1+\alpha_-}\right]\sin\theta, \end{equation}

where ${\rm J}_{1}$ is the Bessel function of the first kind, and $\lambda ^2 = - ({1 + \alpha _-})/({1-\kappa _-^2})$. Note a singularity at $\kappa _-=1$.

Matching: following Larichev & Reznik (Reference Larichev and Reznik1976) we impose the condition of continuity of the co-moving streamfunction, which provides the maximal smoothness of the resulting solution,

(2.13)\begin{equation} \left.\left( \psi_{+} + U y \right) \right|_{r = r_0} = \left. \left( \psi_{-} + U y \right) \right|_{r = r_0} = 0, \end{equation}

which allows us to determine $C_{\pm }$,

(2.14a,b)\begin{equation} C_{+} ={-} \frac{U r_0}{{\rm K}_{1} \left(r_0/\sqrt{1 - \kappa_{+}}\right)}, \quad C_{-} ={-} \frac{U r_0}{(1 + \alpha_{-}) {\rm J}_{1} (\lambda r_0)}. \end{equation}

The magnitude of the modon (streamfunction) is thus primarily driven by $Ur_0$ and modulated by $\alpha _-$. However, the latter is not a tunable parameter, but is such that the matching conditions are satisfied, and drives the radial structure of the streamfunction (through $\lambda$). Continuity of the radial derivative of the streamfunction $\partial _{r} \psi _{+} = \partial _{r} \psi _{-}$ at $r = r_0$ leads to the following transcendental equation for eigenvalues $\lambda$ at given $r_0$ and $\kappa _{+}$, and, hence, $\alpha _{-}$ at a given $\kappa _{-}$:

(2.15)\begin{equation} C_{+}/\sqrt{1 - \kappa^{2}_{+}}\,{\rm K}'_{1} \left(r_0/\sqrt{1 - \kappa_{+}}\right) = C_{-} \lambda {\rm J}'_{1} (\lambda r_0) - \alpha_{-} U/(1+\alpha_-), \end{equation}

where $C_{\pm }$ are given in (2.14a,b), and prime denotes differentiation of a function with respect to its argument. This transcendental equation is solved numerically using standard routines, e.g. with the fsolve function of the SciPy Python library.

The lowest eigenvalue corresponds to a steadily translating vortex dipole, a magnetic modon, the higher eigenvalues correspond to multiple changes of the sign of vorticity inside the separatrix, which are, presumably, less stable and more sensible do dissipation, with more sheared velocity, and will be discarded as is usually done. We fix $r_0 =1$ (i.e. the typical size of the modon is equal to the deformation radius, $Bu=1$) and, thus, get the eigenvalue $\alpha _{-}$ as a function of $\kappa _{-}$ and $\kappa _{+}$. As an example, we present in figure 1 the dependence $\alpha _{-}(\kappa _{-})$ when $\kappa _{+}= 0$. A gap in the curve is visible at $\kappa _-=1$. Indeed, for this value of $\kappa _-$, the interior equation (2.7) gives $\varPsi _-=-\alpha _-(\varPsi _-+Uy)$ which is incompatible with the boundary condition $\varPsi _-+Uy=0$ at $r=r_0$, meaning that there is no solution. We should emphasize that qualitatively different modon configurations result from different choices of the values of parameters $\kappa _{\pm }$: magnetic field confined inside the vortex dipole at $\kappa _{+}= 0$, $\kappa _{-} \neq 0$, hollow bubble, in the sense of magnetic field, propagating through the constant magnetic field at $\kappa _{-}= 0$, $\kappa _{+} \neq 0$, and a magnetized dipole propagating in the adverse, or collinear magnetic field, with respect to the magnetic field at the axis of the dipole, depending on the relative sign of $\kappa _{\pm } \neq 0$. Some of these cases are illustrated in figures 2–3.

Figure 1. Dependence of the first eigenvalue $\alpha _{-}$ on $\kappa _{-}$ at $\kappa _{+}= 0$. Curves for other values of $\kappa _+$ practically coincide with this one, as variations of $\alpha _-$ are dominated by variations of $\kappa _-$ at a given $r_0$. Disruption of the curve reflects the singularity at $\kappa _-=1$.

Figure 2. Magnetohydrodynamic quasi-geostrophic modon without external magnetic field. From left to right: streamfunction, vorticity and magnetic potential. Length in units of $R_{d}$ in this and subsequent figures.

Figure 3. Magnetic potential associated with a MQG modon without internal magnetic field ($\kappa _-=0, \kappa _+=0.02$, (a)), and a regular MQG modon with no tangential discontinuity of the magnetic field at the separatrix ($\kappa _-=\kappa _+=0.02$, (b)).

Before proceeding with numerical simulations, we should make two important remarks. First, we should stress that as the streamfunction $\psi + U y$ in the co-moving frame is continuous across the separatrix, and takes zero value, cf. (2.13), the magnetic potential is also continuous, but its derivative is not if $\kappa _{+} \neq \kappa _{-}$, even though the derivative of $\psi$ is continuous, cf. (2.2). This means that in this case we have a tangential discontinuity of the magnetic field across the separatrix. Such discontinuities are admissible in non-dissipative MHD (e.g. Landau & Lifshitz Reference Landau and Lifshitz1984), which is the parent of (1.6), but will be smeared in the presence of magnetic viscosity (resistivity). We will come back to this point in the discussion of numerical experiments below. As a consequence, we will call the modons with $\kappa _{+} = \kappa _{-}$ regular, and singular otherwise.

Second, as previously mentioned, in the case $\kappa _{+} \neq 0$ we have a configuration with a mean magnetic field which admits (rotational) Alfvèn waves. Unlike magneto-inertia-gravity waves, Alfvèn waves do not have a spectral gap at low frequencies in the presence of rotation (cf. Zeitlin et al. Reference Zeitlin, Lusso and Bouchut2015) and, thus, slow waves of this kind are not filtered out in the QG approximation which corresponds to the fast-time (${\sim }1/f$) averaging. Indeed, by linearizing the system (1.6) about a constant zonal magnetic field $B_{0}$, that is, by taking

(2.16)\begin{equation} A(x, y, t) = B_{0} y + A' (x, y, t), \end{equation}

where $A'$ is a small perturbation, in addition to considering a small streamfunction $\psi$, and substituting in the linearized equation (1.6a) the expression of $A'$ obtained from the linearized equation (1.6b) reduces the MQG equations to a wave equation for $\psi$,

(2.17)\begin{equation} \partial^{2}_{tt}\left( \nabla^2 \psi - \psi \right) - B^{2}_{0} \partial^{2}_{xx} \nabla^2 \psi = 0, \end{equation}

which has harmonic wave solutions $\propto e^{i (\omega t - k x - l y)}$ with the dispersion relation

(2.18)\begin{equation} \omega^{2} = B^{2}_{0} k^{2} \frac{k^{2} + l^{2}}{k^{2} + l^{2} + 1}. \end{equation}

These solutions are rotational Alfvèn waves. In the absence of rotation the non-dimensional $f ({=}1)$ disappears in the denominator of (2.18), and we recover the classical dispersion relation for Alfvèn waves. What is important in the present context, is that absolute values of the phase velocity of these waves lie in the interval $[0, |B_{0}|]$. The non-dimensional velocity of the modon $U$ is always larger than the value of the non-dimensional magnetic field in the far outer region $|B_{+}| = |\partial _{y} C_{+}|_{r \rightarrow \infty } = |\kappa _{+} U|$, as $|\kappa _{+}| <1$, and, hence, is larger than the maximum phase speed of the Alfvèn waves. This allows us to anticipate that the modon does not emit Alfvèn waves and can remain coherent for a long time.

3. Numerical investigation of the stability of magnetic modons

As previously mentioned, the MRSW (respectively MQG) model is structurally close to TRSW (respectively TQG). Our experience with TRSW and TQG modons (Lahaye et al. Reference Lahaye, Zeitlin and Dubos2020) shows that they may be subject to small-scale convective-type instabilities, which are known for vortex solutions in this model (Gouzien et al. Reference Gouzien, Lahaye, Zeitlin and Dubos2017). To check whether this is also the case with MQG modons, and also to see whether sharp gradients of the magnetic fields engender new instabilities, we performed a set of numerical simulations of the MQG equations using a doubly periodic pseudo-spectral code based on the Python library Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), with modon solutions as initial conditions. The typical resolution corresponded to a $512\times 512$ grid with a 3/2 dealiasing factor, and we used a split-explicit fourth-order Runge–Kutta method for temporal integration. Several control runs with higher resolutions were also performed. As is often done in 2DMHD simulations (e.g. Tobias et al. Reference Tobias, Diamond and Hughes2007), we applied the standard Newtonian viscosity and magnetic diffusivity to dissipate energy near the grid scale and ensure numerical stability, in the respective equations of the system (1.6), and with the same coefficient, i.e. we supposed the magnetic Prandtl number to be one. The viscosity and magnetic diffusivity coefficients were set to be $5\times 10^{-4}$ (i.e. $Re=2000$) in most runs. Simulations were run until $t=80$, a non-dimensional unit of time is the time taken by the dipole to travel a distance of one deformation radius (i.e. $U=1$). Notice that with the values of the viscosity which we used, the time of viscous decay at the grid scale, given by $\nu /{{\rm d}\kern0.7pt x}^2$, is about unity. Tests of the impact of the form of dissipation and of the Prandtl number on the results are discussed below in § 3.3. To ensure numerical stability and limit aliasing, the initial profiles of pressure and magnetic potential were smoothed with a Gaussian kernel of width twice the grid spacing. The main purpose of this operation was to smooth out the initial discontinuity in the magnetic field that exists when $\kappa _-\neq \kappa _+$. As will be shown below, in some configurations the modon leaves behind a wake of weak vorticity and magnetic field anomalies. As boundary conditions are periodic, this wake perturbs the modon as it re-enters the domain. In order to limit this effect, we repeated the simulations using a domain that was extended along the $x$-axis by a factor $11/8$, and adding a moving sponge layer (of width $3L/8$, with $L$ the size of the unexpanded domain) with linear damping to suppress the wake. The boundaries of the sponge layer had a $\tanh$ shape with width $0.2$. The maximum amplitude of the damping was 0.5, and its central position and propagation speed were computed using the $x$-position of the barycenter of kinetic energy, updated every $\Delta t=0.1$. To investigate the stability of the modon solutions, and their sensitivity to the values of dynamical parameters, we conducted a set of experiments with different values of the magnitude of the magnetic field. Six values of $\kappa _-$ and seven values of $\kappa _+$, in all combinations, were used: $\kappa _-=0, 0.02, 0.05, 0.1, 0.2, 0.5$ and $\kappa _+=-0.1, -0.05, 0, 0.02, 0.05, 0.1, 0.2$. Note that configurations with large $|\kappa _+-\kappa _-|$ (and, thus, in the present context, negative $\kappa _+$), maximise the inner-outer magnetic field strength differentials and, therefore, the magnitude of the tangential discontinuity of the magnetic field at the seperatrix. Sensitivity to initial conditions was further tested by conducting a second set of experiments with the same values of $\kappa _{\pm }$, adding a small-scale noise superimposed onto the dipolar solution in the initial conditions. The properties of the noise are described in Appendix A. Overall, 70 runs are discussed in the following subsection, and about five times as many were performed in total (see §§ 3.2 and 3.3).

3.1. Initialization with the pure modon solutions

The overall result of these simulations is that in any combination of the values of $\kappa _{\pm }$ each of them should be moderate, typically not exceeding $0.1$ in absolute value, for the modon to keep its form. At higher values the above-mentioned small-scale instabilities arise, most often after a few tens of time units, and destroy the modon. We will exaggeratingly call the modons ‘stable’ in the former and ‘unstable’ in the latter cases. One should keep in mind that these denominations are empirical, being entirely based on the results of numerical simulations of long, but limited duration ($t_{max}=80$). The latter is, however, large compared with the typical eddy time scale of the modon (used for the non-dimensionalisation), which supports the fact that the stable modons are, at least, long-living coherent structures. We present and discuss below examples of stable versus unstable modons in several typical configurations. In all of the snapshots only a part of the domain of the simulations, which extends from $-5$ to $5$ in units of $R_{d}$ in both directions ($\pm$6.875 along the $x$-axis for the simulations with a moving sponge layer), is shown, to make the details of the evolution of the modon clearer. We started by benchmarking the numerical method with a simulation initialized with a ‘non-magnetic’ modon, with $\kappa _{\pm }\equiv 0$ and no magnetic field whatsoever. The result is that, as expected, this modon keeps moving without changing its form and is, thus, stable (not shown).

Modons with no external magnetic field ($\kappa _+=0$), such as the one presented in figure 2, are stable as long as $\kappa _-\leqslant 0.2$, meaning that they exactly conserve their form over the duration of the simulation, travelling at a nearly constant speed (progressively decreasing under the action of energy dissipation). To make the presentation concise, we do not show their evolution. On the contrary, an instability leading to destruction of the modon takes place for higher values of $\kappa _-$, or in the presence of a strong enough external magnetic field, as shown in figure 4 in the case $\kappa _-=\kappa _+=0.2$ (this simulation was extended in time compared with other simulations to show the destabilization of the modon). Figure 5 shows the ‘hollow’ modon (the same type as shown in figure 3a) with $\kappa _+=0.2, \kappa _-=0$, which is unstable with a much faster destabilization compared with the previous case (figure 4).

Figure 4. Evolution of a regular modon with $\kappa _{\pm } = 0.2$. $t = 0$, 80, 90 and 100 (from top to bottom). (a,d,g,j) Streamfunction; (b,e,h,k) vorticity and (c,f,i,l) magnetic potential. Colours are saturated (especially for the vorticity at initial time). Note the evolving labels on the $x$ axes reflecting the propagation of the modon in the periodic domain (values correspond to the $x$-distance travelled by the modon, which position is estimated based on the maximum kinetic energy).

Figure 5. Evolution of a ‘hollow’ modon with $\kappa _{+} = 0.2$. From top to bottom: $t = 0$, $40$ and $80$.

The instability develops as follows: first, vorticity anomalies in the form of narrow filaments arise in the vicinity of the separatrix and spread along it and beyond, forming a wake at the lee side of the modon (see, e.g. figure 4 at $t=80$ and $90$).

The rings of vorticity formed by this process shield the vortices forming the initial dipole, which leads to their separation and slowing down of the modon. With some delay, development of vorticity anomalies along the separatrix is accompanied by ejection of any magnetic field from the core of the vortex pair, still along the separatrix, cf. figure 6. This perturbation then triggers the instability that develops within the dipole starting from the vicinity of the saddle points at the front and rear of the dipole – i.e. at the intersection of the modon axis and outer separatrix. The scales involved in this process are small, as seen in the vorticity field, and are associated with increasing magnitude of the magnetic field anomalies manifesting itself in the increasing magnetic energy (spatial mean of $|\boldsymbol {\nabla } A|^2$). The evolution of the different components of energy, shown in figure 7, clearly exhibits this increase of magnetic energy, as well as a decrease of the potential and kinetic energies – the latter resulting from the enhanced dissipation associated with the production of small-scale filaments by the instability. This implies that the form of the dissipation implemented in the numerical scheme can have an impact upon the details of the evolution. We checked in the cases presented in the paper (corresponding to figures 4, 5, and also figures 11 and 13 discussed below) that numerical convergence is reached by running numerical simulations with double resolution while keeping all other parameters (including the value of the viscosity and $Pr_m$) constant. In all cases, the scenario of the destabilization, and the typical scales involved were not changed, while the time scale of destabilization was only marginally changed (by a few time units at most). These verifications confirm that the reported small-scale instability and the mechanism of destruction of the modon are not numerical artifacts. The impact of the form and magnitude of the dissipation is discussed below in § 3.3.

Figure 6. Destabilization of the modon as seen in the vorticity field and in the anomaly of magnetic potential, highlighting the formation of a sharp gradient in the azimuthal magnetic field in the vicinity of the separatrix and the subsequent destabilization of the modon triggered near the saddle points. (a–d) Regular modon ($\kappa _\pm =0.2$) at intermediate stages compared with figure 4(d–i). (e–h) Hollow modon ($\kappa _+=0.2, \kappa _-=0$) at intermediate stages compared with figure 5. Colourbars are identical in both rows.

Figure 7. Evolution of the total (black), kinetic (blue), potential (orange) and magnetic (green) energy anomaly (see (1.8)), for two unstable cases ($\kappa _+=\kappa _-=0.2$, dash-dotted, and $\kappa _+=0.2, \kappa _-=0$, dashed), and one stable case ($\kappa _+=0, \kappa _-=0.2$, solid). Note the different ranges of the $y$-axis in the two frames.

3.2. Empirical condition for modon stability and impact of initial perturbations

To check the robustness of the modon solutions, we repeated the simulations by superimposing a small-scale noise in the streamfunction and the magnetic potential onto the analytic modon solutions, and used thus perturbed solutions as initial conditions. The characteristics of the noise are given in Appendix A. The overall result of these simulations is that stable regular modons are sufficiently robust to withstand such perturbation and keep coherence for a long time.

We synthesize the time scales of destabilization of the different modon configurations in figure 8. The estimate for this time scale was obtained considering the change (increase) of the slope of kinetic energy dissipation as a function of time, as triggering of the instability is associated with an enhanced decrease of kinetic (and potential) energy (cf. figures 7 and 9). Corrections were made for cases in which the instability only started developing at the end of the simulation, setting $t_{inst}=t_{max}=80$ (this concerns the cases with $(\kappa _+,\kappa _-)=(-0.1, 0.1)$ and $(0.2, 0.2)$). Hollow modons with $\kappa _+ \leqslant 0.1$, which includes the one shown in figure 3(a), are stable (not shown). Likewise, regular modons with $\kappa _\pm \leqslant 0.1$ are stable (not shown), while they are (weakly) unstable for $\kappa _\pm =0.2$ (as visible in figure 4).

Figure 8. Time scale of development of the instability (colours indicate non-dimensional time units) as a function of $\kappa _+$ (columns) and $\kappa _-$ (rows) for the simulations initialized with modon solution. White: stability. Colours: instability, with smaller (darker) to larger (lighter) times of development of the instability. Modons that turn unstable when some noise is superimposed on the initial condition are indicated by ‘P’ labels.

Figure 9. Same as figure 7 but for modons with $\kappa _-=0.2$ and two values of $\kappa _+$, with and without initial perturbation. Configuration $\kappa _\pm =0.2$ without perturbation (dashed) corresponds to the dash-dotted lines of figure 7.

In very few cases, as reported in figure 8, the addition of a perturbation can destabilize the otherwise stable modons (‘P’ labels in the figure). This occurs for modons with dynamical parameters ($\kappa _+, \kappa _-$) in the neighbourhood of cases unstable without initial perturbation. Time evolution of the energy for the regular modon with $\kappa _\pm =0.2$, which is unstable without perturbation (see figure 4), and a singular modon with the same $\kappa _-$ but $\kappa _+=0.1$ – a configuration that turns unstable if initially perturbed – is given in figure 9. It shows that the initial perturbation tends to accelerate the destabilization of the structure.

We can attempt a general conclusion that the greater the values of $\kappa _+$ and/or $\kappa _-$ (in absolute value), the more unstable is the dipole. However, from the simulations with $\kappa _+=0.2$ we infer that the instability develops faster for hollow modons (without internal magnetic field) as compared with the regular modons, while modons without an external magnetic field were found stable for $\kappa _-$ up to 0.20.

We believe this could be associated with the tangential magnetic discontinuity that exists for $\kappa _+\neq \kappa _-$, with the amplitude proportional to the difference between the two, since the initial perturbation seems to develop first in this region. Indeed, in the case of the regular modon, the vortex filaments along the separatrix starts developing only after a peak in the radial direction emerges in the magnetic streamfunction – meaning a sharp gradient of the magnetic field – at the same location (see figure 6a–d).

Thus, there seems to be a competition between the destabilization effects of the magnitudes of $\kappa _-$ and $\kappa _+$ taken separately (for instance, experiments with $\kappa _-=0.5$ exhibited instability at very early times), and of their difference. The details of the underlying mechanisms need further investigation, which is beyond the scope of the present work.

Interestingly, if we superimpose the initial random magnetic noise onto the ‘non-magnetic’ modon with $\kappa _{\pm }\equiv 0$, this noise aggregates in two magnetic dipoles inside the vortex dipole before being homogeneized in nearly circular patches within each pole. This can be understood from the dynamical equations in the limit of a weak magnetic perturbation $A'$: its action on the PV evolution (1.6a) becomes negligible and the magnetic field follows (1.6b), with an additional diffusion term. Magnetic perturbations are then advected along the streamlines and aggregate and/or get diffused. Then, the modon deflects slightly from the rectilinear trajectory, although it keeps its coherence, as follows from figure 10. This deflection is, apparently, due to asymmetry of the ‘insider’ magnetic dipole with respect to reflexions in $y$, resulting from asymmetries in the initial noise.

Figure 10. Evolution of a non-magnetic modon with $\kappa _{\pm } = 0$ under the influence of initial noise. From top to bottom: $t =0, 10$ and $30$. By the end of the simulation, the $y$-position of the modon is about $+1$. To mitigate the impact of diffusivity, as this case is stable, a lower viscosity ($\nu =2\times 10^{-4}$) is used in this simulation.

We should emphasize that the deflection of the modon from its initial trajectory by the residues of the initial perturbation tends to destabilize the modon because its direction of propagation deviates from that of the ambient magnetic field. This could explain the differences in behaviour between the modons with and without initial perturbation reported above, and also that hollow modons are more unstable than their ‘isolated’ counterparts.

3.3. Impact of the form of dissipation and of the Prandtl number

All the results presented above were obtained from numerical simulations using a viscosity term in the PV equation (1.6a), of the form $\nu \nabla ^2(\nabla ^2\psi )$ and diffusion term in the magnetic potential equation (1.6b) of the form $({\nu }/{Pr_m})\nabla ^2 A$, with $\nu =5\times 10^{-4}$ and a magnetic Prandl number $Pr_m=1$. As can be noticed in the behaviour of the kinetic energy and, to a lesser extent, in the magnetic energy in the stable cases (see, e.g. figure 7 with $\kappa _-=0.2, \kappa _+=0$), the loss of energy is substantial over the duration of the simulation. Moreover, as stated previously, the destabilisation of the modon leads to formation of small scales (vorticity filaments and magnetic sheets), meaning that the development and saturation of the instability are likely to be impacted by the dynamics at small scales and, therefore, by the dissipation. To the best of our knowledge, there is no physically based closure scheme for the Reynolds-averaged terms (parametrization) in 2-D or shallow water MHD, which is why we used the simple Newtonian viscosity in the above simulations. To check how sensitive our results are to the form and magnitude of the dissipation, we ran several sets of simulations changing the value $\nu$, but also the form of the viscosity/diffusivity, passing from the standard to hyper-viscosity and changing the magnetic Prandtl number. Recalling the link between the MQG and MRSW models, it would be possible to adapt the dissipation parameterizations developed for the RSW model, (e.g. Gilbert, Riedinger & Thuburn Reference Gilbert, Riedinger and Thuburn2014), but we postpone this development to a future investigation of the modons in the MRSW framework.

Let us first discuss the impact of the value of the viscosity coefficient at fixed $Pr_m=1$. To this aim, we repeated the whole set of numerical simulations (without sponge layers) with a smaller value of $\nu =2\times 10^{-4}$, and made additional runs for several choices of $(\kappa _+, \kappa _-)$ with even stronger viscosity, $\nu =10^{-3}$, or weaker viscosity, $\nu =10^{-4}$. The results show that some modons that appeared stable over the duration of the simulation with $\nu =5\times 10^{-4}$ turn to be unstable at $\nu =2\times 10^{-4}$. This concerns one pair of values of $(\kappa _+,\kappa _-)$ in the case without superimposed noise in the initial condition, and five of them in the perturbed case. However, in most of thus obtained unstable cases the results of the simulations do not appear to be well resolved because an aliasing occurs, as revealed by inspection of the corresponding spectra (not shown). This was also the case in the simulations with even weaker viscosity $\nu =10^{-4}$, while the simulations with stronger viscosity $\nu =10^{-3}$ were suffering from a too rapid dissipation of the modon, smoothing out any growing instability in most investigated cases.

This is why we conducted a set of simulations – for all previously exploited values of $\kappa _+$ and $\kappa _-$, and both with perturbed and unperturbed initial conditions – using a second-order hyper-viscosity (and hyper-diffusivity) instead of the standard viscosity (terms of the form $\nu _2\nabla ^4(\nabla ^2\varPsi )$ and $\nu _2/Pr_m\nabla ^4(A)$ in the PV and magnetic stream function equations, respectively). This closure scheme, although not based on physical arguments, as usual, allows for a sharper cutoff in the wavenumber space and for an extended range of resolved scales, while preserving numerical stability. We used two different values of hyper-viscosity/diffusivity: $\nu _2=10^{-6}$ and $10^{-7}$. As with the Newtonian dissipation, simulations with smaller values ($\nu _2=10^{-8}$) turned out to be not properly resolved. To provide some guidelines, we give here the inverse (hyper-)viscous damping time scale and the cutoff length scale associated with this dissipation. The former is given by $\nu _p/{{\rm d}\kern0.7pt x}^{2p}$, where $p$ is the order of (hyper-)viscosity ($p=1$ for the standard Laplacian viscosity). The latter is based on the typical CFL timestep ${{\rm d}\kern0.7pt x}/U$, and is given by $\delta x^c = (\nu \,{{\rm d}\kern0.7pt x}/U)^{1/2p}$, where we take $U=2$, the typical non-dimensional value of velocity at the centre axis of the modon. The regular viscosity $\nu _1=5\times 10^{-4}$ gives an inverse viscous damping time scale equal to $1.3$, while the hyper-viscosity gives 7 and $0.7$ for $\nu _2=10^{-6}$ and $10^{-7}$, respectively. For the cutoff length scale, we have $\delta x^c=2\times 10^{-3}$ for $\nu _1=5\times 10^{-4}$, against $9\times 10^{-3}$ and $5\times 10^{-3}$ for $\nu _2=10^{-6}$ and $10^{-7}$, respectively.

These simulations gave qualitatively similar results regarding the range of $\kappa _\pm$ values for which the modons are stable or unstable, and a similar scenario of destabilization, although they produce sharper gradients of vorticity and magnetic field, as could be expected. The destabilization of the regular modon with $\kappa _\pm =0.2$ with initial perturbation by the noise is shown in figure 11. In this case, the typical time scale of instability is slightly larger compared with the simulation with standard viscosity. It gets shorter with $\nu _2=10^{-7}$, as shown in figure 12. As visible in this figure, the instability has a more rapid development (sharper decrease of the kinetic and potential energies, and increase of magnetic energy anomaly) once initiated, as compared with the standard viscosity case, favoured by a very small energy dissipation before the onset of the instability, especially for $\nu _2=10^{-7}$. The development of the instability is associated with a deviation of the modon from the axis $y=0$. This deviation is already present in some of the simulations with standard viscosity and diffusion, but is greatly enhanced when an initial noise is present – as mentioned above. Our understanding is that this deviation follows from the onset of the instability (as shown in figure 5), which is asymmetric. (The symmetry breakdown in the configuration initialized without noise is ‘enabled’ by discretization errors. In the real world, it could be triggered by a small-scale noise.)

Figure 11. Evolution of the regular modon with $\kappa _{\pm } = 0.2$, similar to the one shown in figure 4 but with initial noise and hyper-viscosity/diffusivity with $\nu _2=10^{-6}$. From top to bottom: $t = 45$, $50$ and 55.

Figure 12. Comparison of the evolution of the total (black), kinetic (blue), potential (orange) energies and magnetic (green) energy anomaly for the unstable regular modon with $\kappa _\pm =0.2$, and perturbed by the noise, for different dissipation schemes: standard viscosity with $\nu _1=5\times 10^{-4}$ (continuous), hyperviscosity $\nu _2=10^{-6}$ (dashed) and $\nu _2=10^{-7}$ (dash-dotted).

A global comparison of the typical time scales of destabilization of the modons is given in Appendix B. A few more configurations turned to destabilize if this form of dissipation is used, compared with the standard viscosity with $\nu _1=5\times 10^{-4}$: $\kappa _+, \kappa _- = (0.05, 0.2)$ and $(0.1, 0.1)$ for $\nu _2=10^{-6}$ and six more for $\nu _2=10^{-7}$ (see Appendix B). Surprisingly, and related to our previous discussion about the criteria for instability, it appears that the absolute value of $\kappa _\pm$ seems to be more important than the magnitude of the discontinuity: for instance, the regular modons with $\kappa _\pm =0.1$ and $0.2$ are more unstable than their singular counterparts with $\kappa _\mp <\kappa _\pm$ (except for negative $\kappa _+$). This highlights the impact of small-scale dynamics near the viscous dissipation scale on the development of the instability. Its detailed investigation is, again, out of the scope of the present paper. In all configurations, modons with $\kappa _-=0.5$ were destabilizing on a very short time scale, as in the simulations with standard viscosity (performed only in the unextended domain and without a moving sponge layer).

To investigate the impact of the magnetic Prandtl number, we ran additional simulations with $Pr_m=0.1$ and $Pr_m=10$ and compared them to the reference simulations with $Pr_m=1$. We investigated configurations with $\kappa _\pm =0$ and $0.2$ (all four cases), both with the standard Laplacian viscosity $\nu _1=5\times 10^{-4}$ and the second-order hyper-viscosity $\nu _2=10^{-6}$. Random perturbations were superimposed on the initial condition. In general, we observed a destabilizing influence of $Pr_m>1$ (and vice versa). Indeed, the case with $\kappa _-=0.2, \kappa _+=0$, which is stable for $Pr_m=1$, is unstable for $Pr_m=10$ with the standard viscosity and magnetic diffusivity. Likewise, with order 2 viscosity/diffusivity, the modons with $\kappa _+=0.2$ are unstable for $Pr_m\geqslant 1$ and stable for $Pr_m=0.1$, with faster destabilization for $Pr_m=10$ compared with $Pr_m=1$. Figure 13 shows the destabilization of the regular modon with $\kappa _\pm =0.2$ and $Pr_m=10$ (and $\nu _2=10^{-6}$) for comparison with figure 11, where $Pr_m=1$. It exhibits finer patterns (in particular in the magnetic field) and a more rapid development of the instability: all fields are shown at earlier times in the simulation, while the stages of destabilization are similar to the those shown in figure 11(d–i).

Figure 13. Evolution of the regular modon with $\kappa _{\pm } = 0.2$, similar to the one shown in figure 11 but with $Pr_m=10$. From top to bottom: $t = 30$ and 35.