REVIEW 2 major objections 4 minor 2 cited by
A Particle Module for the PLUTO Code: III -- Dust
T0 review · 2 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read A dust integrator based on the exponential midpoint rule removes the stiffness constraint in dust-gas simulations.
desk verdict Dust module for PLUTO with exponential midpoint pusher and curvilinear weights: solid, well-benchmarked, but the 'arbitrary stopping time' claim only covers the particle push, not gas back-reaction. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the exponential midpoint integrator, a drift-kick-drift particle update in which the kick step solves the linear drag term exactly through the variation-of-constants formula with $A = -\operatorname{diag}(1/\tau_s)$ and propagator $h_1 = \tau_s(1 - e^{-\Delta t/\tau_s})$. Because the stiff part of the equation is integrated exactly, the time step is no longer limited by the stopping time; the remaining nonlinear forces are evaluated at the half-time level, giving second-order accuracy and time-reversibility. The secondary object is the set of curvilinear weighting factors obtained from $b$-spline shape functions with a normalization correction, which replace volume-coordinate interpolation and reduce grid noise in cylindrical and spherical domains.
What would settle it
Run the homogeneous deceleration test of Section 4.1 with $\tau_s \ll \Delta t$ (for example $\tau_s=0.001$, $\Delta t=1$) and monitor the gas velocity: the particle velocity decays to the predicted asymptotic value, but if the gas is not accelerated at the rate given by Eq. (50) with the full $1+\epsilon$ drag, the back-reaction part of the central claim is not valid.
Extended reading notes
Core claim
On its own terms, the paper claims that the exponential midpoint rule (Eq. 32) is a superior dust pusher: unlike the semi-implicit scheme of BS10, it does not oscillate and does not under-damp when the stopping time is much smaller than the time step, and it has the correct asymptotic velocity in the stiff limit. Together with the derived curvilinear CIC/TSC weighting factors, this yields a dust-gas scheme with back-reaction that works in Cartesian, cylindrical, and spherical geometries, enabling global protoplanetary-disk simulations in which dust and gas exchange momentum through drag.
Load-bearing premise
The claim of arbitrary stopping-time stability applies to the dust particle push; the gas-dust feedback still uses $\max(\tau_s,\Delta t)$ in the coupling denominator, so when the hydro time step exceeds the stopping time the momentum transferred from tightly coupled grains to the gas is deliberately reduced rather than physical.
Editorial extensions
If this is right
- Tightly coupled dust can be evolved with the same hydro time step, removing the stiff-drag time-step restriction from the particle component.
- Global disk simulations in cylindrical and spherical geometry can include dust back-reaction, not just test particles.
- The streaming instability can be captured in local shearing-box and global disk models, with linear growth rates matching the predicted values at sufficient resolution.
- The new TSC weighting reduces interpolation noise in curvilinear geometries, so fewer particles per cell are needed to represent a smooth dust density.
- Since the exponential scheme reduces to the semi-implicit method as $\tau_s \to \infty$, existing geometric-orbit behavior is retained for loosely coupled grains.
Reading between the lines
- The paper's stiffness claim is about the particle equation; for dust-to-gas momentum feedback in the regime $\tau_s < \Delta t$, the effective drag is capped by the hydro time step, so global simulations with very small grains will underestimate the back-reaction unless the gas step is reduced.
- The same exponential-midpoint construction could be applied to other stiff linear couplings in particle methods, such as Lorentz-force terms or cooling source terms.
- A direct comparison of the exponential pusher with the BS10 semi-implicit pusher in nonlinear streaming-instability runs could show whether removing the oscillatory under-damping changes clump statistics at low Stokes numbers.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper presents a new dust-particle module for the PLUTO code, using an exponential midpoint integrator for the dust push and modified PIC weighting kernels for curvilinear grids. The central claims are that the exponential midpoint rule is second-order, time-reversible, has bounded energy errors, and remains stable and asymptotically correct for arbitrarily small stopping times, and that the resulting hybrid scheme supports dust back-reaction in Cartesian, cylindrical, and spherical geometries. The paper validates the implementation with analytic deceleration and epicyclic tests, orbital tests in cylindrical and spherical coordinates, radial-drift and vertical-oscillation tests, a rigidly rotating disk with feedback, and local and global streaming-instability simulations. The global streaming-instability growth rates at the highest resolution (6144x384) approach the linear-theory values from Kowalik et al. (2013).
Significance. If the claims were fully established, the module would be a valuable community tool: it would allow global protoplanetary-disk simulations to evolve tightly coupled dust without the Δt ≲ τ_s restriction, using conservative feedback and improved curvilinear weighting. The benchmark suite is genuinely extensive, including analytic reference solutions for deceleration, epicyclic motion, orbital motion, radial drift, vertical settling, and the rigid-disk feedback problem, and the streaming-instability growth rates are checked against external linear-theory results rather than against outputs of the same model. The exponential midpoint update itself is well motivated and clearly superior to the semi-implicit scheme in the stiff single-particle limit. However, as detailed below, the load-bearing claim about arbitrarily small stopping times is not established for the coupled gas-dust system with back-reaction, because the gas-side feedback is capped in the predictor step.
major comments (2)
- [Abstract; §3.1 (Eq. 19); §3.2.2 (Eq. 32)] The claim that the exponential midpoint integrator 'remains stable in the limit of arbitrarily small particle stopping times yielding the correct asymptotic solution' is only demonstrated for the single-particle push, not for the coupled gas-dust system with back-reaction. In the predictor step, Eq. (19) replaces τ_s by max(τ_s, Δt), following BS10, so for τ_s << Δt the gas receives a drag acceleration of order Δv/Δt rather than Δv/τ_s. The corrector, Eq. (21), conserves total momentum but cannot repair the half-time gas velocity used to push the particles, so the external impulse is partitioned incorrectly between the two species. A concrete homogeneous example makes this explicit: with a constant acceleration F acting only on the dust, ε = ρ_D/ρ_g = 1, and Δt >> τ_s, the exact coupled solution after one step has v_g ≈ v_D ≈ FΔt/2 and relative velocity O(Fτ_s), whereas the algorithm as described gives v_g ≈ FΔt and v_D ≈ Fτ_s, i.e., a relative velocity O(FΔt). This directly limits the headline claim to test-particle dynamics or to Δt ≲ τ_s whenever back-reaction is dynamically important. The authors' §5 lists fluid-side stiffness as an open concern, but the abstract and summary do not carry this restriction.
- [§4 (Benchmarks)] No benchmark in Section 4 exposes the stiff back-reaction error identified above. The deceleration test in §4.1 is stiff (τ_s = 0.02 with Δt = 1) but has no external acceleration, so total momentum conservation alone can hide the incorrect gas-dust momentum partition. The rigid-disk feedback test in §4.7 uses τ_s = 1 with a hydro timestep set by the CFL condition, so the feedback term is not stiff, and the global streaming-instability test in §4.8 uses St ≈ 1.2. I ask for a stiff coupled benchmark that includes an external or background acceleration with τ_s << Δt and back-reaction, such as a homogeneous two-fluid terminal-velocity problem or a radial-drift setup with very small Stokes number. The abstract and Section 5 should be reworded so that the 'arbitrarily small stopping time' claim is restricted to the particle update, or the gas-side integration should be made stiff-accurate as well.
minor comments (4)
- [§4.4] The text says computations use time step sizes Δt = 0.1, 0.02, 0.04, but the panels and convergence discussion refer to 0.1, 0.01, and 0.001; the listed values are inconsistent.
- [§4.4] For the tilted orbits, the sentence 'the time steps have been halved, that is, Δt = 5×(10^{-2}, 10^{-3}, 10^{-3})' contains a duplicated entry; the intended sequence is presumably 5×10^{-2}, 5×10^{-3}, 5×10^{-4}.
- [§5] There is a typo in the summary: 'ensamble' should read 'ensemble'.
- [§4.1, Eq. (52)] The error definition err = (min_p err_p + max_p err_p)/2 is unusual; please state why the midpoint of the extremal particle errors is used rather than, for example, the mean or L2 norm over particles.
Circularity Check
No significant circularity: the exponential integrator and curvilinear weights are derived from exact ODE/geometry principles and validated against external analytic and linear-theory benchmarks.
full rationale
The paper's central claims are supported by independent external checks: the particle-gas deceleration test is compared to the closed-form solution Eq. (51); the epicyclic and orbital tests use analytic Keplerian/epicyclic solutions; the local streaming instability is measured against the linear eigenmodes of Youdin & Goodman (2005) as tabulated in YJ07 and BS10; and the global streaming instability uses the independent linear dispersion solver of Kowalik et al. (2013). The exponential midpoint integrator's stiff-limit behavior follows from the exact variation-of-constants formula (Eq. 29), so its asymptotic correctness is a mathematical property of the construction, not a fitted input disguised as a prediction. The new curvilinear weighting factors are defined by projecting b-spline shape functions (Eqs. 38-46) onto grid volumes and are tested against traditional volume weighting and analytic drift/settling solutions; their derivation is definitional in the benign sense of constructing a scheme, not circular. Self-citations to Paper I and the PLUTO code inherit the PIC deposition framework and hydro solver infrastructure, but the load-bearing novelty (exponential midpoint pusher and curvilinear weights) does not reduce to those citations. The paper itself notes in Section 3.1 and Section 5 that fluid-side stiffness remains an open issue because Eq. (19) caps the stopping time by max(τ_s, Δt); this narrows the 'arbitrarily small stopping time' claim to the dust push and is a genuine correctness caveat, but it is not a circular step. No equation in the manuscript is equivalent to its inputs by construction in a way that fakes a prediction.
Assumptions & free parameters
assumptions (6)
- domain assumption Locally isothermal equation of state p = rho c_s^2 and neglect of the gas energy equation; drag work and frictional heating are not evolved.
- domain assumption Epstein drag regime with linear drag force f_D = -rho_D (v_g - v_D)/tau_s.
- ad hoc to paper The gas back-reaction stiffness is controlled by replacing tau_s with max(tau_s, dt) in Eq. (19).
- standard math Standard theory of exponential integrators (variation of constants) and Störmer-Verlet geometric integration.
- domain assumption Shearing-box approximation with q = 3/2 and the orbital advection (FARGO) decomposition.
- domain assumption Global streaming instability test uses the local shearing-sheet dispersion relation and neglects vertical stratification and radial boundary effects.
Cite this review
Pith. "Pith review of A Particle Module for the PLUTO Code: III -- Dust." pith.science (2026). https://pith.science/paper/HCQUAPT5
@misc{pith2026190810793,
author = {Pith},
title = {Pith review of: A Particle Module for the PLUTO Code: III -- Dust},
year = {2026},
howpublished = {\url{https://pith.science/paper/HCQUAPT5}},
note = {Machine review of arXiv:1908.10793}
}
read the original abstract
The implementation of a new particle module describing the physics of dust grains coupled to the gas via drag forces is the subject of this work. The proposed particle-gas hybrid scheme has been designed to work in Cartesian as well as in cylindrical and spherical geometries. The numerical method relies on a Godunov-type second-order scheme for the fluid and an exponential midpoint rule for dust particles which overcomes the stiffness introduced by the linear coupling term. Besides being time-reversible and globally second-order accurate in time, the exponential integrator provides energy errors which are always bounded and it remains stable in the limit of arbitrarily small particle stopping times yielding the correct asymptotic solution. Such properties make this method preferable to the more widely used semi-implicit or fully implicit schemes at a very modest increase in computational cost. Coupling between particles and grid quantities is achieved through particle deposition and field-weighting techniques borrowed from Particle-In-Cell simulation methods. In this respect, we derive new weight factors in curvilinear coordinates that are more accurate than traditional volume- or area-weighting. A comprehensive suite of numerical benchmarks is presented to assess the accuracy and robustness of the algorithm in Cartesian, cylindrical and spherical coordinates. Particular attention is devoted to the streaming instability which is analyzed in both local and global disk models. The module is part of the PLUTO code for astrophysical gas-dynamics and it is mainly intended for the numerical modeling of protoplanetary disks in which solid and gas interact via aerodynamic drag.
Figures
Figures from the paper (10 more)
Forward citations
Cited by 2 Pith papers
-
Developing a Non-Newtonian Fluid Model for Dust, for Application to Astrophysical Flows
Collisionless dust in turbulent gas is derived as a 6D anisotropic Maxwell fluid whose rheological stress tensor is dynamically important in accretion discs.
-
A Staggered Semi-Analytic Method for Simulating Dust Grains Subject to Gas Drag
SSA, a staggered semi-analytic integrator, tracks dust under linear gas drag with 2nd-order accuracy in the non-stiff regime and near-exact terminal velocity behavior in the stiff regime, allowing time steps up to 10^...
Reference graph
Works this paper leans on
-
[1]
L., P \'e rez , L
ALMA Partnership , Brogan , C. L., P \'e rez , L. M., et al. 2015, , 808, L3
2015
-
[2]
M., Huang , J., P \'e rez , L
Andrews , S. M., Huang , J., P \'e rez , L. M., et al. 2018, , 869, L41
2018
-
[3]
P., Garufi , A., et al
Avenhaus , H., Quanz , S. P., Garufi , A., et al. 2018, , 863, 44
2018
-
[4]
Bai , X.-N., & Stone , J. M. 2010, , 190, 297
work page 2010
-
[5]
Balsara , D. S., Tilley , D. A., Rettig , T., & Brittain , S. D. 2009, , 397, 24
work page 2009
-
[6]
Ben \' tez-Llambay , P., Krapp , L., & Pessah , M. E. 2019, , 241, 25
2019
-
[7]
Birdsall, C., & Langdon, A. 2004, Plasma Physics via Computer Simulation, Series in Plasma Physics and Fluid Dynamics (Taylor & Francis)
work page 2004
- [8]
Show all 52 references
-
[9]
1990, Journal of Computational Physics, 87, 171
Colella , P. 1990, Journal of Computational Physics, 87, 171
1990
-
[10]
M., & Matthews , P
Cox , S. M., & Matthews , P. C. 2002, Journal of Computational Physics, 176, 430
2002
-
[11]
Gottlieb , S., & Shu , C. W. 1998, Mathematics of Computation, 67, 73
1998
-
[12]
2006, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations; 2nd ed
Hairer, E., Lubich, C., & Wanner, G. 2006, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations; 2nd ed. (Dordrecht: Springer), doi:10.1007/3-540-30666-8
2006 doi
-
[13]
F., Gammie , C
Hawley , J. F., Gammie , C. F., & Balbus , S. A. 1995, , 440, 742
1995
-
[14]
1981, Progress of Theoretical Physics Supplement, 70, 35
Hayashi , C. 1981, Progress of Theoretical Physics Supplement, 70, 35
1981
-
[15]
2010, Acta Numerica, 19, 209
Hochbruck , M., & Ostermann , A. 2010, Acta Numerica, 19, 209
2010
-
[16]
2005, , 634, 1353
Johansen , A., & Klahr , H. 2005, , 634, 1353
2005
-
[17]
2007, , 662, 627
Johansen , A., & Youdin , A. 2007, , 662, 627
2007
-
[18]
2013, , 434, 1460
Kowalik , K., Hanasz , M., W \'o lta \'n ski , D., & Gawryszczak , A. 2013, , 434, 1460
2013
-
[19]
Laibe , G., & Price , D. J. 2014, , 440, 2136
2014
-
[20]
2012, Journal of Computational Physics, 231, 795
Lapenta , G. 2012, Journal of Computational Physics, 231, 795
2012
-
[21]
J., Hewett , D
Larson , D. J., Hewett , D. W., & Langdon , A. B. 1995, Computer Physics Communications, 90, 260
1995
-
[22]
2017, , 607, A74
Liu , Y., Henning , T., Carrasco-Gonz \'a lez , C., et al. 2017, , 607, A74
2017
-
[23]
Marble , F. E. 1970, Annual Review of Fluid Mechanics, 2, 397
1970
-
[24]
2012, , 545, A134
Meheut , H., Meliani , Z., Varniere , P., & Benz , W. 2012, , 545, A134
2012
-
[25]
M.-A., Mackey , J., Langer , N., et al
Meyer , D. M.-A., Mackey , J., Langer , N., et al. 2014, , 444, 2754
2014
-
[26]
M.-A., Mignone , A., Kuiper , R., Raga , A
Meyer , D. M.-A., Mignone , A., Kuiper , R., Raga , A. C., & Kley , W. 2017, , 464, 3229
2017
-
[27]
2014, Journal of Computational Physics, 270, 784
Mignone , A. 2014, Journal of Computational Physics, 270, 784
2014
-
[28]
2007, , 170, 228
Mignone , A., Bodo , G., Massaglia , S., et al. 2007, , 170, 228
2007
-
[29]
2018, , 859, 13
Mignone , A., Bodo , G., Vaidya , B., & Mattia , G. 2018, , 859, 13
2018
-
[30]
M., & Muscianisi , G
Mignone , A., Flock , M., Stute , M., Kolb , S. M., & Muscianisi , G. 2012 a , , 545, A152
2012
-
[31]
2012 b , , 198, 7
Mignone , A., Zanni , C., Tzeferacos , P., et al. 2012 b , , 198, 7
2012
-
[32]
2010, Journal of Computational Physics, 229, 3916
Miniati , F. 2010, Journal of Computational Physics, 229, 3916
2010
-
[33]
1986, , 67, 375
Nakagawa , Y., Sekiya , M., & Hayashi , C. 1986, , 67, 375
1986
-
[34]
P., Gressel , O., & Umurhan , O
Nelson , R. P., Gressel , O., & Umurhan , O. M. 2013, , 435, 2610
2013
-
[35]
2013, SIAM Journal on Scientific Computing, doi:10.1137/050635018
Pelanti, M., & Leveque, R. 2013, SIAM Journal on Scientific Computing, doi:10.1137/050635018
2013 doi
-
[36]
Picogna , G., Stoll , M. H. R., & Kley , W. 2018, , 616, A116
2018
-
[37]
Pinte , C., Dent , W. R. F., M \'e nard , F., et al. 2016, , 816, 25
2016
-
[38]
P., & Keppens , R
Porth , O., Xia , C., Hendrix , T., Moschou , S. P., & Keppens , R. 2014, , 214, 4
2014
-
[39]
Ruyten , W. M. 1993, Journal of Computational Physics, 105, 224
1993
-
[40]
2019, Journal of Computational Physics, 382, 27
Shen, X., & Leok, M. 2019, Journal of Computational Physics, 382, 27
2019
-
[41]
Stoll , M. H. R., & Kley , W. 2016, , 594, A57
2016
-
[42]
2001, , 557, 990
Takeuchi , T., & Artymowicz , P. 2001, , 557, 990
2001
-
[43]
2016, , 589, A10
Thun , D., Kuiper , R., Schmidt , F., & Kley , W. 2016, , 589, A10
2016
-
[44]
2018, , 865, 144
Vaidya , B., Mignone , A., Bodo , G., Rossi , P., & Massaglia , S. 2018, , 865, 144
2018
-
[45]
J., Meliani , Z., Keppens , R., & Decin , L
van Marle , A. J., Meliani , Z., Keppens , R., & Decin , L. 2011, , 734, L26
2011
-
[46]
Verboncoeur , J. P. 2001, Journal of Computational Physics, 174, 421
2001
-
[47]
Weidenschilling , S. J. 1977, , 180, 57
1977
-
[48]
Whipple , F. L. 1972, in From Plasma to Planet, ed. A. Elvius , 211
1972
-
[49]
2016, , 224, 39
Yang , C.-C., & Johansen , A. 2016, , 224, 39
2016
-
[50]
2007, , 662, 613
Youdin , A., & Johansen , A. 2007, , 662, 613
2007
-
[51]
N., & Goodman , J
Youdin , A. N., & Goodman , J. 2005, , 620, 459
2005
-
[52]
M., Rafikov , R
Zhu , Z., Stone , J. M., Rafikov , R. R., & Bai , X.-n. 2014, , 785, 122
2014
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.