© The Authors 2026

1. Introduction

The magneto-rotational instability (MRI) is a crucial mechanism for sustaining magneto-hydrodynamic (MHD) turbulence in accretion disks, leading to enhanced angular momentum transport (Shakura & Sunyaev 1973; Balbus & Hawley 1998). The instability arises from the interaction between a weak magnetic field and a sheared background flow. Its significance in accretion disk modelling was first recognized in the early 1990s by Balbus & Hawley (1991), building on earlier work by Chandrasekhar (1960) and Velikhov (1959). The relevance lies in the fact that, while a differentially rotating disk is hydrodynamically stable according to the Rayleigh stability criterion (Rayleigh 1917), the presence of a magnetic field destabilizes any configuration in which the angular velocity profile decreases radially outwards.

Due to the general nature of the instability and its rapid growth rate, the MRI is often considered a key mechanism in the context of compact binary mergers (Shibata 2015; Radice & Hawke 2024; Abac et al. 2025). It is frequently invoked to explain the amplification of magnetic field strengths from the 10⁸–10¹¹ G typically observed in old neutron star systems (Lorimer 2008) to the ‘magnetar-level’ range of 10¹⁵–10¹⁶ G, hence aiding in the collimation of relativistic jets and the production of associated gamma-ray bursts (Ruiz et al. 2016; Ciolfi et al. 2019; Ciolfi 2020; Combi & Siegel 2023; Mattia et al. 2023; Jiang et al. 2025; Kalinani et al. 2025). Other mechanisms are also believed to play a role: magnetic winding, driven by differential rotation, is expected to linearly amplify the toroidal magnetic field, while the Kelvin-Helmholtz instability (KHI)–which develops in the shearing layer at the slip-line between merging neutron stars–has been shown to be highly efficient in amplifying the magnetic field in the core of the remnant (Price & Rosswog 2006; Anderson et al. 2008; Kiuchi et al. 2015, 2018).

Nonetheless, the MRI is expected to play a crucial role in this scenario (Duez et al. 2006; Siegel et al. 2013; Aguilera-Miret et al. 2023; Miravet-Tenés et al. 2022), as magnetic winding leads to linear growth; the KHI is believed to be less efficient in the envelope of the remnant. Recent evidence also suggests that an αΩ-type dynamo operates in the post-merger remnant, contributing to magnetic field amplification at large scales (Kiuchi et al. 2024; Most 2023). Given its characteristics, the MRI is the natural candidate for sustaining turbulence development in the remnant’s envelope and inducing the large-scale dynamo. A key feature of the MRI is the presence of a critical wavelength, λ_MRI ≈ 2πv_A/Ω, which scales with the Alfvén speed and corresponds to the wavelength at which the instability’s growth rate is maximized. As a result, discussions of the MRI in the context of binary mergers often focus on this fastest-growing mode (Siegel et al. 2013; Kiuchi et al. 2018; Ciolfi et al. 2019; Miravet-Tenés & Pessah 2025). Demonstrating that this critical wavelength is adequately resolved in simulations ensures that the simulation has sufficient resolution to capture the underlying dynamics. Notably, although this criterion originates from linear MRI theory, it has proven useful for empirically assessing the non-linear MRI development, as demonstrated by convergence studies in non-linear disk simulations (e.g. Hawley et al. 2013, 2011).

However, it is important to remember that the MRI is derived under certain assumptions, which raises questions about the validity of its predictions when those assumptions are violated. For instance, Celora et al. (2023) demonstrate that the MRI phenomenology is highly sensitive to the assumption of global axisymmetry. Despite the MRI being based on a plane-wave analysis–suggesting that the instability is local and independent of the large-scale conditions–the assumption of global axisymmetry is hard-wired in the derivation. In fact, Celora et al. (2023) show that relaxing this assumption of global symmetry results in a drastic change in the instability’s behaviour: the magnetic field acts to suppress or slow the unstable modes in an anisotropic manner. In this sense, a word of caution is in order.

Another important piece of the puzzle comes from the long-lived binary neutron star merger simulations presented in Aguilera-Miret et al. (2023, 2024). These simulations failed to provide clear evidence of poloidal growth after the vigorous amplification due to the Kelvin-Helmholtz unstable phase, despite the authors demonstrating that they could adequately resolve the supposed wavelength of the fastest-growing mode throughout the simulation. It is also worth mentioning that such simulations employ large-eddy simulation (LES) schemes (Viganò et al. 2019a; Carrasco et al. 2020; Viganò et al. 2020) designed to resolve the small-scale dynamics and its impact on magnetic field amplification–see Celora et al. (2024a,b, 2021) for recent progress towards the first covariant formulation for relativistic LESs (Abac et al. 2025).

While this may initially seem surprising, it is relatively easy to come up with a possible explanation. The simulations show a magnetic field configuration that remains turbulent even at late times (t ≳ 200 ms after merger), characterized by a small-scale dominated, randomly oriented field. In this scenario, the perturbations will likely undergo a random-walk process, which must be taken into account when calculating the effective Alfvén speed (Kiuchi et al. 2018; Palenzuela et al. 2022). As a consequence, the estimate for the Alfvén velocity is significantly reduced, which in turn lowers the wavelength of the fastest-growing mode, making the instability even more challenging to resolve numerically. Moreover, in contrast with protoplanetary disks, in such a turbulent highly rotating scenario, the advective timescales can easily be faster than or comparable to the growth timescales, thus challenging the very basic assumption about having a clear background.

The goal of this work was to examine more deeply the insights provided by the above simulations and, for the first time, explore how MRI conditions are influenced by a complex magnetic field. The standard MRI derivation in fact assumes a weakly magnetized fluid configuration, where the magnetic field is regular enough that we can neglect magnetic field gradients. In this respect, a few works have explored different assumptions: for instance, Doğan & Pekünlü (2012) considered a diamagnetic Ohm’s law but still with relatively simple configurations. In this study we focused in particular on the original MRI and the spectrally unstable axisymmetric modes. Notice also that depending on the net magnetic flux, non-axisymmetric modes may be more relevant to the MRI dynamics and developments (see e.g. Balbus & Hawley 1992; Rincon et al. 2007; Squire & Bhattacharjee 2014; Gogichaishvili et al. 2017; Held & Mamatsashvili 2022). We leave the exploration of the impact of complex field geometries on the development of transiently growing non-axisymmetric modes to future studies. We also present a comprehensive diagnostic of the generalized instability in the post-merger environment, achieved by directly post-processing the numerical output from said simulations. This is evidently motivated by our interest in understanding the fineprints of the magnetic field amplification in the post-merger environment, though the results we present are also of broader significance. As we demonstrate, magnetic gradients can significantly affect the instability, potentially suppressing it entirely if they are large enough. The relevance of this extends beyond just post-merger environments.

The paper is laid out as follows: In Sect. 2 we derive the generalized MRI, first discussing the analytical set-up in Sect. 2.1 and then focusing on the (most relevant) case with radial magnetic field gradients in Sect. 2.2. In Sect. 3 we discuss the generalized instability results for two magnetic field toy models inspired by accretion disk models. Finally, we present a detailed diagnostic of the generalized instability in a compact binary merger simulation in Sect. 4, and draw our conclusions in Sect. 5.

2. Impact of magnetic gradients on the axisymmetric MRI

As a first step, we performed a linear perturbation analysis on a magnetized fluid configuration, considering gradients in both the angular velocity and the magnetic field. In particular, we focused on the axisymmetric MRI modes driving the formation of so-called channel modes (Goodman & Xu 1994) and eventually leading to a magnetized turbulent state as these are destroyed by parasitic instabilities (see Miravet-Tenés & Pessah 2025, and references within). First, we verified that the results of the standard MRI are recovered when the magnetic field is constant. Then, we extended the analysis by assuming non-vanishing magnetic field gradients, which modify the characteristic timescale and wavelength of the instability.

2.1. Analytical set-up and WKB expansion

We considered a plane-wave analysis of perturbations on top of a non-homogeneous background fluid, which allowed us to account for the impact of differential rotation (and entropy gradients in the most general case) on the time-evolution and propagation of MHD modes. Effectively this means that every quantity A is split into background contribution plus fluctuations, and the latter are written in terms of a Wentzel–Kramers–Brillouin (WKB)-like expansion of the form

$$\begin{array}{c}\hfill A=a+\delta a=a+\overline{\delta }\left({\displaystyle \sum _{q=0}}{\epsilon }^q{a}_q\right)e^{i\theta (x)/\epsilon },\end{array}$$(1)

where x stands for the relevant (spatial) coordinates describing the system. The book-keeping parameter $\overline{\delta }$ measures the relative magnitude of the background versus perturbations, while ε ≈ λ/L is the ratio between λ, the typical wavelength of the waves, and L, which is the typical length-scale over which the wave amplitude, polarization, and wavelength vary (Anile 1989; Thorne & Blandford 2017). Using this expansion for each quantity entering the MHD equations, the terms of order $\mathcal{O}({\overline{\delta }}^0,{\epsilon }^0)$ are background terms, while those of order $\mathcal{O}({\overline{\delta }}^1,{\epsilon }^0)$ represent linear perturbations. Along with this WKB-type ansatz the (mathematical) calculation is based on the assumption that some of the background gradients are ‘steep’ in the sense that ${\partial }_{x}a=\mathcal{O}({\overline{\delta }}^0,{\epsilon }^0)$, while others are not, in the sense that ${\partial }_{x}a=\mathcal{O}({\overline{\delta }}^0,{\epsilon }^1)$. Steep background gradients–relative to a fast varying quantity–will consequently affect the perturbation equations (as well as changing the notion of background). By expanding the phase θ(x)/ε = θ(0)/ε + k_xx + …, we can then work out a dispersion relation where such gradients will appear. See Celora et al. (2023) for more details of this multi-scale WKB expansion in the context of the MRI.

We next observed that background gradients enter explicitly in the perturbation equations only when they appear as part of a non-linear term in the MHD equations (e.g. A∂_xB). Since the MHD equations are linear in the time-derivatives, temporal gradients do not explicitly modify the perturbation equations or the dispersion relation. In contrast, spatial gradients–both of the magnetic field (which is neglected in the standard derivation) and of the background velocity profile (which is considered in the standard derivation)–do play a direct role. This does not imply, however, that temporal gradients are irrelevant. The assumption that the background is stationary, or at most evolving slowly, is crucial; without this assumption, the separation of background and perturbation would not make sense. Assuming that the background evolves on a slow timescale T, any unstable mode must grow on a much shorter timescale τ ≪ T (or at least τ ≲ T) to be relevant. If this condition is not satisfied, the instability would not have sufficient time to grow before the background changes.

