Pith. sign in

REVIEW 3 major objections 4 minor 18 references

Numerical evolution of the resistive relativistic magnetohydrodynamic equations: a minimally implicit Runge-Kutta scheme

T0 review · 3 major / 4 minor · reviewed 2026-08-10 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read This paper argues that the stiff Ohmic source in resistive relativistic MHD can be advanced by a minimally implicit Runge-Kutta update whose implicit step is a single algebraic division, giving explicit-like cost and no added recovery…

desk verdict A promising, clearly explained RRMHD stiff-source scheme whose headline stability claim rests on an omitted linear-stability derivation; worth refereeing, but the missing analysis and the unexplained order reduction need to be addressed. read the letter →

arxiv 2502.00990 v1 pith:WBIL6HFI submitted 2025-01-19 physics.comp-ph astro-ph.HE

Isabel Cordero-Carrión , Samuel Santos-Pérez , Clara Martínez-Vidallach This is my paper · ORCID
classification physics.comp-phastro-ph.HE MSC 35B3535F5035L6035Q3535Q7565M0665M1276W05
keywords ResistiverelativisticmagnetohydrodynamicsMinimallyimplicitRunge-Kutta(MIRK)StiffsourcetermsImplicit-explicit(IMEX)methodsPrimitivevariablerecoveryFinite-differenceSelf-similarcurrentsheetCircularlypolarizedAlfvénwaves
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper proposes Minimally Implicit Runge-Kutta (MIRK) schemes that make the stiff Ohmic source term in resistive relativistic magnetohydrodynamics (RRMHD) cost about the same as an explicit update. The key move is to evaluate implicitly only the electric- and magnetic-field components multiplied by the conductivity, which are conserved variables themselves, so the implicit step is a single algebraic division instead of an iterative primitive-variable recovery. First- and second-order versions are constructed, with coefficients chosen from wave-propagation PIRK methods and a linear stability analysis in the ideal infinite-conductivity limit. The schemes are tested on a self-similar current sheet and on a circularly polarized Alfvén wave, and the authors argue that existing explicit codes can be adapted to them with minimal changes. The central claim, read sympathetically, is that physical resistivity can be added to relativistic MHD simulations without the extra cost or convergence risk of implicit-explicit (IMEX) recovery loops.

What carries the argument

The central object is the MIRK (Minimally Implicit Runge-Kutta) update for the electric field: in the evolution equation $\partial_t E^j = S_E^j - \bar\sigma [E^j + (v\times B)^j - (v_l E^l)v^j]$, only the bracketed $E$- and $B$-components are evaluated at the new stage, while $\bar\sigma = \sigma W$ and $v$ stay at their explicit stage values. For the first-order scheme the update becomes $E^j|_{n+1} = (E^j|_n + \Delta t [S_E^j|_n + \bar\sigma|_n E^l|_n (v^j|_n v_l|_n - \delta^j_l) - \bar\sigma|_n (v|_n \times B|_{n+1})^j])/(1 + \Delta t\bar\sigma|_n)$, an explicit-looking step with effective time step $\Delta t/(1+\Delta t\bar\sigma)$; the second-order scheme applies the same idea in two stages with $\bar\sigma$ evaluated at the old and first-stage values. The remaining coefficients are fixed by PIRK coefficients for wave-like behavior and by a linear stability analysis in the $\sigma\to\infty$ limit: the ideal limit requires a zero eigenvalue of multiplicity at least three, so $E\to -v\times B$, and the other eigenvalue must stay inside the unit disk for any $v$, which selects $c_1 = -1/\sqrt{2}$ for the second-order method.

What would settle it

Assemble the amplification matrix of the MIRK update for the full RRMHD system with the stated coefficients at $\sigma\to\infty$ or at a large finite $\sigma$, for generic velocity directions, and check whether every eigenvalue lies in the unit disk for $\Delta t$ up to the CFL limit; any admissible velocity and step size with spectral radius above 1 would directly contradict the claimed stability, since the paper does not display this calculation.

