{"id":"7d259f28-a2e0-487c-b947-80e3fd4986a4","arxiv_id":"2501.05917","paper_version":1,"verdict":"REJECT","confidence":"MODERATE","novelty_score":2.0,"correctness_risk":"high","formal_verification":"none","parameter_count":4,"one_line_summary":"An explicit finite difference solver for 3-D Poisson-Nernst-Planck equations is presented, with stability criteria that restrict reliable simulation to short time windows.","lead":"This paper applies an explicit finite difference scheme to solve 3-D Poisson-Nernst-Planck equations for ion transport in electrolytes, including Butler-Volmer surface reactions. The authors derive stability conditions and find that coupling the equations makes stable simulation difficult beyond a short time window.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The central stability criteria (Eqs. 30 and 36) rely on an invalid additive combination of amplification factors and a real-valued treatment of the imaginary migration term, so they do not constrain the actual explicit PNP scheme.","rationale":"The reader's weakest assumption identifies exactly the load-bearing weakness: the stability analysis in Section 2.1.4 treats the migration potential as a constant multiplier and combines the N-P and Poisson amplification factors by addition. My reading confirms this is not a minor technical slip but a fundamental breakdown of the von Neumann method for the coupled PNP system. The migration term, when Fourier-analyzed with a spatially varying potential, either fails to diagonalize (if Φ varies) or produces an imaginary amplification factor (if Φ is constant), yet the paper's stability conditions assume a real λ in [-1,1]. The additive composition in Eq. (43) has no valid basis for coupled time-stepping schemes. Therefore the derived inequalities (30) and (36), which form the paper's main practical output, are unsupported. The paper itself acknowledges unreliable results after 0.3 seconds and inconsistent time-step guidance, which is consistent with this flaw. Because the central claim of enabling precise simulation hangs on these stability criteria, and because no code, data, or benchmark validation is supplied to independently confirm the framework, the reject verdict is appropriate. A concrete linearized stability test would settle the matter definitively.","tokens_in":13128,"tokens_out":3683,"duration_ms":34128,"concrete_test":"Perform a linearized von Neumann analysis of the explicit PNP update (Eqs. 11-14) around a uniform steady state with a constant background electric field E0. Replace Φ by Φ0 - E0·x, linearize the migration flux in c, and compute the exact amplification matrix G(k) for the coupled (c, Φ) system (or the equivalent one-step operator if Φ is updated by Gauss-Seidel each step). Then check whether any eigenvalue of G has modulus greater than 1 for parameter sets that satisfy Eqs. (30) and (36), for example the paper's own parameters with Δx=0.1, D_z=1e-2, and Δt chosen so that the μ terms satisfy Eq. (30). If such an eigenvalue exists, the derived stability conditions are not sufficient. As a complementary check, run a 1-D two-ion PNP simulation with those parameters and observe whether high-wavenumber perturbations grow within a few time steps.","verdict_should_be":"REJECT","load_bearing_attack":"The paper's central deliverable is the explicit FDM stability criteria. Section 2.1.4 Fourier-analyzes the migration flux by substituting a Fourier mode for c while keeping Φ as a fixed grid function (Eqs. 20-26). This is not a von Neumann analysis: for a non-constant Φ, the finite-difference operator is not translation-invariant, so Fourier modes do not diagonalize it; if Φ is treated as a constant multiplier, the migration term contributes an imaginary amplification factor i μ_E sin(kΔx), not a real one. The subsequent 'total amplification factor' λ = λ_NP + λ_Φ (Eq. 43) is additive, but the coupled PNP system's amplification operator is the composition (or product) of the N-P and Poisson update operators, not the sum; even for scalar factors, |λ1+λ2|≤1 is neither necessary nor sufficient for the product to satisfy |λ|≤1. The Case 1/2 analysis (Eqs. 28-30) treats λ as real and within [-1,1], which fails once the migration term is included. Consequently, the inequalities (30) and (36) do not follow from the discretization, and the claim that 'it is more impossible to simulate successfully a real electrochemical system' is not established. This also undermines the recommended time-step and grid guidelines (e.g., 'the time step size cannot be less than 50' in Section 2.2.1) and thus the paper's stated ability to 'precisely simulate' concentration profiles.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes an explicit finite-difference scheme for solving the three-dimensional Poisson-Nernst-Planck (PNP) equations with Butler-Volmer boundary conditions. It derives stability inequalities for the Nernst-Planck diffusion and migration terms and for the Poisson equation, and uses these to argue that explicit PNP coupling is severely restricted, leading to practical guidance on time steps and grid sizes. The paper includes simulations of single-ion and two-ion (1:1 electrolyte) systems and discusses the influence of grid resolution and time-step count. Its central deliverable is the set of stability criteria in Eqs. (30), (36), and (43).","tokens_in":13443,"tokens_out":3077,"duration_ms":30268,"significance":"If the stability criteria were correct, they would provide quantitative guidance for choosing time steps and grid sizes in explicit FDM solutions of 3D PNP equations, a practically relevant problem in electrochemical modeling. The paper also demonstrates a step-by-step discretization and applies Butler-Volmer boundary conditions, which could be useful as a starting point for similar implementations. However, the central derivation is built on an invalid Fourier analysis of the migration term and an incorrect additive combination of amplification factors for the coupled system. Because the paper's main conclusions and practical recommendations rest on these steps, the contribution as presented is not reliable. The manuscript also contains internal inconsistencies in the boundary conditions and important presentation gaps, so the claimed 'precise simulation' capability is not substantiated.","major_comments":[{"comment":"The stability analysis of the migration term treats the electric potential Φ as a constant while substituting a Fourier mode for the concentration. For a non-constant Φ, the finite-difference migration operator is not translation-invariant, so Fourier modes do not diagonalize the update. If Φ is treated as a constant multiplier, the migration term contributes an imaginary amplification factor i μ_E sin(kΔx), not the real expression written in Eq. (26). Consequently, the real-valued derivation of cases 1 and 2 leading to inequalities (28)–(30) does not govern the actual discretized scheme.","section":"§2.1.4, Eqs. (20)–(26)"},{"comment":"The total amplification factor for the coupled PNP system is defined additively as λ_total = λ_N−P + λ_Φ. The coupled update is a composition of the Nernst-Planck and Poisson updates; even for scalar factors, |λ1 + λ2| ≤ 1 is neither necessary nor sufficient for stability of the product/composition. Thus the claim that the coupled stability condition is 'more restrictive' and the conclusion that 'it is more impossible to simulate successfully a real electrochemical system' do not follow from the analysis presented.","section":"§2.1.5.2, Eq. (43)"},{"comment":"The Poisson equation is elliptic and does not involve time stepping, yet an 'amplification factor' and a stability condition (36) are derived for it and then used as a constraint on the spatial grid. No derivation shows how this condition relates to the temporal stability of the coupled explicit scheme; absent such a link, Eq. (36) does not support the time-step recommendations made later in Section 2.2.1.","section":"§2.1.4.2, Eqs. (35)–(36)"},{"comment":"The boundary conditions for the concentration are internally contradictory: Eq. (44) imposes zero Neumann conditions at the lateral and top boundaries, while Eq. (45) imposes zero Dirichlet conditions (c = 0) at the same lateral boundaries. Both cannot hold simultaneously, and the boundary treatment used in the simulations is therefore not well-defined. Additionally, Eq. (46) mixes the exchange current density J0, concentration gradient, potential gradient, and mobility ν in a single expression whose dimensional consistency is not evident as printed.","section":"§2.1.6, Eqs. (44)–(45)"}],"minor_comments":[{"comment":"Equation (40) is garbled and missing operator symbols, making the Gauss-Seidel update formula hard to parse; it should be rewritten with clear placement of Φ and ρ terms.","section":"§2.1.5.1, Eq. (40)"},{"comment":"The phrase 'the time step size cannot be less than 50' is ambiguous and dimensionally incomplete; a time step size should be expressed with units (e.g., seconds), and this recommendation is not supported by the preceding stability analysis.","section":"§2.2.1"},{"comment":"The figures lack axis labels and color scales, and the text does not specify which spatial slice or time instant is shown. This makes it difficult to verify the claimed concentration profiles.","section":"Figures 3–6"},{"comment":"Equation (12) contains apparent typographical errors, such as a misplaced superscript and the notation 'c_c' and 'cα' mixed with 'cα' superscripts, which should be corrected for clarity.","section":"§2.1.3.3, Eq. (12)"},{"comment":"The references are inconsistently formatted (some with page ranges, some without), and the citations for the Nernst-Planck and Poisson equations [9,10] are to computational studies rather than to the original equation formulations; the authors should provide canonical references.","section":"References"}],"recommendation":"reject","confidential_remarks":"The manuscript's central stability criteria rest on a flawed Fourier analysis and an incorrect additive combination of amplification factors. These are not presentation issues but load-bearing errors: the main conclusions about the impossibility of explicit PNP coupling and the practical parameter recommendations would need to be rederived with a proper coupled-system analysis (e.g., matrix or energy methods). Given the scope of the paper, this amounts to a fundamental revision that goes beyond routine major revision. I therefore recommend rejection, though the authors may wish to resubmit a corrected analysis if they can provide a rigorous stability framework for the explicit PNP discretization."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The central stability analysis does not hold up, and the paper's main claim—that its stability criteria enable precise explicit 3D PNP simulation—is unsupported. The von Neumann treatment of the migration term treats Φ as a fixed multiplier, which is only valid for constant potential; for the non-uniform potentials relevant here, Fourier modes do not diagonalize the operator. Adding the N-P and Poisson amplification factors (Eq. 43) is not a valid way to analyze the coupled system; the update operators compose, not sum. The resulting inequalities (30) and (36) therefore don't follow.\n\nThere is something useful here. The paper is a clear, readable walkthrough of the FDM discretization for the 3D PNP equations with a Butler-Volmer boundary condition, and it correctly reproduces the textbook diffusion-only stability limit. The discussion of why explicit coupling becomes restrictive is a legitimate question. But the execution is too flawed to support a research claim. Boundary conditions are contradictory—Eq. (44) imposes zero concentration gradient at the lateral walls while Eq. (45) imposes zero concentration at the same walls. The authors' own simulations become unstable after 0.3 s, and the claim that 'the time step size cannot be less than 50' is incoherent. No code, data, or benchmark comparison is provided, so 'precise simulation' is not demonstrated. The conclusion that explicit PNP coupling is 'more impossible' is an overreach.\n\nThis reads like a methods tutorial that got ahead of itself. The topic—stability of explicit PNP schemes—is worth studying, and the paper cites appropriate prior FDM work. But as it stands, the load-bearing analysis is wrong, and I wouldn't want to referee it until that's fixed. I recommend a desk reject for this version, with an invitation to resubmit after a proper coupled stability analysis and actual validation.","headline":"The paper's stability criteria don't follow from its own analysis, and without them the central claim of a validated explicit 3D PNP solver collapses.","tokens_in":13937,"tokens_out":2526,"would_cite":false,"duration_ms":24719,"reading_group":"no","serious_thinker":"no","would_accept_peer_review":false},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper derives explicit finite-difference stability conditions for the 3-D Poisson-Nernst-Planck equations and shows that, under those conditions, explicit time-stepping of the coupled system is practically infeasible for realistic…","keywords":["Poisson-Nernst-Planck equations","finite difference method","concentration profile","stability analysis","3-D simulation","Butler-Volmer boundary condition","electrochemical transport","explicit time-stepping"],"falsifier":"Run an explicit 3-D PNP code with parameters satisfying Eqs. (30) and (36) but with a steep potential gradient near the electrode, and check whether oscillations or divergence appear before the nominal stable time horizon; if they do, the additive combination of amplification factors is disproven.","tokens_in":12894,"feed_emoji":"⚡","tokens_out":6683,"duration_ms":62490,"temperature":0.7,"pith_summary":"The paper tries to put 3-D Poisson-Nernst-Planck (PNP) simulation on firmer footing by giving a step-by-step finite-difference recipe and, crucially, the stability conditions that an explicit scheme must satisfy. It derives two inequalities, Eqs. (30) and (36), from Fourier analysis of the diffusion, migration, and Poisson updates, and argues that combining them makes explicit coupling of the full PNP system nearly impossible for realistic batteries or concentrated electrolytes. To get usable results, it treats the Poisson equation as an elliptic boundary-value problem solved by Gauss-Seidel iteration and applies Butler-Volmer kinetics at the electrode surface. If the stability analysis is right, modelers can check their time step and grid choice against closed-form limits before simulating, and will know why many existing 3-D PNP attempts diverge.","feed_headline":"Explicit 3-D electrolyte simulations need impractically strict limits","feed_subtitle":"Safe-step criteria from a finite-difference analysis rule out many realistic battery and electrochemical simulations.","key_machinery":"The machinery is the explicit finite-difference discretization of the 3-D Nernst-Planck and Poisson equations on a uniform cubic grid, together with Fourier-mode (von Neumann) stability analysis. The central objects are the amplification factors for the diffusion part, migration part, and Poisson part, whose combination yields the safe-step inequalities. Then the Poisson equation is removed from the time-stepping loop and solved as an elliptic problem by Gauss-Seidel iteration, which sidesteps the stricter Poisson stability constraint. Boundary conditions are handled with no-flux Neumann conditions on lateral and top faces and the Butler-Volmer electrode reaction at the bottom surface.","core_discovery":"The central discovery is a closed-form stability criterion for explicit finite-difference solution of the 3-D PNP system. By inserting Fourier modes into the discretized Nernst-Planck equation and treating the electric potential as a constant coefficient in the migration term, the paper obtains an amplification factor whose magnitude must stay below one; this yields Eq. (30), linking diffusion coefficients, time step, grid spacing, and migration strengths. A second Fourier analysis of the discretized Poisson update gives Eq. (36), a constraint on the grid spacing alone. Since the coupled-system amplification factor is taken as the sum of the two individual factors, the combined bound is stricter than either alone. The paper concludes that explicit time-stepping of the fully coupled PNP equations is generally not viable for realistic electrochemical systems, and that solving Poisson as an elliptic equation with Gauss-Seidel is the practical path; with that choice it produces 3-D concentration profiles for single-ion and 1:1 electrolyte systems under finite-layer and Butler-Volmer boundary conditions.","pith_inferences":["Because the stability derivation freezes the potential inside the migration term, the true stability limit in steep concentration or potential gradients is likely stricter than Eq. (30); a practical rule of thumb would be to stay well inside the predicted window.","The same amplification-factor construction could be adapted to implicit or semi-implicit time stepping, where the Poisson constraint would likely relax; this is a natural next test the paper does not run.","The reported instability at longer times, even inside the nominal stable parameter set, suggests the linear stability analysis misses nonlinear feedback between concentration and potential; checking against a highly resolved reference solution would quantify how much the practical window shrinks."],"forward_implications":["Anyone using an explicit finite-difference scheme for 3-D PNP can pre-check the time step and grid against Eqs. (30) and (36) and avoid runs that diverge.","Realistic battery parameters, such as high concentration, strong electric fields, and large domains, fall outside the stable window, which explains why fully explicit PNP coupling tends to fail in practice.","Treating the Poisson equation as an elliptic problem rather than time-stepping it is not just a convenience; under the paper's analysis it is the difference between feasible and infeasible simulation.","The parameter studies give practical guidance: roughly 50 to 80 time steps with moderate grid resolution balance stability and accuracy for the test cases.","The framework can serve as a benchmark or training-data generator for machine-learning models of ion transport, as the paper itself suggests."],"supporting_citations":[{"why":"Supplies the original 1-D Nernst-Planck current-distribution calculation in porous electrodes that the 3-D model extends.","marker":"[4]"},{"why":"Prior finite-difference PNP solver enforcing free-energy constraints; the paper's stability analysis targets the class of schemes this work represents.","marker":"[7]"},{"why":"Prior finite-difference treatment of a modified PNP equation with WKB comparison, used as a baseline for FDM-based PNP modeling.","marker":"[8]"},{"why":"Source of the Nernst-Planck flux equation that the discrete scheme is built from.","marker":"[9]"},{"why":"Source of the 3-D Poisson equation linking potential to charge density in the model.","marker":"[10]"}],"fun_headline_variants":["Explicit 3D PNP simulations face stability dead end","Stability criterion kills explicit 3D ion simulations","3D electrolyte modeling: explicit solvers impractical","Explicit 3D PNP: too strict for realistic electrochemistry","New stability bound dooms explicit 3D PNP solvers"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The proof assumes the electric potential can be treated as a fixed coefficient while Fourier-analyzing the migration term, and that the overall amplification factor is the sum of the Nernst-Planck and Poisson factors; if potential gradients vary sharply in space or time, those assumptions break and the safe-step inequalities no longer guarantee stability.","fun_headline_variants_meta":{"raw":{"variants":["Explicit 3D PNP simulations face stability dead end","Stability criterion kills explicit 3D ion simulations","3D electrolyte modeling: explicit solvers impractical","Explicit 3D PNP: too strict for realistic electrochemistry","New stability bound dooms explicit 3D PNP solvers"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000217,"raw_usage":{"total_tokens":1414,"prompt_tokens":899,"completion_tokens":515,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":515,"completion_tokens_details":{"reasoning_tokens":431}},"tokens_in":515,"tokens_out":515,"duration_ms":5121,"temperature":1.0,"reasoning_tokens":431,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T21:05:53.921806+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run an explicit 3-D PNP code with parameters satisfying Eqs. (30) and (36) but with a steep potential gradient near the electrode, and check whether oscillations or divergence appear before the nominal stable time horizon; if they do, the additive combination of amplification factors is disproven.","supporting_citations":[{"cited_title":"Newman, C.W","cited_arxiv_id":null,"evidence_quote":"Supplies the original 1-D Nernst-Planck current-distribution calculation in porous electrodes that the 3-D model extends."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Prior finite-difference PNP solver enforcing free-energy constraints; the paper's stability analysis targets the class of schemes this work represents."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Prior finite-difference treatment of a modified PNP equation with WKB comparison, used as a baseline for FDM-based PNP modeling."},{"cited_title":"Dolatabadi, A","cited_arxiv_id":null,"evidence_quote":"Source of the Nernst-Planck flux equation that the discrete scheme is built from."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Source of the 3-D Poisson equation linking potential to charge density in the model."}],"review_version":1}