2.2. Linear analysis of the MHD equations

We next derived the generalized instability in the presence of a non-homogeneous magnetic field. We focused on the impact of adding radial magnetic field gradients since, as we will see shortly, they can have a crucial impact on the instability results. Our starting point was the non-relativistic ideal MHD equations

$$\begin{array}{cc}& {\partial }_t\rho +\mathrm{\nabla }·(\rho \mathbf{v})=0,\hfill \end{array}$$(2a)

$$\begin{array}{cc}& {\partial }_t\mathbf{v}+\mathbf{v}·\mathrm{\nabla }\mathbf{v}=-\frac{1}{\rho }\mathrm{\nabla }P-\mathrm{\nabla }\mathrm{\Phi }+\frac{1}{4\pi \rho }[\mathbf{B}·\mathrm{\nabla }\mathbf{B}-\mathrm{\nabla }\left(\frac{B^2}{2}\right)],\hfill \end{array}$$(2b)

$$\begin{array}{cc}& {\partial }_t\mathbf{B}=\mathrm{\nabla }\times (\mathbf{v}\times \mathbf{B}),\hfill \end{array}$$(2c)

where ρ,  v,  P,  Φ, and B indicate the mass density, fluid velocity, thermodynamic pressure, gravitational potential and magnetic field, respectively. We then performed the analysis using the local co-rotating frame construction presented in Celora et al. (2023), who explicitly show that the axisymmetric MRI can be reformulated this way. In short, this means setting up (and working with) a local Cartesian frame (comoving with the background flow) such that the x-axis always points in the radial direction, the y-axis in the azimuthal one, and the z-axis is aligned with the global axis of rotation. Clearly, such a construction is not necessarily possible in absence of a global sense of rotation (and hence a notion of axis of symmetry). In this regard, while this formulation is equivalent to the original, it has the advantage of explicitly highlighting the key role that axisymmetry plays in the result.

We started by assuming a circular velocity profile, $\mathbf{v}=\mathrm{\Omega }(R)R\widehat{\phi }$, and then considered an observer co-orbiting with the fluid on some given orbit, R₀. The observer velocity is then ${\mathbf{v}}_{\mathit{obs}}=\mathrm{\Omega }(R_0)R\widehat{\phi }={\mathrm{\Omega }}_{0}R\widehat{\phi }$ while the background velocity as measured in the co-rotating frame is¹

$$\begin{array}{cc}\hfill \mathbf{v}\prime =(\mathrm{\Omega }-{\mathrm{\Omega }}_0)& R\widehat{\phi }=\frac{\partial \mathrm{\Omega }}{\partial R}{|}_{R_0}(R_0+x)x\widehat{\phi }\hfill \\ \hfill & ⟹\mathbf{v}\prime =S_{0}x\widehat{y}+\mathcal{O}(x^2),S_0=\frac{\partial \mathrm{\Omega }}{\partial logR}{|}_{R_0}.\hfill \end{array}$$(3)

As a result, when writing the linearized MHD equations in the co-rotating frame, the relevant background velocity gradients to consider are ∂_iv_j = S₀δ_ixδ_jy. When it comes to the magnetic field, the one measured by the co-rotating observer is the same as in the inertial frame, as one can easily see using Lorentz transformations and taking the Newtonian limit. In essence, we assumed the background field to satisfy

$$\begin{array}{c}\hfill \mathbf{B}=\mathbf{B}(x)=B^y(x)\widehat{y}+B^z(x)\widehat{z},\end{array}$$(4)

namely we imposed only radial gradients in this section and assumed no radial field–noting that this is parallelled in the original MRI calculation to disentangle the instability from winding, and is later shown not to change the instability criteria.

We performed a plane-wave analysis with perturbations proportional to exp(ik_jx^j − iωt) under the following simplifications: consider only axisymmetric modes (k_y = 0), assume that the flow is incompressible (δρ = 0) and consider wavelengths small compared to the scale height, allowing us to neglect pressure gradients entirely². These simplifications are possible because the key ingredients to reproduce the MRI are a non-zero magnetic field and a sheared background. Finally, we considered the gravitational potential as being externally sourced, so that³δΦ = 0. After imposing all these conditions, the linearized equations for the perturbations can be written as

$$\begin{array}{cc}& k_x\delta v^x+k_z\delta v^z=0,\hfill \end{array}$$(5a)

$$\begin{array}{cc}& -i\omega \delta v^x-2{\mathrm{\Omega }}_0\delta v^y+({\partial }_x{v}_A^y)\delta v_A^y+({\partial }_x{v}_A^z)\delta v_A^z\hfill \\ \hfill & +i{k}_x(v_A^y\delta v_A^y+v_A^z\delta v_A^z)-i\delta v_A^x{k}_z{v}_A^z=0.\hfill \end{array}$$(5b)

$$\begin{array}{cc}& -i\omega \delta v^y+\frac{{\kappa }^2}{2\mathrm{\Omega }}\delta v^x-({\partial }_x{v}_A^y)\delta v_A^x-i{k}_z{v}_A^z\delta v_A^y=0,\hfill \end{array}$$(5c)

$$\begin{array}{cc}& -i\omega \delta v^z-({\partial }_x{v}_A^z)\delta v_A^x+i{k}_z{v}_A^y\delta v_A^y=0,\hfill \end{array}$$(5d)

$$\begin{array}{cc}& -i\omega \delta v_A^x-i{v}_A^z{k}_z\delta v^x=0,\hfill \end{array}$$(5e)

$$\begin{array}{cc}& -i\omega \delta v_A^y+({\partial }_x{v}_A^y)\delta v^x-i{v}_A^z{k}_z\delta v^y-s_0\delta v_A^x=0,\hfill \end{array}$$(5f)

$$\begin{array}{cc}& -i\omega \delta v_A^z+({\partial }_x{v}_A^z)\delta v^x-i{v}_A^z{k}_z\delta v^z=0,\hfill \end{array}$$(5g)

where we introduced the Alfvén velocity and the square of the epicyclic frequency:

$$\begin{array}{c}\hfill v_A^i=\frac{B^i}{\sqrt{4\pi \rho }},{\kappa }^2=2\mathrm{\Omega }\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}logR}+4{\mathrm{\Omega }}^2.\end{array}$$(6)

We have also highlighted in blue the additional terms originating from gradients in the background magnetic field.

For the next steps we followed the strategy adopted in the original derivation. First, we rewrote all perturbations in terms of δv^z, and introduced

$$\begin{array}{c}\hfill {\partial }_x{v}_A^z=a{\partial }_x{v}_A^y,\end{array}$$(7)

where a can be positive or negative, to slim down the notation. Doing this, and solving a two-equation linear system for δv_A^y and δv^y we arrive at

$$\begin{array}{cc}& \delta v^x=-\frac{k_z}{k_x}\delta v^z,\hfill \end{array}$$(8a)

$$\begin{array}{cc}& \delta v^y=i\left(\frac{\omega }{{\omega }^2-p^2}\right)\frac{k_z}{k_x}(\frac{{\kappa }^2}{2\mathrm{\Omega }}-S_0\frac{p^2}{{\omega }^2})\delta v^z,\hfill \end{array}$$(8b)

$$\begin{array}{cc}& \delta v_A^x=\frac{k_z}{k_x}\frac{p}{\omega }\delta v^z,\hfill \end{array}$$(8c)

$$\begin{array}{cc}& \delta v_A^y=-i\frac{k_z}{k_x}\frac{1}{{\omega }^2-p^2}[2p\mathrm{\Omega }-\frac{{\partial }_x{v}_A^y}{\omega }({\omega }^2-p^2)]\delta v^z\hfill \\ \hfill & =\frac{1}{v_A^y}\left(\frac{\omega }{k_z}-ia\frac{{\partial }_x{v}_A^y}{k_x}\frac{p}{\omega }\right)\delta v^z,\hfill \end{array}$$(8d)

$$\begin{array}{cc}& \delta v_A^z=-\left(\frac{p}{\omega }-ia\frac{{\partial }_x{v}_A^y}{\omega }\frac{k_z}{k_x}\right)\delta v^z,\hfill \end{array}$$(8e)

where we have introduced the shortcut p ≡ k_zv_A^z. Substituting these into the x-component of the Euler equation above (cf. Eq. (5b)) and rearranging terms, we arrive at the following dispersion relation:

$$\begin{array}{c}\hfill a_0{\omega }^4+a_2{\omega }^2+a_4=0,\end{array}$$(9)

where

$$\begin{array}{cc}& a_0=1-i\frac{{\partial }_x{v}_A^y}{v_A^y}\frac{k_x}{k_z^2},\hfill \\ \hfill & a_2=-[2p^2+{\kappa }^2\frac{k_z^2}{k^2}-ia({\partial }_x{v}_A^y)\frac{k_x{k}_z}{k^2}p-i\frac{({\partial }_x{v}_A^y)}{v_A^y}\frac{k_x}{k^2}{p}^2\hfill \\ \hfill & +a\frac{{({\partial }_x{v}_A^y)}^2}{v_A^y}\frac{k_z}{k^2}p-a^2{({\partial }_x{v}_A^y)}^2\frac{k_z^2}{k^2}],\hfill \\ \hfill & a_4=[p^4+2\mathrm{\Omega }{s}_0\frac{k_z^2}{k^2}{q}^2-ia({\partial }_x{v}_A^y)\frac{k_x{k}_z}{k^2}{p}^3+a\frac{{({\partial }_x{v}_A^y)}^2}{v_A^y}\frac{k_z}{k^2}{p}^3\hfill \\ \hfill & -a^2{({\partial }_x{v}_A^y)}^2\frac{k_z^2}{k^2}{p}^2].\hfill \end{array}$$

This dispersion relation is a polynomial of fourth order as in the standard MRI analysis, but with much more involved coefficients that depend not only on the angular velocity gradients but also on the magnetic field ones.

2.3. Recovering the standard MRI: Constant magnetic fields

As a first step, let us show that the usual MRI can be recovered from the more general dispersion relation Eq. (9). This provides a sanity check and also serves as a comparison for the more general results that we will discuss shortly. To recover the standard MRI we disregard the terms originating from the gradients of the magnetic field in Eq. (9) (i.e. the blue terms). We also considered purely vertical modes (that is set k_x = 0), which are the ones that grow the fastest. With these simplifications, we arrived at the following dispersion relation:

$$\begin{array}{c}\hfill {\omega }^4-(2p^2+{\kappa }^2){\omega }^2+p^2(p^2+2\mathrm{\Omega }{S}_0)=0,\end{array}$$(10)