Watch

Extended reading notes

Core claim

The paper's central claim is that the stiffness of Ohm's law in resistive relativistic MHD can be handled by treating the conductivity-weighted electric and magnetic field components implicitly while keeping every other term explicit; because those components are already conserved variables, the implicit stage reduces to solving a diagonal linear system, i.e. a division. The authors present a first-order scheme and a second-order two-stage scheme, with coefficients selected by PIRK wave-propagation coefficients and by a linear stability analysis in the ideal $\sigma\to\infty$ limit that fixes the remaining parameters, including $c_1=-1/\sqrt{2}$ for the second-order method. They validate the schemes on a self-similar current sheet at $\sigma=10^3$ and on a circularly polarized Alfvén wave at $\sigma=10^8$, recovering the exact and ideal-MHD solutions, and they argue that the computational cost stays comparable to explicit methods with no additional primitive-variable recovery loops.

Load-bearing premise

The scheme's stability rests on a linear stability analysis in the perfect-conductor limit that the paper asserts without displaying, and on that analysis continuing to hold at finite conductivity and for the full set of coupled fluid and field equations.

Editorial extensions

If this is right

  • Existing explicit RRMHD codes can incorporate the stiff Ohm term by replacing the electric-field update with the MIRK formula, a division by $1+\Delta t\bar\sigma$, without changing the rest of the time-stepping structure.
  • The number of primitive-variable recoveries per timestep stays the same as in an explicit scheme, which removes the main extra cost and the convergence risk that IMEX methods introduce.
  • At high conductivity the electric field is driven toward the ideal value $-v\times B$ between stages, so the ideal MHD limit is recovered as $\bar\sigma\to\infty$.
  • The same minimally implicit construction, making implicit only the terms $G_i(U)U^i$ that multiply a large parameter while evaluating $G_i$ explicitly, extends to other hyperbolic balance laws with stiff sources.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Extending beyond the paper: because the implicit solve is a single diagonal division, the same MIRK update should transfer almost unchanged to multidimensional and shock-dominated RRMHD, but the decisive test would be a 2D or 3D problem with non-smooth data, which the paper explicitly leaves to future work.
  • The paper does not use the term asymptotic-preserving, but its effective time step $\Delta t/(1+\Delta t\bar\sigma)$ formally drives $E$ toward $-v\times B$ as $\bar\sigma\to\infty$; this suggests the method could serve as a building block for schemes that switch between resistive and ideal regimes, provided the explicitly treated source terms are also controlled.
  • The observed instability of the circularly polarized Alfvén wave test at CFL=0.8, attributed by the authors to the explicitly evolved source terms, implies that the practical barrier to large timesteps may sit in $S_E$, $S_B$, and $S_Y$ rather than in the conductivity term; a natural extension would be to apply the same minimal-implicitness idea to those terms, which the paper does not pursue.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The manuscript introduces first- and second-order minimally implicit Runge-Kutta (MIRK) schemes for the resistive relativistic magnetohydrodynamic (RRMHD) equations, following Komissarov's augmented-system formulation. In these schemes only the electric-field components appear implicitly, and the resulting algebraic solve is analytic, so the method is claimed to retain a computational cost comparable to explicit schemes while remaining stable at high conductivity. The paper derives the two schemes, motivates the coefficient choices through a linear stability analysis in the infinite-conductivity limit, and validates the methods on two one-dimensional smooth tests: the self-similar current sheet and circularly polarized Alfvén waves. The current-sheet test agrees well with the exact solution, and the Alfvén-wave test approaches the ideal MHD limit at sigma = 10^8, though the measured convergence order of the second-order method is first order in that test.