which is a real, quartic algebraic equation for ω. From this equation it is straightforward to check that, in the limit p² ≫ κ² ∼ Ωs₀, the solutions are ω² = p² and there is no instability.

The solution for ω² of the full dispersion relation is given by

$$\begin{array}{c}\hfill {\omega }_{\pm }^2=\frac{1}{2}[2p^2+{\kappa }^2\pm \sqrt{{\kappa }^4+16{\mathrm{\Omega }}^2{p}^2}].\end{array}$$(11)

Notice that the ω_±²-solutions are real since the discriminant is always positive. However, if one of these solutions is negative, the corresponding two ω-solutions will be complex conjugates, and hence one is unstable.

To derive the instability criterion, one can then study the marginal case ω² ≈ 0, in which one can neglect the ω⁴ term and solve the dispersion relation as

$$\begin{array}{c}\hfill \frac{{\omega }^2}{p^2}=\frac{p^2+2\mathrm{\Omega }{S}_0}{2p^2+{\kappa }^2}.\end{array}$$(12)

We could restrict the analysis to ‘small wavenumbers’, p² ≪ κ² ∼ Ωs₀, as we have already found that in the opposite limiting case there is no instability. The denominator is always positive, since the Rayleigh criterion (ensuring hydrodynamic stability) is precisely κ² > 0 (Rayleigh 1917). It is then easy to check that the sign of ω² is dictated by the sign of S₀. In essence, we have recovered the standard MRI criterion, which states that, if the background angular velocity decreases radially outwards (i.e. S₀ < 0), then the system is unstable. It is worth emphasizing that unstable modes can still exist when p ∼ κ, although in this case the precise instability criterion is no longer wavenumber-independent.

The fastest vertical growing mode is found by considering that the ω₋²-solution is a function of k and solving for the stationary points, namely

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\omega }_{-}^2}{\mathit{dk}}=0⟹k_{\mathit{MRI}}={\left(\frac{B^z}{\sqrt{4\pi \rho }}\right)}^{-1}{({\mathrm{\Omega }}^2-\frac{{\kappa }^4}{16{\mathrm{\Omega }}^2})}^{1/2}=\frac{2\pi }{{\lambda }_{\mathit{MRI}}},\end{array}$$(13)

which allow us to define λ_MRI. We could also work out the corresponding growth rate:

$$\begin{array}{c}\hfill {\omega }_{\mathit{MRI}}^2={\omega }_{-}^2(k_{\mathit{MRI}})=-\frac{1}{4}{S}_0^2⟹{\tau }_{\mathit{MRI}}=\frac{2}{|S_0|},\end{array}$$(14)

which sets the timescale of the instability τ_MRI as the e-folding time of the fastest-growing mode. Although all these results are obviously well known in the literature (Balbus & Hawley 1998), it is worth stressing a point that is relevant in particular for numerical simulations: since the wavenumber of the fastest-growing mode scales with the inverse of the background magnetic field, the instability occurs at larger wavelengths and is easier to resolve when the background magnetic field is stronger (i.e. λ_MRI ∝ B^z). In contrast, the precise value of the magnetic field does not influence the instability criterion nor the growth rate of the fastest-growing mode.

Before we move on to the case with radial gradients, it is important to recall that a fundamental assumption of the linear analysis is that unstable wavelengths should be much smaller than the length-scales over which pressure, density and gravity fields vary. This is relevant because it implies that wavenumbers smaller than the minimum value k_min = 2π/L_disk, where, say, L_disk = min(R, H) and H is the scale-height while R is the disk size, cannot be excited in practice. The instability then operates only when the wavenumber exceeds this threshold k > k_min, which represents a consistency condition ensuring that the mode fits within the system. On the other hand, it is suppressed in regimes where p ≫ κ, implying that, for a given background magnetic field, there exist a maximum wavenumber k_max for which the MRI is active. The net effect, a well-established result, is that sufficiently strong magnetic fields suppress the instability by preventing both conditions from being simultaneously satisfied. Similarly, the growth timescale should be much smaller than the timescale of variation of the background field. While this likely poses no problems in scenarios like protoplanetary disks, we have to be very careful in the case of the fast-evolving dynamics of a binary merger remnant.

2.4. Extended MRI with radial gradients of the magnetic field

We solved the dispersion relation by retaining the terms that originate from background gradients in the magnetic field. In this case, the dispersion relation is still a quartic algebraic equation but no longer real since the coefficients are generically complex-valued. Nevertheless, we could consider purely vertical modes as in the previous section: by setting k_x = 0, the dispersion relation reduces to a real quartic algebraic equation again,

$$\begin{array}{c}\hfill {\omega }^4-[2p^2+{\kappa }^2+f]{\omega }^2+p^2[p^2+2\mathrm{\Omega }{S}_0+f]=0,\end{array}$$(15)

where all the information about the gradients of the magnetic field is contained in the new quantity defined as

$$\begin{array}{c}\hfill f\equiv ({\partial }_x{v}_A^y)({\partial }_x{v}_A^z)\frac{v_A^z}{v_A^y}-{({\partial }_x{v}_A^z)}^2.\end{array}$$(16)

To understand the impact of the new terms, we followed the same strategy adopted in the standard analysis. Solving Eq. (15) for ω² we arrive at

$$\begin{array}{c}\hfill {\omega }_{\pm }^2=\frac{1}{2}[2p^2+{\kappa }^2+f\pm \sqrt{{({\kappa }^2+f)}^2+16{\mathrm{\Omega }}^2{p}^2}]\end{array}$$(17)

and observe that, as before, the discriminant is positive: the ω²-solutions are real, implying that there is an instability if one of them is negative. Following the same procedure as in the previous subsection, we focused on the ω₋² solutions as this family contains the (potentially) unstable modes with larger growth rates. We could then compute the wavenumber of the fastest-growing mode and the corresponding maximum growth rate, namely

$$\begin{array}{cc}\hfill k_{\mathit{eMRI}}& ={\left(\frac{B^z}{\sqrt{4\pi \rho }}\right)}^{-1}{({\mathrm{\Omega }}^2-\frac{{({\kappa }^2+f)}^2}{16{\mathrm{\Omega }}^2})}^{1/2}=\frac{2\pi }{{\lambda }_{\mathit{eMRI}}},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill {\omega }_{\mathit{eMRI}}^2& =-\frac{1}{4}{(s_0-\frac{f}{2\mathrm{\Omega }})}^2⟹{\tau }_{\mathit{eMRI}}=\frac{2}{|s_0-f/2\mathrm{\Omega }|}.\hfill \end{array}$$(19)

Regarding the instability criterion, we first observe that the regime where p² ≫ f, κ² ∼ Ωs₀ leads again to the same equation as in the standard MRI case, namely

$$\begin{array}{c}\hfill {\omega }^4-2p^2{\omega }^2+p^4={({\omega }^2-p^2)}^2=0⟹\omega =\pm ({\mathbf{v}}_A·\mathbf{k}).\end{array}$$(20)

Therefore, in this regime there is no instability and the modes are essentially Alfvén waves⁴.

To derive the modified instability criterion, we again considered the marginal case, ω² ≈ 0, and restricted to the regime⁵p² ≪ f, κ²:

$$\begin{array}{c}\hfill \frac{{\omega }^2}{p^2}=\frac{p^2+2\mathrm{\Omega }{S}_0+f}{2p^2+{\kappa }^2+f}\simeq \frac{2\mathrm{\Omega }{S}_0+f}{{\kappa }^2+f}.\end{array}$$(21)

To better understand the impact of the additional gradients on the instability, we considered different choices that lead to increasing complexity in the analysis. In the first case we considered radial gradients only in the azimuthal magnetic field (i.e. ∂_xv_A^z = 0, ∂_xv_A^y ≠ 0). With this setting, the function f vanishes trivially, implying that the phenomenology is the same as in the standard MRI.

In the second case we included radial gradients only in the vertical magnetic field (i.e. ∂_xv_A^z ≠ 0, ∂_xv_A^y = 0). In this case, f = −(∂_xv_A^z)² < 0, and from Eq. (19) it is clear that the timescale of the instability is longer than in the standard MRI case. As for the instability criterion, from Eq. (21) we immediately obtain⁶

$$\begin{array}{cc}\hfill |f|<{\kappa }^2⟶{\left(\frac{\partial v_A^z}{\partial R}\right)}^2& <2\mathrm{\Omega }\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}logR}+4{\mathrm{\Omega }}^2\hfill \\ \hfill & ⟹{\left(\frac{\partial v_A^z}{\partial R}\right)}^2<{\mathrm{\Omega }}^2\text{for a Keplerian disk}.\hfill \end{array}$$(22)

In this set-up the instability window gets smaller, as magnetic field gradients slow down the instability and can even stabilize the system if they are large enough. Consequently, there is an increase in τ_eMRI, which means that the instability grows more slowly, compared to the standard case.

Finally, we considered the general case, where both magnetic field components have non-vanishing radial gradients (i.e. ∂_xv_A^z ≠ 0, ∂_xv_A^y ≠ 0). This case might be better understood by introducing the (signed) characteristic scales over which the background magnetic field varies:

$$\begin{array}{c}\hfill \frac{{\partial }_x{v}_A^y}{v_A^y}=\frac{1}{L_y},\frac{{\partial }_x{v}_A^z}{v_A^z}=\frac{1}{L_z}⟹f=-{({\partial }_x{v}_A^z)}^2(1-\frac{L_z}{L_y}).\end{array}$$(23)

This is convenient as it allows us to distinguish three scenarios. A first scenario is when L_y ≫ L_z, that is, when gradients in the azimuthal component are much smaller than those in the vertical field. In this situation we effectively get back to the phenomenology just discussed for the case where f = −(∂_xv_A^z)². A second scenario is when L_z ≫ L_y, namely vertical gradients are negligible when compared to azimuthal ones. In this case, we can rewrite f = (v_A^z)²/L_yL_z, and observe that its sign is dictated only by the relative orientation of the radial gradients in the two components. The third, most general scenario is the one in which the two lengths are comparable L_y ≈ L_z, and the sign of f is dictated by the relative orientations as well as the magnitudes of the gradients of the two magnetic field components.

After presenting all the possible scenarios, the generalized instability can be characterized as follows depending on the sign of the f-function:

1. If f > 0: the instability window is defined by f < 2Ω|S₀| (for a Keplerian disk f < 3Ω²), while λ_eMRI grows together with a decrease in τ_eMRI.

2. If f < 0: the instability window is given by |f|< κ² (for a Keplerian disk |f|< Ω²), while τ_eMRI grows. The wavelength λ_eMRI can either decrease or increase.

Although λ_eMRI and τ_eMRI can shift towards larger or smaller values depending on the sign of f, we find that radial gradients generically reduce the instability window and can make the system stable for sufficiently large values, of the order of 𝒪(Ω²) for a Keplerian disk.

Notice that we directly considered the case where S₀ < 0, as this is the regime where the instability appears in the standard MRI setting. It is interesting to notice that the case S₀ > 0 can be unstable in the presence of radial magnetic field gradients, provided that

$$\begin{array}{c}\hfill f<0,{\kappa }^2-4{\mathrm{\Omega }}^2\le |f|\le {\kappa }^2.\end{array}$$(24)

Note that, in the practical cases presented in the following sections, we have not found any region satisfying these quite restrictive conditions.

Finally, while it is natural to consider the regime p² ≪ f,  κ² to derive a simple instability criterion and to compare with the standard scenario, it is worth revisiting the issue of consistency with the basic assumption of linear analysis–namely, that unstable wavelengths should be smaller than the relevant background variation scales. We can distinguish two qualitatively different regimes. In absence of significant gradients, compatibility with the basic assumption of linear analysis together with the fact that the instability is suppressed when $p^2\gg {\kappa }^2$, implies that very strong magnetic fields can inhibit the instability (as discussed previously). In the presence of strong background gradients, the instability is suppressed when p² ≫ κ², f. In this case, if f ≳ κ², the relevant maximum wavenumber becomes k_max ≃ 1/L_B, where L_B is the characteristic length-scale of magnetic field variations. Here, consistency with the basic assumption of linear analysis does not lead to any additional condition depending on magnetic field strength. Ultimately, this simply reflects the fact that, in the presence of gradients, an additional characteristic frequency, $\sqrt{f}$, enters the problem, leading to a richer and more varied instability phenomenology.

2.5. MRI with generic gradients

We have just demonstrated that radial magnetic field gradients can significantly impact the MRI criterion and its phenomenology. This raises the natural question of whether we can extend the analysis to study the effect of, for example, vertical magnetic field gradients. In principle, we could aim to consider the most general setting, though this would inevitably complicate the resulting dispersion relation, making it less intuitive and harder to interpret. Instead, the strategy we adopted is to examine a series of different settings and analyse each of them individually. We omit a detailed analysis here (which can be found in Appendix A) and instead provide a summary of the results in Table 1. This is useful as it clearly highlights the central role that radial magnetic field gradients play in generalizing the standard MRI criteria to complex field configurations. As a result, we focus in the following sections on exploring the consequences of the generalized criterion we have just derived.

2.6. MRI-driven dynamos and axisymmetric MRI

Throughout this work, and in the preceding discussion of this section, we have focused on the influence of magnetic field gradients on axisymmetric, spectrally unstable MRI modes. These modes are particularly relevant in scenarios where, say, the accretion disk is threaded by a weak mean vertical magnetic field, such that there is a net poloidal magnetic flux. In such cases, axisymmetric MRI modes become unstable, giving rise to so-called channel modes, which are eventually disrupted by secondary parasitic instabilities leading to the magnetized turbulent state (Goodman & Xu 1994; Miravet-Tenés & Pessah 2025). It is important to note, however, that MRI-driven dynamos have also been extensively investigated in alternative set-ups, such as the toroidal magnetic flux and the zero-net-flux configuration (see e.g. Fromang & Papaloizou 2007; Kapyla & Korpi 2011; Herault et al. 2011; Mamatsashvili et al. 2013; Gogichaishvili et al. 2017; Held & Mamatsashvili 2022). These studies are designed to explore the (long-term) sustainability of a turbulent magnetized state in the absence of a mean poloidal magnetic field, as might be the case in certain accretion disk environments. In particular, it is known that a subcritical, self-sustaining dynamo process is possible in such configurations. This process fundamentally relies on the non-modal (transient) amplification of non-axisymmetric modes (Rincon et al. 2007; Mamatsashvili et al. 2013; Rincon 2019). Our present study does not address the influence of additional gradients on the development of such non-axisymmetric modes, which would require a more global approach than the local WKB analysis adopted here.

3. Application to analytical models

In this section we explore the phenomenology of the extended MRI criteria using simple analytical models for the magnetic field configuration, which are inspired from the planetary and accretion disks literature (e.g. Pelletier & Pudritz 1992).

3.1. Vertical magnetic field with radial power-law decay

We considered a simple magnetic field configuration with a purely vertical magnetic field with a constant and a decaying components, namely

$$\begin{array}{c}\hfill B^z=B_0+\frac{k}{R^n}=B_0+\frac{B_1}{x^n},x=\frac{R}{R^{\ast }},\end{array}$$(25)

where n sets the power-law decay and R^* is some reference radial distance. The redefinition in terms of x ensures that B₁ has the same dimensions as B₀, and also that the magnetic field at R^* is the same across all models with different values of n. In this configuration the characteristic scale of the magnetic field gradients reads

$$\begin{array}{c}\hfill L_z=-\frac{x{R}^{\ast }}{n}(\frac{B_0}{B_1}{x}^n+1)\approx -2\frac{R^{\ast }}{n},\end{array}$$(26)

where, in the second step, we have assumed: (i) that the magnitudes of the decaying and constant components are comparable, and (ii) that x ≈ 1. This expression shows that L_z is primarily set by R^* and scales inversely with the power-law exponent n–larger values of n correspond to steeper magnetic field gradients. In the limit where B₀ ≫ B₁, the radially varying component becomes negligible compared to the constant offset, and the characteristic scale becomes much larger, |L_z|≫1, accordingly. We assumed that the velocity field follows a Keplerian profile, Ω = c_ΩR^−3/2, and introduced the dimensionless parameter

$$\begin{array}{c}\hfill C=\frac{B_1^2{R}^{\ast }}{4\pi \rho c_{\mathrm{\Omega }}^2},\end{array}$$(27)

which quantifies the relative strength of the radially varying magnetic component with respect to the shear in the background flow. We considered values C ≲ 0.1, corresponding approximately to a situation in which the total magnetic energy is less than 10% of the kinetic energy, for which the magnetic feedback on the hydrodynamics can be safely neglected.

Considering the definition of f, Eq. (16), it is evident that f < 0, in which case the generalized instability criterion reads

$$\begin{array}{c}\hfill \frac{-f-{\kappa }^2}{{\mathrm{\Omega }}^2}=n^{2}C{x}^{1-2n}-1<0.\end{array}$$(28)

It is straightforward to compute the timescale of the instability τ_eMRI, Eq. (19), and compare it to the standard MRI result τ_MRI, Eq. (14):

$$\begin{array}{c}\hfill \frac{{\tau }_{\mathit{MRI}}}{{\tau }_{\text{eMRI}}}=\frac{|2\mathrm{\Omega }{s}_0-f|}{|2\mathrm{\Omega }{s}_0|}=|1-\frac{1}{3}{n}^{2}C{x}^{1-2n}|.\end{array}$$(29)

We also compared the wavelength of the fastest-growing mode, obtaining

$$\begin{array}{c}\hfill \frac{{\lambda }_{\mathit{MRI}}}{{\lambda }_{\mathit{eMRI}}}=\sqrt{|1-\frac{1}{15}[1-{(1-n^{2}C{x}^{1-2n})}^2]|}.\end{array}$$(30)

It is not surprising that B₀ does not enter into the formula for the generalized instability criterion nor that of the timescale of the instability, as these depend only on the gradients of the magnetic field. This is analogous to the standard MRI case, where the instability criterion makes no reference to the background magnetic field (Balbus & Hawley 1998). We also stress that we here focused on the wavelength and growth rate of the fastest-growing mode, as this provides a straightforward basis for comparison with the standard case. It is important to note, however, that this mode is most relevant only when its wavelength remains sufficiently smaller than the characteristic size of the system. If this condition is not satisfied, slower-growing modes with more suitable wavelengths may dominate the dynamics. Nevertheless, the comparison remains meaningful, since the function f, which captures the effect of radial gradients, enters the dispersion relation in a ‘wavenumber-independent’ manner (cf. Eq. (17)). As a result, its impact can be extrapolated to modes beyond the fastest-growing one.

It is illustrative to plot the quantities that we just derived considering different values of n and C, starting with the generalized instability criterion, displayed in the left panel of Fig. 1. As shown in the figure, the impact on all three plotted quantities increases with larger values of the constant C and steeper power-law decays–parameters that together determine the strength of the gradients. In particular, for n = 2, the gradients are never strong enough to suppress the instability (given the values of C considered), although they do lead to changes of a few percent in the characteristic wavelength and up to 20% in the growth timescale. This changes for higher values of n (and larger values of C), where the gradients become strong enough to suppress the instability in certain regions of the domain. Most notably, we observe a significantly stronger effect on the growth timescale, with reductions reaching up to 30% in the most extreme cases.

Fig. 1. Extended MRI analysis for the analytical model given by a vertical magnetic field with radial power-law decay in Eq. (25). Left: Generalized instability criterion. Middle: Comparison between the generalized and standard MRI for the wavelength of the fastest-growing mode. Right: Timescale of the instability. We have considered the values with n = 2,  4 and C = 0.01,  0.025,  0.1 in the model defined by Eq. (25). The wavelength and timescale are shown only when the instability is active according to the generalized criterion, a transition that is marked with black dots.

3.2. Vertical and toroidal magnetic field

A second case that is worth considering is motivated by the simplified analytical models of disks (Galli et al. 2006; Jafari 2019). Protoplanetary disks are an excellent observational probe of the magnetic field topology: polarization measurements of high-resolution millimetre observations (e.g. Beltrán et al. 2019; Maury et al. 2022; Huang et al. 2024) provide the magnetic field direction on the plane of the sky, which can be often approximated as a ‘hourglass’ configuration, i.e. a vertical field pinched on the mid-plane. Neglecting the turbulent contributions, which are actually likely relevant in several cases (Huang et al. 2025), the field can be approximated by non-vanishing vertical and azimuthal components.

To investigate the effects of non-vanishing radial gradients in both magnetic field components, we introduced two key quantities: the ratio (r) between the characteristic gradient scales in the vertical and azimuthal components, and a dimensionless parameter, $\stackrel{\sim }{C},$ which characterizes the strength of the radial gradient in the vertical Alfvén speed relative to the background shear, namely