Significance. If the central claims hold, the MIRK approach is a useful contribution to numerical RRMHD: it avoids the repeated primitive-variable recoveries required by IMEX schemes, and the analytic inversion of the implicit operator makes it easy to retrofit into existing explicit codes. The paper is transparent about the algebraic form of both schemes, and the two tests are checked against independent exact or reference solutions, which is a genuine strength. The potential significance is, however, moderated by the fact that the method is validated only on smooth one-dimensional problems and that the stability assertion rests on an unshown analysis. The practical claim of 'stable numerical evolutions at explicit-like cost' is plausible but not fully established by the evidence presented.

major comments (3)
  1. [Section 4, after Eqs. (16) and (22)] The coefficient choices c3=1, c6=1/2, c4=(1-c1)^2/(2c1), and c1=-1/sqrt(2) are said to follow from a linear stability analysis of the evolution system in the infinite-conductivity limit, but neither the linearized system, the amplification matrix, nor the eigenvalue calculation is shown. This is load-bearing because the claimed stability of the schemes and the selection of all ci coefficients depend on that analysis. The authors should provide the full amplification matrix and the eigenvalue bounds, including the statement that the remaining eigenvalue is bounded by 1 in absolute value independently of the velocity field, and should state any assumptions on the discretization of the spatial derivatives used in the analysis.
  2. [Section 5.2, Table 2] Table 2 shows that the second-order MIRK method converges only with order approximately 1.0 in the CP Alfvén wave test (p values near 0.87-0.98 for the finer resolutions), while the method is presented as second-order. The text states that the authors were not able to find the reason for this reduction. This is a central issue: the method's advertised order is not realized in the coupled hydrodynamic test. The authors should either identify the source of the order reduction (e.g., the Kreiss-Oliger dissipation term, the CFL value, or the recovery procedure) or explicitly present the second-order method as first-order in practice for coupled evolutions.
  3. [Section 5.1 and Section 5.2] The paper claims in the abstract that the MIRK methods 'are able to deal with stiff terms producing stable numerical evolutions,' but the stability demonstration is limited to two smooth one-dimensional tests with CFL numbers at or below 0.8 (and 0.3 for the Alfvén wave test). Section 5.1 itself shows that at large CFL values oscillations appear, which are attributed to the explicitly treated source terms, and Section 6 concedes that non-smooth and more complex simulations are needed. The stability claim should therefore be qualified to the tested regime, and the conditions under which the method is stable (including restrictions on CFL and on the explicitly treated source terms) should be stated precisely.
minor comments (4)
  1. [Figures 6 and 7] The legends in Figures 6 and 7 label the reference solution as 'Exact solution v=0' even though the simulations are performed with vx=0.1; the legend should be corrected to avoid confusion about which solution is being compared.
  2. [Section 5.2, after Eq. (36)] The text says that for the convergence study the power of the factor Delta x in the Kreiss-Oliger term is increased by one unit and the coefficient is set to 0.1, but the exact modified form of the dissipation term is not written; please specify the expression used in the convergence runs.
  3. [Section 6] The statement that 'MIRK methods have the same computational cost as explicit ones' is stronger than what is demonstrated, since no wall-clock timing or operation-count comparison is reported and the implicit solve, although analytic, still introduces additional algebraic operations; the cost claim should be stated as 'the same order of cost per time step'.
  4. [Section 1 and Section 4] The text refers to 'terms in orange' in Eq. (9) and the corresponding explanation, but in a monochrome printout or in some rendering of the manuscript this color coding is not visible; please add an explicit mathematical designation of the implicitly treated terms.

Circularity Check

0 steps flagged · score 2.0 of 10

No circular reduction: the MIRK coefficients come from a linear-stability analysis (stated but not displayed) and the two numerical tests use independent exact/reference solutions; the only self-citations are prior PIRK method papers and are not load-bearing.

full rationale

The derivation chain is not circular. The first- and second-order MIRK coefficients are selected from boundedness conditions in the infinite-conductivity limit, the ideal-regime constraint E = -v x B, and a linear-stability analysis that is asserted but not shown in Section 4 (between Eqs. (16) and (22)); none of these conditions is fitted to the numerical solutions presented later. The two tests are independent benchmarks: the current-sheet test is compared to the analytic erf solution in Eq. (27), and the CP Alfven test is compared to the exact ideal-MHD limit solution specified in Eqs. (29)-(30), so the stiff-limit behavior is checked against an external reference rather than predicted from fitted data. The self-citations to the authors' prior PIRK papers [16,17] are used only to set the wave-like coefficients c2 and c5; they are prior published method constructions, and the remaining coefficients c1, c3, c4 and c6 are additionally fixed by the paper's own stability requirements, so these citations are not load-bearing. The paper's own admitted gaps - the omitted linear-stability calculation, the unexplained first-order convergence of the second-order method in Section 5.2, and the statement in Section 6 that non-smooth and more complex tests are still needed - are support and validation limitations, not circular reductions. The score is therefore a low 2, reflecting only the minor non-load-bearing self-citation and the unshown stability analysis, with no circular step identified.

Assumptions & free parameters 2 free parameters · 4 assumptions · 0 invented entities

The only hand-chosen numbers are the Runge-Kutta coefficients, selected by a stability analysis at infinite conductivity, and the artificial dissipation coefficient epsilon used in the Alfven test. The schemes introduce no new physical entities. The central claim depends on the unshown stability analysis and on the assumption that the explicitly treated source terms are non-stiff.

free parameters (2)
  • MIRK coefficients c1...c6 = 1st order: c1=0, c2=0, c3=1; 2nd order: c1=-1/sqrt(2), c2=1-sqrt(2)/2, c3=1, c4=(1-c1)^2/(2c1), c5=(sqrt(2)-1)/2, c6=1/2
    Chosen via a linear stability analysis at infinite conductivity (Section 4); presented as derived, not fitted to test data.
  • Kreiss-Oliger dissipation coefficient epsilon = 0.01 (0.1 for convergence runs)
    Artificial dissipation added to the hydrodynamic sector in the Alfven wave test (Eq. 36); value chosen by hand, affects stability and convergence measurements.
assumptions (4)
  • domain assumption Komissarov's augmented RRMHD system with scalar potentials phi, psi and damping kappa (Eqs. 1-5) is the correct evolution system.
    Adopted from [12] without derivation; the method is built on this model.
  • domain assumption Ohm's law with J^i = sigma*W*[E^i + (v x B)^i - (E^j v_j)v^i] + q v^i.
    Assumed in Section 2; all source terms in the stiff part derive from this.
  • domain assumption The recovery procedure of [15] using a quartic for W and bisection gives the unique primitive variables.
    Used in the Alfven wave test at every stage; failure of convergence would break the simulation.
  • domain assumption In the ideal limit sigma->infinity, E^i = -(v x B)^i.
    Used in Section 4 to set the zero-eigenvalue condition and select coefficients; physical limit of the model.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Numerical evolution of the resistive relativistic magnetohydrodynamic equations: a minimally implicit Runge-Kutta scheme." pith.science (2026). https://pith.science/paper/WBIL6HFI

@misc{pith2026250200990,
  author       = {Pith},
  title        = {Pith review of: Numerical evolution of the resistive relativistic magnetohydrodynamic equations: a minimally implicit Runge-Kutta scheme},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/WBIL6HFI}},
  note         = {Machine review of arXiv:2502.00990}
}
read the original abstract

We present the Minimally-Implicit Runge-Kutta (MIRK) methods for the numerical evolution of the resistive relativistic magnetohydrodynamic (RRMHD) equations, following the approach proposed by Komissarov (2007) of an augmented system of evolution equations to numerically deal with constraints. Previous approaches rely on Implicit-Explicit (IMEX) Runge-Kutta schemes; in general, compared to explicit schemes, IMEX methods need to apply the recovery (which can be very expensive computationally) of the primitive variables from the conserved ones in numerous additional times. Moreover, the use of an iterative process for the recovery could have potential convergence problems, increased by the additional number of required loops. In addition, the computational cost of the previous IMEX approach in comparison with the standard explicit methods is much higher. The MIRK methods are able to deal with stiff terms producing stable numerical evolutions, minimize the number of recoveries needed in comparison with IMEX methods, their computational cost is similar to the standard explicit methods and can actually be easily implemented in numerical codes which previously used explicit schemes. Two standard numerical tests are shown in the manuscript.