$$\begin{array}{c}\hfill r=\frac{L_z}{L_y},\stackrel{\sim }{C}={\left({\partial }_R{v}_A^z\right)}^2\frac{{R^{\ast }}^3}{c_{\mathrm{\Omega }}^2}⟹f=\stackrel{\sim }{C}(r-1).\end{array}$$(31)

We also emphasize that, since the characteristic scales carry a sign (cf. Eq. (23)), their ratio can be positive when the two magnetic field components vary in the same radial sense (either both increasing or both decreasing), and negative when they vary in opposite senses. In the case of a vertical magnetic field obeying a radial power-law profile, as considered in the previous section, this quantity is related to the parameter C via $\stackrel{\sim }{C}\approx n^{2}C$. We considered the range $\stackrel{\sim }{C}\in [0,0.24]$, which corresponds to C ≈ 0.06 for n = 2. The ratios of the modified MRI characteristic scales to their standard counterparts are given by

$$\begin{array}{cc}\hfill \frac{{\lambda }_{\text{MRI}}}{{\lambda }_{\text{eMRI}}}& =\sqrt{|1-\frac{1}{15}[{\stackrel{\sim }{C}}^2{(r-1)}^2+2\stackrel{\sim }{C}(r-1)]|},\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill \frac{{\tau }_{\text{MRI}}}{{\tau }_{\text{eMRI}}}& =|1+\frac{1}{3}\stackrel{\sim }{C}(r-1)|.\hfill \end{array}$$(33)

The instability window and the corresponding ratios of characteristic scales are shown in Fig. 2. The origin of each panel corresponds to the standard MRI, located in a regime where magnetic field gradients are small compared to the velocity shear. A key feature is the line r = 1, where the radial gradient scales of the magnetic field in the two components are equal. Along this line, the ratio of scales is exactly one, regardless of the magnitude of the gradients, as seen in Eq. (32). Deviations from unity in the plotted ratios remain modest for small values of f, and become more pronounced with stronger gradients and/or significant mismatches between L_z and L_y. The impact of magnetic gradients on the instability window is particularly evident in the lower half-plane, where instability is suppressed at lower values of $\stackrel{\sim }{C}$, especially in the presence of strong mismatches. In contrast, the instability is never fully suppressed in the upper half-plane for the parameter values considered. The timescale ratio also displays a clear asymmetry: the standard MRI timescale is larger than the extended one in the upper half-plane–indicating that gradients enhance the growth of unstable modes–while it is smaller in the lower half-plane, where gradients reduce the growth rate. As in the previous example, the reduction remains moderate, with deviations in the lower half-plane typically limited to 20–30%. Finally, the wavelength ratio is close to one across most of the parameter space explored. However, for sufficiently strong gradients and mismatches, the characteristic wavelength becomes noticeably shorter in the upper half-plane.

Fig. 2. Extended MRI analysis for the analytical model given by a vertical and toroidal magnetic field. From left to right: Analytic instability window, ratio of the wavelength of the fastest-growing mode, and ratio of the timescale of the instability in the generalized case vs the standard MRI result. The wavelength and timescale of the fastest-growing mode are shown only where the instability is active according to the generalized criterion.

4. Application to a numerical simulation of binary neutron star remnants

In this section we investigate the impact of magnetic gradients on the generalized MRI phenomenology associated with a more extreme scenario: the long-lived, hyper-massive neutron star produced by a binary merger. The numerical solution was obtained by solving the Einstein equations coupled to the general relativistic MHD equations (Andersson et al. 2017, 2022) through the coalescence of a binary neutron star system. In particular, we considered the outcome of the simulations by Palenzuela et al. (2022), Aguilera-Miret et al. (2023) and Aguilera-Miret et al. (2024), which were able to at least partially resolve the small-scale dynamics and lead to a converging saturation of the magnetic energy, dominated by winding and small scales, starting from a much smaller magnetic field seed. These results were possible by a combination of high-resolution, high-order numerical schemes, and the use of advanced techniques to capture the sub-grid scales present during the turbulent regime. It is important to highlight that, while the analysis described in the previous sections was carried out within the non-relativistic framework, the simulation we considered was performed using full general relativity. This approach remains reasonable nonetheless, as the MRI is derived within the non-relativistic framework before being applied to the context of binary mergers.

4.1. Numerical set-up and analysis

The numerical simulation was performed using the code MHDUET, which is automatically generated by the platform SIMFLOWNY (Arbona et al. 2013, 2018) to run on the SAMRAI infrastructure (Hornung & Kohn 2002; Gunney & Anderson 2016) with parallelization and adaptive mesh refinement. It uses fourth-order-accurate operators for the spatial derivatives in the Einstein equations, supplemented with sixth-order Kreiss-Oliger dissipation; a high-resolution shock-capturing method for the fluid, based on the Lax-Friedrich flux splitting formula (Shu 1998) and the fifth-order reconstruction method MP5 (Suresh & Huynh 1997); a fourth-order Runge-Kutta scheme with sufficiently small time step Δt ≤ 0.4 Δx (where Δx is the grid spacing); and an efficient and accurate treatment of the refinement boundaries when sub-cycling in time (McCorquodale & Colella 2011; Mongwane 2015). A complete assessment of the accuracy and robustness of the numerical methods implemented can be found in Palenzuela et al. (2018), Viganò et al. (2019b). The binary neutron star is evolved in a cubic domain of size [ − 1228, 1228] km using eight refinement levels, each with twice the resolution of the previous. The finest level covers the region where the density is ρ ≥ 5 × 10¹² g cm⁻³. In this region, we achieve a maximum resolution of Δx_min = 60 m during the first 110 ms after the merger, which is decreased to Δx_min = 120 m later on.

For our post-process analysis, and due to quasi-axisymmetry of the solution at late times, it is more convenient to compute averages of the relevant fields along the azimuthal direction φ. Therefore, we decomposed a generic field U(R, z, φ), which can be either a scalar or a component of a tensor field, as the average $\overline{U}(R,z)$ plus a residual δU(R, z, φ). We represent 2D maps of these averaged fields as functions of the cylindrical coordinates (R, z). To be more specific, the averages were calculated using

$$\begin{array}{cc}\hfill \overline{U}(R,z)& =\frac{{\int }_0^{2\pi }U(R,\phi ,z){\chi }^{-1/2}Rd\phi }{{\int }_0^{2\pi }{\chi }^{-1/2}Rd\phi },\hfill \end{array}$$(34)

where the cylindrical coordinates are centred on the system centre of mass–the cylindrical axis z is the same as the one relative to the Cartesian coordinates used throughout the simulation (i.e. perpendicular to the orbital plane)–and χ = γ^−1/3 is the conformal factor of the spatial metric⁷γ_ij in the standard 3+1 decomposition (Mizuno & Rezzolla 2025; Palenzuela 2020; Gourgoulhon 2012). These azimuthal averages are performed considering only the high-density regions with ρ ≥ 5 × 10¹⁰ g cm⁻³, which contain the bulk of the remnant mass and are less sensitive to artefacts arising from the low inertia treatment (simulations use a floor value 6 × 10⁴ g cm⁻³ for density).

In Fig. 3 we plot the 2D maps of the magnitude of the various magnetic field components at t = 40 ms. In particular, we observe that the magnetic field components present large gradients, which suggests that they need to be considered to study the conditions for MRI. These gradients are in both vertical and radial directions and particularly evident for the azimuthal and vertical components. The variations seem more ordered close to the orbital plane, while they are more stochastic in the less dense regions, where also the azimuthal variations (averaged out by definition) are more important (see Appendix B for a detailed analysis). In Fig. 4 we show the angular velocity (left panel) and its (logarithmic) derivative (middle panel). As it is expected in hyper-massive neutron star remnants, the angular velocity decays both with the cylindrical radii and with the absolute value of the z coordinate (i.e. as we move away perpendicularly from the equatorial plane).

Fig. 3. 2D maps of the azimuthally averaged magnetic field magnitudes at t≃ 40 ms after merger. From left to right: 2D maps of the magnitude of the radial, azimuthal, and vertical component of the magnetic field. We see the different components present large gradients, which we take into account in the generalized MRI criterion.

Fig. 4. Instability window calculation at t≃ 40 ms after merger. From the angular frequency (Ω) we computed its gradient (s0) and mask those points where s0 > 0, namely those for which the standard MRI would be inactive. Then we evaluated the updated instability criterion, considering at once the two cases f > 0 and f < 0. In the right panel, blue points are unstable according to the updated MRI criterion, and red are those for which the magnetic gradients are large enough to switch off the instability; all coloured points are MRI-unstable according to the standard criterion.

4.2. Comparison of the standard versus extended MRI analysis

Here we describe in detail the procedure for calculating a new instability window using the extended MRI criteria, as well as the predictions for λ_eMRI and τ_eMRI. The required steps are represented in Fig. 4 for an illustrative time snapshot at t = 40 ms. First, using the angular velocity we could easily compute S₀ (cf. Eq. (3)) and we mask those points where S₀ > 0, namely those for which the standard MRI would be inactive. Then we evaluated the updated instability criterion, considering at once the two cases f > 0 and f < 0 (cf. Items and ). The points in the right panel that are blue are unstable according to the new generalized MRI criterion, whereas the red ones are those for which the magnetic gradients are large enough to shut down the instability. All coloured points in this panel would be MRI-unstable according to the standard criterion. Notice that including the new branch of the instability S₀ ≥ 0 has no impact at all in this plot, since our solution never satisfies the quite restrictive conditions κ² − 4Ω² ≤ |f|≤κ².

We next compared the predictions for λ_eMRI and λ_MRI, displayed in Fig. 5. In the standard MRI case (upper plot), we observe that λ_MRI can reach values of 10³ − 10⁴ m, especially at large radii, R ≳ 15 km, while it is considerably smaller, λ_MRI ∼ 10 m, at R ≲ 15 km. The extended MRI case is considered in the lower panel. First of all, one can note that the potentially unstable regions are reduced compared to the standard criterion. The regions excluded by the criteria extension are mainly the ones having the largest λ_MRI. The potentially unstable region is mainly composed by a pseudo-torus, with a maximum ∼20km in height and extended from ∼7 to ∼17km in the orbital plane. In these regions, λ_eMRI reaches values comparable to the standard case, λ_MRI, being of the order of 100 m or smaller at small radii R ≲ 12km. Formally, there are extended potential MRI-unstable regions, at larger R and z, but their λ_eMRI is larger than in the main pseudo-torus. A wavelength that is too large may preclude MRI, as the linear analysis is valid only when λ_eMRI is much smaller than the characteristic length-scales of velocity and magnetic field variations.