Figures

Figures reproduced from arXiv: 2502.00990 by the authors.

Figure 1
Figure 1. Self similar current sheet test. Numerical values [PITH_FULL_IMAGE:figures/full_fig_p011_1.png] view at source ↗
Figure 2
Figure 2. Self similar current sheet test. Numerical values [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗
Figure 3
Figure 3. Self similar current sheet test. Zoom on the numeri [PITH_FULL_IMAGE:figures/full_fig_p013_3.png] view at source ↗
Figures from the paper (9 more)
Figure 4
Figure 4. Figure 4: Self similar current sheet test. Numerical result [PITH_FULL_IMAGE:figures/full_fig_p014_4.png]
Figure 5
Figure 5. Figure 5: Self similar current sheet test. Numerical values [PITH_FULL_IMAGE:figures/full_fig_p015_5.png]
Figure 6
Figure 6. Figure 6: Self similar current sheet test. Numerical values [PITH_FULL_IMAGE:figures/full_fig_p015_6.png]
Figure 7
Figure 7. Figure 7: Self similar current sheet test. Numerical values [PITH_FULL_IMAGE:figures/full_fig_p016_7.png]
Figure 8
Figure 8. Figure 8: Self similar current sheet test. Numerical result [PITH_FULL_IMAGE:figures/full_fig_p016_8.png]
Figure 9
Figure 9. Figure 9: CP Alfv´en waves test. Numerical values for [PITH_FULL_IMAGE:figures/full_fig_p018_9.png]
Figure 10
Figure 10. Figure 10: CP Alfv´en waves test. Numerical values for [PITH_FULL_IMAGE:figures/full_fig_p019_10.png]
Figure 11
Figure 11. Figure 11: CP Alfv´en waves test. Zoom on the numerical value [PITH_FULL_IMAGE:figures/full_fig_p020_11.png]
Figure 12
Figure 12. Figure 12: CP Alfv´en waves test. Numerical results for seve [PITH_FULL_IMAGE:figures/full_fig_p022_12.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

18 extracted references · 4 canonical work pages

  1. [1]

    S. A. Balbus, J. F. Hawley, Instability, turbulence, and enhance d trans- port in accretion disks, Reviews of Modern Physics 70 (1998) 1–53. doi:10.1103/RevModPhys.70.1

    doi:10.1103/revmodphys.70.1 1998

  2. [2]

    O. M. Blaes, A. S. Balbus, Local shear instabilities in weakly ionized, weakly magnetized disks, The Astrophysical Journal 421 (1994) 1 63. doi:10.1086/173634

    doi:10.1086/173634 1994

  3. [3]

    J.-P. D. Villiers, J. F. Hawley, H. H. Krolik, Magnetically driven accre tion flows in the kerr metric. i. models and overall structure, The Astro physical Journal 599 (2003) 1238–1253. doi:10.1086/379509

    doi:10.1086/379509 2003

  4. [4]

    P. C. Fragile, O. M. Blaes, P. Anninos, J. D. Salmonson, Global gen eral rel- ativistic magnetohydrodynamic simulation of a tilted black hole accret ion disk, The Astrophysical Journal 668 (2007) 417–429. doi:10.1086/521092. 24

    doi:10.1086/521092 2007

  5. [5]

    Gabler, P

    M. Gabler, P. Cerd´ a-Dur´ an, J. A. Font, E. Mueller, N. Stergio ulas, Magneto-elastic oscillations and the damping of crustal shear mode s in magnetars, Monthly Notices of the Royal Astronomical Society 41 0 (2011) L37–L41. doi:10.1111/j.1745-3933.2010.00974.x

    arXiv 2011

  6. [6]

    V. M. Kaspi, A. M. Beloborodov, Magnetars, Annual Re- view of Astronomy and Astrophysics 55 (2017) 261–301. doi:10.1146/annurev-astro-081915-023329

    doi:10.1146/annurev-astro-081915-023329 2017

  7. [7]

    Komissarov, O

    S. Komissarov, O. Porth, Numerical simulations of jets, New Ast ronomy Reviews 92 (2021) 101610. doi:10.1016/j.newar.2021.101610

    arXiv 2021

  8. [8]

    J. M. Mart ´ ı, Numerical simulations of jets from active galactic nu clei, Galaxies 7 (2019) 24. doi:10.3390/galaxies7010024

    work page doi:10.3390/galaxies7010024 2019

Show all 18 references
  1. [9]

    Perucho, Dissipative processes and their role in the evolution o f radio galaxies, Galaxies 7 (2019) 70

    M. Perucho, Dissipative processes and their role in the evolution o f radio galaxies, Galaxies 7 (2019) 70. doi:10.3390/galaxies7030070

    2019 doi

  2. [10]

    H. C. Spruit, Dynamo action by differential rotation in a stably st rat- ified stellar interior, Astronomy and Astrophysics 381 (2002) 923– 932. doi:10.1051/0004-6361:20011465

    2002 doi

  3. [11]

    Palenzuela, L

    C. Palenzuela, L. Lehner, O. Reula, L. Rezzolla, Beyond ideal MH D: towards a more realistic modelling of relativistic astrophysical plasma s, Monthly Notices of the Royal Astronomical Society 394 (2009) 172 7–1740. doi:10.1111/j.1365-2966.2009.14454.x

    2009

  4. [12]

    S. S. Komissarov, Multidimensional numerical scheme for resist ive relativis- tic magnetohydrodynamics, Monthly Notices of the Royal Astrono mical Society 382 (2007) 995–1004. doi:10.1111/j.1365-2966.2007.12448.x

    2007

  5. [13]

    Ant´ on, O

    L. Ant´ on, O. Zanotti, J. A. Miralles, J. M. Mart ´ ı, J. M. Ib´ a˜ nez, J. A. Font, J. A. Pons, Numerical 3+1 general relativistic magnetohydrodyna mics: a local characteristic approach, The Astrophysical Journal 637 ( 2006) 296. 25

    2006

  6. [14]

    J. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerd´ a-Dur´ an, Comp utational general relativistic force-free electrodynamics, Astronomy and Astrophysics 647 (2021) A57. doi:10.1051/0004-6361/202038907

    2021 doi

  7. [15]

    Dumbser, O

    M. Dumbser, O. Zanotti, Very high order pnpm schemes on unstructured meshes for the resistive relativistic mhd equations, Journal of Com puta- tional Physics 228 (2009) 6991–7006. doi:10.1016/j.jcp.2009.06.009

    2009 doi

  8. [16]

    Cordero-Carri´ on, P

    I. Cordero-Carri´ on, P. Cerd´ a-Dur´ an, Partially implicit runge-kutta meth- ods for wave-like equations (arxiv:1211.5930) (2012)

    2012 arXiv

  9. [17]

    Cordero-Carri´ on, P

    I. Cordero-Carri´ on, P. Cerd´ a-Dur´ an, Partially implicit runge-kutta meth- ods for wave-like equations, Advances in Differential Equations and Appli- cations, SEMA SIMAI Springer Series (Springer International Pub lishing Switzerland, Switzerland) 4 (2014)

    2014

  10. [18]

    palenzuela

    Private discussion with c. palenzuela. 26 -1 -0.8 -0.6 -0.4 -0.2 0 x -4 -3 -2 -1 0 1 2 3 4 5 order taking first two resolutions taking last two resolutions 0 0.2 0.4 0.6 0.8 1 x 0 0.5 1 1.5 2 2.5 3 3.5order 1st order MIRK method 2nd order MIRK method 0.52 0.525 0.53 0.535 0.54...

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.