Fig. 5. Comparison of the wavelength of the fastest-growing mode in the standard vs generalized case, at t≃ 40 ms after merger. Top: Standard MRI case, all coloured points are unstable to the MRI. Bottom: Generalized case. We mask the points where the instability condition is not met. Where the generalized instability is active, we observe λeMRI to reach values close to but slightly smaller than the standard case. Most interestingly, λmax is of the order of 100 m or smaller at small radii, r ≲ 12 km.

Finally, we could compare the values of τ_eMRI and τ_MRI, which are displayed in Fig. 6. The τ_MRI obtained from the standard MRI is a fraction of a millisecond for most of the region R ≲ 30km, which makes the instability sufficiently fast to have an impact on the post-merger dynamics. The one calculated with the generalized MRI analysis is displayed in the lower panel, where we mask the points in which the instability condition is not satisfied. When the generalized instability is active, we observe that τ_eMRI is close to or slightly larger than 1ms, particularly at large radii R ≳ 15 km. This is important, because such timescales are comparable to the dynamical timescales of the remnant. In other words, the basic assumption that the velocity and the large-scale (background) field change much more slowly than the instability growth time is not guaranteed, which might further hamper the possibility of developing MRI.

Fig. 6. Comparison of the timescale of the fastest-growing mode in the standard vs generalized case, at t≃ 40 ms after merger. Top: Standard MRI case. All coloured points are unstable to the MRI, and we observe the τMRI is a fraction of a millisecond at almost all radii ≲30 km, thus making the instability sufficiently fast to impact the post-merger dynamics. Bottom: Updated case. We mask the points where the instability condition is not met. Where the generalized instability is active, we observe that τeMRI is close to or slightly larger than 1ms, particularly at large radii, r ≳ 15km. This means that the instability takes longer to develop according to the generalized criterion, and might be too slow to significantly impact the post-merger dynamics.

4.3. Evolution of the MRI quantities

So far, we have illustrated the analysis at a given, representative time. In order to see the evolution, we show how the instability window, wavelength, and growth timescale evolve in time, in Fig. 7. We see that at t ≃ 20 ms the instability is possible only at small radii 7 ≲ R ≲ 12 km, and close to the equatorial plane. As time evolves (e.g. at t ≳ 50 ms after merger), the region extends to larger radii, but still close to the orbital plane, becomes more and more suited for the instability to develop. Away from these preferred regions, there are other locations in the R − z plane where MRI could potentially be active, although these appear to be somewhat spatially disconnected. For a region to be susceptible to MRI activity, it must be sufficiently extended to accommodate multiple wavelengths of the unstable modes.

Fig. 7. Evolution of the generalized 2D instability window. All coloured points are MRI-unstable according to the standard MRI criterion, while only the blue ones are according to the updated criterion. We see that at 20ms the instability is potentially active only at small radii, 7 ≲ R ≲ 12km, and close to the equatorial plane. As time evolves, we observe that the region at larger radii, but still close to the equatorial plane, also becomes more and more suited for the instability to develop, e.g. at t ≃ 50ms after merger.

The evolution of the modified τ_eMRI is displayed in Fig. 8. Focusing on the region close to the equatorial plane – the one that becomes progressively more suited for the instability to develop – we observe that the timescale of the instability decreases steadily, and is of the order of a fraction of a millisecond at late times. This means that the central region becomes more prone to the instability, but also that the instability would grow on shorter timescales. This suggests that the MRI only acquires an important role at late times t ≃ 100ms.

Fig. 8. Evolution of the timescale of the fastest-growing mode, τeMRI, of the extended MRI. All coloured points meet the generalized instability criteria.

The evolution of λ_eMRI is plotted in Fig. 9. Focusing again in the region close to the equatorial plane, we can observe that λ_eMRI decreases progressively, reaching values of the order of 1 − 10m at late times. Consistently with what we just said about the decreasing growth time, this is promising since λ_eMRI is much smaller than the magnetic and velocity field variation length-scales. However, numerically, it poses hard challenges to resolve it, due to an increasingly demanding spatial resolution over a larger region.

Fig. 9. Evolution of the wavelength of the fastest-growing mode, λeMRI, of the extended MRI. All coloured points meet the generalized instability criteria. Focusing on the region close to the equatorial plane – the one that becomes progressively more suited for the instability to develop – we see that λeMRI becomes progressively smaller and smaller, of the order of 1 − 10 m at late times. This means that in the central region the characteristic wavelength of the instability is small, thus making the instability more difficult to resolve.

All these qualitative results can be summarized quantitatively in Fig. 10, where the following averaged ratios over the 2D maps are represented as a function of time: the area in which the extended MRI is active (𝒜_eMRI) over the area in which the standard MRI is active (𝒜_MRI); the ratio λ_MRI/λ_eMRI; and the ratio τ_MRI/τ_eMRI. As we discussed previously, the instability region with the generalized criterion is smaller than with the standard one. The ratio λ_MRI/λ_eMRI is close to one because the values are similar when both instability criteria (standard and extended) are satisfied. As we noticed earlier though, λ_MRI takes on the largest values only when the standard MRI is active. This means that, effectively, it is harder to capture numerically the MRI according to the extended criterion. Even most relevant is the result of the ratios τ_MRI/τ_eMRI, which shows that the growth rate is significantly shorter when magnetic field gradients are considered in the analysis. One can include also the ratio of the areas to estimate an effective ratios τ_MRI/τ_eMRI, which is even smaller than the previous, especially at early times. These results suggest that the magnetic field gradients might be suppressing (or at least slowing down) the axisymmetric MRI.

Fig. 10. Averaged ratios of relevant quantities as a function of time. We show the time evolution of three ratios: the area in which the extended MRI is active (𝒜eMRI) over the area in which the standard MRI is active (𝒜MRI); the ratio λMRI/λeMRI; and the ratio τMRI/τeMRI. Interestingly, the low value of τMRI/τeMRI (especially when multiplied by 𝒜MRI/𝒜eMRI to estimate the effective ratio) indicates that the magnetic field gradients might be diminishing the MRI.

5. Outlook and conclusions

The MRI is considered a key mechanism to induce turbulence in accretion disks and is believed to be a key process mediating the emergence of large-scale, poloidal magnetic fields in the post-merger remnants of compact binary mergers (Duez et al. 2006; Siegel et al. 2013; Shibata 2015; Ruiz et al. 2016; Ciolfi et al. 2019; Ciolfi 2020; Radice & Hawke 2024; Abac et al. 2025). These simulations show a complex magnetic field topology, which differs significantly from the simplified conditions under which MRI was originally formulated. In particular, certain simulations of long-lived remnants fail to provide evidence of poloidal magnetic field amplification after the vigorous Kelvin-Helmholtz amplification phase–while clearly showing the impact of linear winding in the toroidal component (Palenzuela et al. 2022; Aguilera-Miret et al. 2023, 2024).

In this study we tackled this question by investigating, for the first time, how complex magnetic field topologies, specifically magnetic field gradients, influence the axisymmetric MRI. The strategy adopted here is parallelled by the analysis of Celora et al. (2023), who show that the conventional MRI criteria and results crucially depend on the assumption of an axisymmetric background. After having explored the impact of several magnetic field background configurations, we find that radial magnetic field gradients can have a notable impact on the MRI: the ‘instability window’ is reduced and both the wavelength and growth rate of the MRI’s fastest-growing mode are altered.

Next, we examined the generalized conditions for the MRI using two simplified models inspired from the accretion and protoplanetary disks described in the literature (Pelletier & Pudritz 1992; Galli et al. 2006; Jafari 2019). This exercise serves a dual purpose: firstly, it demonstrates that while our initial focus is on binary neutron star applications, our findings have broader implications. Secondly, it helps in developing an intuition for how additional magnetic gradients influence the instability.

Finally, we investigated the generalized MRI conditions in the post-merger environment by postprocessing a simulation of a long-lived binary neutron star merger remnant. In particular, the simulation considered effectively captures magnetic field growth from an initial realistic field strength of approximately 10¹¹ G (Palenzuela et al. 2018; Viganò et al. 2019b). Our findings indicate that the (generalized) MRI can be active primarily in regions of small radii (5 ≲ R ≲ 15 km) during the early stages (t ≲ 40ms). As time progresses, the region close to the equatorial plane (−10 ≲ z ≲ 10 km) and large radii (R ≳ 15 km) also become susceptible to MRI development. However, the e-folding timescale of the instability remains consistently close to or slightly above 1ms across all radii and up to t ≳ 50 ms after mergers. This timescale decreases to a fraction of a millisecond (the typical value according to conventional MRI) only at late times, t ≳ 100 ms. Additionally, we observed that the wavelength of the fastest-growing MRI mode is generally shorter compared to the standard case: of the order of 10 − 100 m at most in regions where the MRI would potentially be active, namely small radii at early times, or all radii but close to the equatorial plane at late times. Notably, the value at small radii is consistently around 1 − 10 m throughout the simulation, particularly at late times. Such small values are, on the one hand, promising for the applicability of linear analysis, as a key requirement for its validity is that the dominant unstable MRI modes have wavelengths much shorter than the characteristic scales over which the magnetic and velocity fields vary. On the other hand, accurately resolving such short wavelengths presents significant challenges for current computational resources–see also the discussion in Appendix C, where we report on the empirical quality criteria commonly used to quantify the ability of a simulation to capture MRI-driven MHD turbulence. A final caveat we feel compelled to raise concerns the magnetic non-axisymmetric indicators discussed in Appendix B, which remain close to unity across much of the domain even at late times, 𝒪(100) ms. While we caution that non-axisymmetric indicators being close to unity does not necessarily imply the absence of growth of an axisymmetric large-scale poloidal field, this nevertheless further challenges the expected influence of the axisymmetric MRI. It also suggests that the characteristic timescales and length-scales of magnetic field variations are significantly shorter than commonly assumed. This caveat could also apply to certain protoplanetary disks where no organized large-scale magnetic field is observed (Huang et al. 2024, 2025).

In summary, the results presented in this paper offer a coherent explanation for the lack of growth in the poloidal component (after the Kelvin-Helmholtz unstable phase) in the simulation examined here. Most notably, the longer timescales associated with the generalized instability, combined with its later onset, would suggest that the axisymmetric MRI’s influence on post-merger dynamics is significantly reduced. In turn, this indicates that the role of the MRI in affecting the post-merger dynamics should be approached with care in modelling efforts. Finally, we caution that setting an unrealistically high initial large-scale magnetic field, instead of starting with a weak field amplified by KHI and winding, biases the post-merger configuration, and the magnetic gradients in particular (Aguilera-Miret et al. 2024; Gutiérrez et al. 2025).

Acknowledgments

TC is an ICE Fellow and is supported through the Spanish program Unidad de Excelencia Maria de Maeztu CEX2020-001058-M. CP’s work was supported by the Grant PID2022-138963NB-I00 funded by the Spanish Ministry of Science and Innovation, MCIN/AEI/10.13039/501100011033/FEDER, UE. DV is supported through the Spanish program Unidad de Excelencia Maria de Maeztu CEX2020-001058-M as well as by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Starting Grant – IMAGINE- No. 948582, PI: DV). RA-M is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under the Germany Excellence Strategy - EXC 2121 ‘Quantum Universe’ – 390833306. The authors thankfully acknowledge RES resources provided by BSC in MareNostrum, projects AECT-2024-1-0010, AECT-2024-2-0004 and AECT-2024-3-0007.

References

Appendix A: Impact of more general gradients on the MRI

As we have seen that radial gradients in the magnetic field can have a crucial impact on the MRI, in this appendix we consider the impact of additional gradients on the instability. Specifically, we consider the case where we included vertical gradients in the background magnetic field as well as in the velocity (on top of differential rotation). In essence this means we took

$$\begin{array}{cc}& \mathbf{v}=[\mathrm{\Omega }(R,z)-\mathrm{\Omega }(R_0,z_0)](R_0+x)\widehat{y},\hfill \end{array}$$(A.1a)

$$\begin{array}{cc}\hfill & \mathbf{B}=B^y(z)\widehat{y}+B^z,{\partial }_i{B}^i={\partial }_z{B}^z=0,\hfill \end{array}$$(A.1b)

so that the gradients term entering the WKB expansions are

$$\begin{array}{cc}& {\partial }_i{v}_j={\delta }_{j2}{\delta }_{i1}{S}_0+{\delta }_{j2}{\delta }_{i3}{S}_1,S_1=\frac{R}{z}\left(\frac{\partial \mathrm{\Omega }}{\partial logz}\right),\hfill \end{array}$$(A.2a)

$$\begin{array}{cc}\hfill & {\partial }_i{B}_j={\delta }_{i3}{\delta }_{j2}({\partial }_z{B}^y),\hfill \end{array}$$(A.2b)

while S₀ is defined as in Eq. (3). We also used the no-monopoles constraint to kill gradients in the vertical magnetic field component.

Performing the linear analysis under the same simplifying assumptions as in Sect. 2.2, the perturbation equations corresponding to the background configuration under consideration read

$$\begin{array}{cc}& k_x\delta v^x+k_z\delta v^z=0,\hfill \end{array}$$(A.3a)

$$\begin{array}{cc}\hfill -i& \omega \delta v^x-2{\mathrm{\Omega }}_0\delta v^y+i{k}_x(v_A^y\delta v_A^y+v_A^z\delta v_A^z)-i\delta v_A^x{k}_z{v}_A^z=0,\hfill \end{array}$$(A.3b)

$$\begin{array}{cc}\hfill -i& \omega \delta v^y+\frac{{\kappa }^2}{2\mathrm{\Omega }}\delta v^x-i{k}_z{v}_A^z\delta v_A^y+S_1\delta v^z-a\delta v_A^z=0,\hfill \end{array}$$(A.3c)

$$\begin{array}{cc}\hfill -i& \omega \delta v^z+i{k}_z{v}_A^y\delta v_A^y+a\delta v_A^y=0,\hfill \end{array}$$(A.3d)

$$\begin{array}{cc}\hfill -i& \omega \delta v_A^x-i{v}_A^z{k}_z\delta v^x=0,\hfill \end{array}$$(A.3e)

$$\begin{array}{cc}\hfill -i& \omega \delta v_A^y-i{v}_A^z{k}_z\delta v^y-S_0\delta v_A^x+a\delta v_A^z-S_1\delta v_A^z=0,\hfill \end{array}$$(A.3f)

$$\begin{array}{cc}\hfill -i& \omega \delta v_A^z-i{v}_A^z{k}_z\delta v^z=0,\hfill \end{array}$$(A.3g)

where we have introduced a = ∂_zv_A^y to slim down the notation. As before, we first rewrote the various perturbed quantities in terms of δv^z:

$$\begin{array}{cc}& \delta v^x=-\frac{k_z}{k_x}\delta v^z,\hfill \end{array}$$(A.4a)

$$\begin{array}{cc}\hfill & \delta v^y=i\left(\frac{\omega }{{\omega }^2-p^2}\right)\frac{k_z}{k_x}[\frac{{\kappa }^2}{2\mathrm{\Omega }}-S_0\frac{p^2}{{\omega }^2}-S_1\frac{k_x}{k_z}(1-\frac{p^2}{{\omega }^2})]\delta v^z,\hfill \end{array}$$(A.4b)

$$\begin{array}{cc}\hfill & \delta v_A^x=\frac{k_z}{k_x}\frac{p}{\omega }\delta v^z,\hfill \end{array}$$(A.4c)

$$\begin{array}{cc}\hfill & \delta v_A^y=\frac{\omega }{k_z{v}_A^y-ia}\delta v^z,\hfill \end{array}$$(A.4d)

$$\begin{array}{cc}\hfill & \delta v_A^z=-\frac{p}{\omega }\delta v^z,\hfill \end{array}$$(A.4e)

and substituting these into the x-component of the Euler equation to arrive at

$$\begin{array}{cc}\hfill (1+\frac{k_x^2}{k_z}\frac{v_A^y}{v_A^y-ia})& {\omega }^4-[2p^2+p^2\frac{k_x^2}{k_z^2}+{\kappa }^2-2\mathrm{\Omega }{S}_1\frac{k_x}{\mathit{kz}}+\frac{k_x^2}{k_z}\frac{v_A^y{p}^2}{v_A^y{k}_z-ia}]{\omega }^2\hfill \\ \hfill & +p^2[2\mathrm{\Omega }{S}_0-2\mathrm{\Omega }{S}_1\frac{k_x}{k_z}+p^2(1+\frac{k_x^2}{k_z^2})]=0.\hfill \end{array}$$(A.5)

It is easy at this point to see that if we consider vertical modes only (k_x = 0, as in Sect. 2), the dispersion relation reduces to the standard MRI case we discussed in Sect. 2.3. This means that adding vertical gradients in the background velocity and/or magnetic field, separately or combined, does not impact on the usual MRI phenomenology and criteria.

Without showing the details of the calculation, let us also mention that one can make progress in the case where we included radial and vertical gradients in the velocity, and also radial gradients in the magnetic field. In this last case in fact the dispersion relation can be reduced to the form

$$\begin{array}{c}\hfill {\omega }^4-(2p^2+{\kappa }^2+f){\omega }^2+[p^4+2\mathrm{\Omega }{S}_0{p}^2+2i\mathrm{\Omega }{S}_{1}p({\partial }_x{v}_A^z)+f{p}^2]=0,\end{array}$$(A.6)

while the function f is defined as in Eq. (16). While the equation is now complex, thus preventing a simple analytical solution, it is easy to see that when considering vertical gradients in the velocity as well as radial gradients in both the velocity and the magnetic field, the solution to the approximated dispersion relation is identical to the case with radial gradients only.

All the different cases analysed here and in Sect. 2.2 are conveniently summarized in Table 1, which makes it evident the central role played by radial magnetic field gradients in the generalized instability.

Appendix B: A closer look at the topological indicators

The aim of this appendix is to discuss more in detail the indications we have from the simulations presented in Aguilera-Miret et al. (2023, 2024). We focus in particular on the two topological indicators that are introduced there, and can offer useful insights into the amount of turbulence in the post-merger environment. By means of a Fourier-type decomposition, we argue that the behaviours observed when these indicators are evaluated on the simulations output is supportive of the generalized MRI conditions we studied in this work.

Fig. B.1. Evolution of the total magnetic non-axisymmetric indicators, Imtot. We observe the non-axisymmetric indicator is Im ∼ 1 even at late times (t ≳ 150) ms, particularly at large radii and close to the equatorial plane. We focus here on the total indicator, which is dominated by the azimuthal contribution. When plotting the poloidal indicators, we observe that they are ≈1 at all radii and vertical distances from the equatorial plane even at late times, t ≃ 150 ms.

The non-axisymmetry indicators are introduced as follows:

$$\begin{array}{c}\hfill I_K^i=\frac{\delta E^i}{{\overline{E}}_K^i+\delta E^i},I_m=\frac{\delta E_m^i}{{\overline{E}}_m^i+\delta E_m^i},\end{array}$$(B.1)

where

$$\begin{array}{cc}\hfill {\overline{E}}_K^i& ={\overline{U}}_K^i{\overline{U}}_i^K,\delta E_K^i=\delta U_K^i\delta U_i^K,\hfill \end{array}$$(B.2a)

$$\begin{array}{cc}\hfill {\overline{E}}_m^i& ={\overline{U}}_m^i{\overline{U}}_i^m,\delta E_m^i=\delta U_m^i\delta U_i^m,\hfill \end{array}$$(B.2b)

and the two vector fields considered are $U_K^i=\sqrt{\rho }{v}^i$ and U_mⁱ = Bⁱ. These vector fields (components) are decomposed into an average along the azimuthal direction (cf. Eq. (34) in main text) and an integrated residual,

$$\begin{array}{cc}\hfill \delta U^i(R,z)& =\frac{{\int }_l|U^i-{\overline{U}}^i|{\chi }^{-1/2}Rd\phi }{{\int }_l{\chi }^{-1/2}Rd\phi }.\hfill \end{array}$$(B.3)

Evaluating these indicators on the simulations output—considering either the dominant azimuthal component only or the summed contributions of the individual components—the authors find the kinetic axisymmetric indicator to be I_K ≃ 0.05 at most—thus showing a high degree of axisymmetry—while the magnetic one I_m ∼ 1. This is true both at early (t ≳ 20ms) and late times (t ≳ 150ms) and at all radii R ≳ 10km. At the same time though, the residual energies are similar in values δE_K^φ ≃ δE_m^φ at all radii, and until late times t ∼ 150 ms.

To unpack the information encoded in the indicators, we introduced a Fourier decomposition in the azimuthal direction,

$$\begin{array}{c}\hfill U^{\phi }=c_0(R,Z)+{\displaystyle \sum _{n=1}^{\infty }}{a}_n(R,z)cos(n\phi )+{\displaystyle \sum _{n=1}^{\infty }}{b}_n(R,z)sin(n\phi ),\end{array}$$(B.4)

and rewrote the averaged energy as ${\overline{E}}^{\phi }={(c_0)}^2$. As for the residual, we can write the following inequality,

$$\begin{array}{c}\hfill \delta U^{\phi }\le {\displaystyle \sum _{n=1}^{\infty }}\frac{2}{\pi }(|a_n|+|b_n|),\end{array}$$(B.5)

so that for the purpose of the following argument we can write

$$\begin{array}{c}\hfill \delta E^{\phi }\simeq {\displaystyle {\left(\sum |a_n|+|b_n|\right)}^2},I\simeq {[1+\frac{{(c_0)}^2}{{{\displaystyle \left(\sum |a_n|+|b_n|\right)}}^2}]}^{-1}.\end{array}$$(B.6)

In terms of this decomposition we can make the following observations

1. As I_K ≪ 1 the fluid motion is mostly axisymmetric, and for the purpose of a linear analysis the background velocity at t ≳ 20ms can be obtained using azimuthal averages.

2. Since I_m ≃ 1 but δE_K ≃ δE_m, we see that the difference in the axisymmetric indicators really comes from the fact that the averaged kinetic energy is dominating over the magnetic one—the background is in a ‘kinematic regime’ in a loose sense.

3. Since I_m ≃ 1, δE_K ≃ δE_m at all radii and up until late times (cf. Fig. 5b of Aguilera-Miret et al. 2024), the azimuthal average of the magnetic field varies (in the radial direction) on the same scales as the residuals in the velocity.

In short then, the indicators suggest that in order to analyse the role of the MRI in the simulations, we first need to generalize the standard criteria to the case where we considered (radial) gradients in the magnetic field on top of differential rotation, as discussed in Sect. 2.

In this work we often considered 2D maps of the relevant averaged quantities, showing them as a function of the cylindrical coordinates (R, z). As such, we here show in Fig. B.1 the evolution of the total magnetic non-axisymmetric indicator at some representative times—which is dominated by the azimuthal component. The figure clearly shows that the non-axisymmetric indicator is close to 1 in large parts of the domain even at late times 𝒪(100) ms. Effectively, this raises an additional caveat on the applicability of the MRI criteria in the simulation considered here. We also note a similar point may apply to some protoplanetary disks where no ordered large-scale field is observed (Huang et al. 2024, 2025).

Appendix C: Empirical quality factors

In this appendix we present several quantitative diagnostics that have been established as useful empirical criteria for assessing the adequacy of a simulation in capturing MRI-driven MHD turbulence (Hawley et al. 2011, 2013). Specifically, we focus on the vertical and azimuthal quality factors Q_z and Q_φ, defined as the ratio of the corresponding characteristic MRI wavelength (λ_MRI) to the grid spacing in the vertical and azimuthal directions, respectively. Following Kiuchi et al. (2018), we also examined the ratio of Maxwell stress to magnetic pressure α_mag, as well as the ratio between the radial and toroidal magnetic components $\mathcal{R}=B_R^2/B_{\phi }^2$. These quantities are plotted as functions of radius at some representative snapshots after the merger.

As shown in Fig. C.1, the vertical quality factor remains around Q_z ∼ 8–10 at radii larger than 10–15km, while the azimuthal quality factor reaches Q_φ ≳ 30. Although the vertical quality factor suggests the resolution is only marginally sufficient, the higher values of Q_φ are expected to compensate for the modest Q_z, consistent with previous convergence studies. We also show the radial profile of α_mag, finding that it saturates at approximately 0.25-0.3 shortly after the merger and remains roughly constant thereafter. Overall, the computed global quantities are consistent with –though slightly lower than– those reported in the literature when simulations are performed at sufficiently high resolution (see e.g. Kiuchi et al. 2018). In particular, we note that the two main differences between our simulations and those in the literature are: (i) we started with a significantly weaker (and more realistic) initial magnetic field, and (ii) we employed LES techniques with the gradient sub-grid-scale model to capture more accurately the magnetic field amplification during the turbulent phase.

Fig. C.1. Quantities of empirical criteria.Top: Fastest-growing MRI vertical wavelength (λMRIz), its corresponding quality factor (Qz), and the ratio $\mathcal{R}=B_R^2/B_{\phi }^2$. Bottom: Fastest-growing MRI toroidal wavelength (λMRIφ), toroidal quality factor (Qφ), and the ratio of Maxwell stress to magnetic pressure (αmag). All these quantities are shown as a function of distance (in km) for different times after the merger.

We also stress that in this appendix the quality factors Q_z and Q_φ as a function of the radial distance R have been computed using the standard MRI wavelength λ_MRI to enable a clearer comparison with results reported in the literature. It is reasonable to expect that generalized quality criteria based on the extended MRI wavelengths will also provide useful information on the adequacy of a simulation, since λ_eMRI reduces to λ_MRI in the absence of significant gradients and these criteria are derived from linear analysis. The situation is different for the remaining diagnostics considered, α_mag and ℛ, which are directly linked to the non-linear development of MHD turbulence (Hawley et al. 2011). As the impact on the characteristic wavelength is moderate compared to that on the timescale, we expect the quality factors to be moderately smaller, which may help explain the slightly lower values observed for the computed global quantities. Finally, we note that smaller-scale magnetic fields, such as those generated during a Kelvin–Helmholtz amplification phase, are more susceptible to numerical dissipation, which is also known to impact on the efficiency of MRI-driven dynamos in local simulations (e.g. Fromang & Papaloizou 2007).

All Tables

All Figures

Fig. 1. Extended MRI analysis for the analytical model given by a vertical magnetic field with radial power-law decay in Eq. (25). Left: Generalized instability criterion. Middle: Comparison between the generalized and standard MRI for the wavelength of the fastest-growing mode. Right: Timescale of the instability. We have considered the values with n = 2,  4 and C = 0.01,  0.025,  0.1 in the model defined by Eq. (25). The wavelength and timescale are shown only when the instability is active according to the generalized criterion, a transition that is marked with black dots.

In the text

Fig. 2. Extended MRI analysis for the analytical model given by a vertical and toroidal magnetic field. From left to right: Analytic instability window, ratio of the wavelength of the fastest-growing mode, and ratio of the timescale of the instability in the generalized case vs the standard MRI result. The wavelength and timescale of the fastest-growing mode are shown only where the instability is active according to the generalized criterion.

In the text

	Fig. 3. 2D maps of the azimuthally averaged magnetic field magnitudes at t≃ 40 ms after merger. From left to right: 2D maps of the magnitude of the radial, azimuthal, and vertical component of the magnetic field. We see the different components present large gradients, which we take into account in the generalized MRI criterion.
In the text

Fig. 4. Instability window calculation at t≃ 40 ms after merger. From the angular frequency (Ω) we computed its gradient (s0) and mask those points where s0 > 0, namely those for which the standard MRI would be inactive. Then we evaluated the updated instability criterion, considering at once the two cases f > 0 and f < 0. In the right panel, blue points are unstable according to the updated MRI criterion, and red are those for which the magnetic gradients are large enough to switch off the instability; all coloured points are MRI-unstable according to the standard criterion.

In the text

Fig. 5. Comparison of the wavelength of the fastest-growing mode in the standard vs generalized case, at t≃ 40 ms after merger. Top: Standard MRI case, all coloured points are unstable to the MRI. Bottom: Generalized case. We mask the points where the instability condition is not met. Where the generalized instability is active, we observe λeMRI to reach values close to but slightly smaller than the standard case. Most interestingly, λmax is of the order of 100 m or smaller at small radii, r ≲ 12 km.

In the text

Fig. 6. Comparison of the timescale of the fastest-growing mode in the standard vs generalized case, at t≃ 40 ms after merger. Top: Standard MRI case. All coloured points are unstable to the MRI, and we observe the τMRI is a fraction of a millisecond at almost all radii ≲30 km, thus making the instability sufficiently fast to impact the post-merger dynamics. Bottom: Updated case. We mask the points where the instability condition is not met. Where the generalized instability is active, we observe that τeMRI is close to or slightly larger than 1ms, particularly at large radii, r ≳ 15km. This means that the instability takes longer to develop according to the generalized criterion, and might be too slow to significantly impact the post-merger dynamics.

In the text

Fig. 7. Evolution of the generalized 2D instability window. All coloured points are MRI-unstable according to the standard MRI criterion, while only the blue ones are according to the updated criterion. We see that at 20ms the instability is potentially active only at small radii, 7 ≲ R ≲ 12km, and close to the equatorial plane. As time evolves, we observe that the region at larger radii, but still close to the equatorial plane, also becomes more and more suited for the instability to develop, e.g. at t ≃ 50ms after merger.

In the text

	Fig. 8. Evolution of the timescale of the fastest-growing mode, τeMRI, of the extended MRI. All coloured points meet the generalized instability criteria.
In the text

Fig. 9. Evolution of the wavelength of the fastest-growing mode, λeMRI, of the extended MRI. All coloured points meet the generalized instability criteria. Focusing on the region close to the equatorial plane – the one that becomes progressively more suited for the instability to develop – we see that λeMRI becomes progressively smaller and smaller, of the order of 1 − 10 m at late times. This means that in the central region the characteristic wavelength of the instability is small, thus making the instability more difficult to resolve.

In the text

Fig. 10. Averaged ratios of relevant quantities as a function of time. We show the time evolution of three ratios: the area in which the extended MRI is active (𝒜eMRI) over the area in which the standard MRI is active (𝒜MRI); the ratio λMRI/λeMRI; and the ratio τMRI/τeMRI. Interestingly, the low value of τMRI/τeMRI (especially when multiplied by 𝒜MRI/𝒜eMRI to estimate the effective ratio) indicates that the magnetic field gradients might be diminishing the MRI.

In the text

Fig. B.1. Evolution of the total magnetic non-axisymmetric indicators, Imtot. We observe the non-axisymmetric indicator is Im ∼ 1 even at late times (t ≳ 150) ms, particularly at large radii and close to the equatorial plane. We focus here on the total indicator, which is dominated by the azimuthal contribution. When plotting the poloidal indicators, we observe that they are ≈1 at all radii and vertical distances from the equatorial plane even at late times, t ≃ 150 ms.

In the text

Fig. C.1. Quantities of empirical criteria.Top: Fastest-growing MRI vertical wavelength (λMRIz), its corresponding quality factor (Qz), and the ratio $\mathcal{R}=B_R^2/B_{\phi }^2$. Bottom: Fastest-growing MRI toroidal wavelength (λMRIφ), toroidal quality factor (Qφ), and the ratio of Maxwell stress to magnetic pressure (αmag). All these quantities are shown as a function of distance (in km) for different times after the merger.

In the text