{"id":"86f61f39-e01e-400c-a5f5-3e00c6546a21","arxiv_id":"1908.07975","paper_version":4,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A staggered finite difference scheme, based on half-step shifts and an interpolated relaxation term, yields stable second-order accurate simulations of Poynting-Thomson-Zener waves and is much faster than COMSOL for the elastic case.","lead":"This paper builds a new type of fast numerical recipe for waves in soft solids, by placing the grid points for different physical quantities half a step apart in space and time. The goal is to make simulations of rheological materials, like rocks or plastics, both accurate and cheap to run.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The advertised stability at the optimal setting ^C=1 depends on an unproven boundary-condition extension of the von Neumann analysis; the exceptional mode kΔx=π is not shown to be excluded by the pulse/free-boundary problem.","rationale":"The central deliverable is a finite-difference scheme whose usefulness rests on the stability and accuracy claims in Sections 4–6. The most load-bearing condition is that the scheme is stable at the recommended operating point ^C=1, because every PTZ demonstration (Section 5.2) uses that value. The von Neumann analysis proves only ^C<1 for the infinite periodic problem; equality is passed with a boundary-condition rule-of-thumb. The Hooke case exposes a concrete reason equality is delicate: the mode kΔx=π is a defective eigenvalue at C=1, and stability then hinges on whether the boundary conditions exclude it. No such exclusion is shown for the pulse/free-boundary problem, and the paper explicitly labels the extension as a rule-of-thumb rather than a theorem. I did not find an independent fatal flaw. The second-order claim for α=1/2 in the PTZ case is not written out, but the discretization (32) is the trapezoidal/midpoint rule for the rheological ODE and is locally second-order in the same sense as the two displayed proofs; the omission is exposition, not a detected error. The COMSOL comparison is about run time and qualitative behavior in one elastic case; it does not carry the central correctness claim. The parameter count is one and the scheme is simple enough that the stability question is the main risk. Since the identified concern is a gap in proof rather than a demonstrated failure, and since extensive numerical experiments at ^C=1 appear stable, the reader's CONDITIONAL verdict should stand unchanged.","tokens_in":17859,"tokens_out":14704,"duration_ms":141137,"concrete_test":"Perform a normal-mode (Godunov–Ryabenkii) stability analysis of the finite-interval problem with the actual boundary conditions (stress pulse (68) at x=0, zero stress at x=X), at ^C=1, for both the Hooke and PTZ versions. Specifically, substitute the kΔx=π mode into the boundary conditions and compute the determinant of the boundary-mode matrix; if the determinant vanishes, the mode is admissible and the scheme is only marginally stable, contradicting the ^C≤1 claim. If instead all boundary modes satisfy |ξ|≤1 (with no Jordan-block growth), the rule-of-thumb is validated for this problem.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's stability proof is a von Neumann analysis for plane waves on an infinite periodic domain. The transition from the strict condition ^C<1 (Eq. (60)) to the claimed ^C≤1 (Eq. (67)) is justified only by a stated rule-of-thumb (Section 4, after Eq. (43)): that boundary conditions cannot destabilize modes that are stable on the infinite line, and may only remove bad modes. This is exactly the step the main simulations rely on, since Section 5.2 runs the PTZ scheme at ^C=1. The Hooke analysis shows why the step is non-trivial: at C=1 the mode S=1 (kΔx=π) has a double eigenvalue ξ=-1 with geometric multiplicity one (Section 4.1, after Eq. (44)); the authors note that boundary conditions may prohibit this mode but do not analyze whether the pulse/free-boundary problem does. For PTZ at ^C=1 the strict Jury inequality (60) fails at S=1, and the paper's extension to equality is not derived. If that exceptional mode is admissible for the actual boundary conditions, or if the finite interval admits a weakly unstable boundary mode, the headline stability guarantee for the reported simulations is unsupported. The paper itself flags this as a rule-of-thumb, so the gap is explicit rather than hidden.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents a staggered finite difference scheme for the one-dimensional Poynting–Thomson–Zener (PTZ) rheological model, building on the symplectic Euler method and extending it to a dissipative continuum system. The authors derive the scheme from a spacetime-staggered arrangement of stress, strain, and velocity, analyze its von Neumann stability for both the Hooke and PTZ cases, and study dissipative and dispersive errors analytically and numerically. They then compare the scheme with COMSOL finite element simulations on a Hookean wave-propagation problem, reporting large run-time and accuracy advantages. The central claims are that the scheme is second-order accurate for α=1/2, that stability is governed by the Courant condition based on the fastest wave speed, and that with α=1/2 and the Courant number at the boundary of the stability region the scheme has linear dispersion branches and suppressed dissipative error.","tokens_in":18105,"tokens_out":1915,"duration_ms":21101,"significance":"If the claims hold, the scheme is a useful, simple, and fast alternative to standard finite element tools for linear rheological wave problems, with the attractive feature that its design follows from the spacetime structure of the governing equations rather than from fitted parameters. The paper contributes a clean, self-contained derivation: no constants are fitted, the Hooke-case stability analysis is rigorous and explicit, and the PTZ stability conditions are derived in closed form via both Jury and Routh–Hurwitz criteria. The explicit demonstration that the thermodynamic condition τ̂>τ emerges from the numerical stability analysis is a nice conceptual point. The COMSOL comparison, though limited to a linear elastic setup, is concrete and reproducible in its setup.","major_comments":[{"comment":"The transition from the strict stability condition Ĉ<1 derived in Eq. (60)/(66) to the claimed Ĉ≤1 in Eq. (67) is not proven. The paper explicitly labels this as a rule-of-thumb extension: the von Neumann analysis treats plane waves on an infinite periodic domain, while the actual simulations in Section 5.2 use a finite sample with a stress pulse at one end and a free boundary at the other. The exceptional mode with S=1 (kΔx=π) is precisely where the Jury inequality (60) fails, and the paper does not show that this mode is inadmissible for the pulse/free-boundary problem. Since the headline numerical results are run at Ĉ=1, the advertised stability of the scheme for the reported boundary-value problem rests on an unproven assumption. Please either prove that the exceptional mode is excluded by the boundary conditions, perform a boundary-mode analysis, or revise the stability claim to Ĉ<1 and rerun the key simulations accordingly.","section":"Section 4.2.2, Eq. (66)–(67)"},{"comment":"The claim that α=1/2 renders the PTZ update (33) second-order accurate is asserted without proof. The sentence \"Second order accuracy of (33) for α=1/2 is then straightforward to verify\" is not sufficient, especially because the update combines a finite difference ratio with an interpolation in σ and ε. The paper's stated advantage over the first-order symplectic Euler method depends on this accuracy claim, so a local truncation error derivation for the rheological update should be included or explicitly referenced. Without it, the accuracy comparison in Section 5.2 and the error discussion in Section 6.2 lack a rigorous basis.","section":"Section 3, Eq. (32)–(33)"}],"minor_comments":[{"comment":"There is a typo in the sentence \"this affects only one mode, S=1, k=π/k\"; this should read \"kΔx=π\".","section":"Section 4.1, after Eq. (46)"},{"comment":"The caption of Figure 8 does not indicate whether the upper and lower rows correspond to C=1 and C=1/2, respectively, as stated in the text; the caption should be self-contained.","section":"Section 6.1, Figure 8"},{"comment":"The comparison with COMSOL would be more convincing if the exact COMSOL settings (mesh element order, solver tolerances, time-stepping parameters) were listed in a table, since the runtime differences depend strongly on these choices.","section":"Section 7"},{"comment":"Reference [21] is cited as \"in preparation\" and \"under review\"; if the manuscript is being finalized, this reference should either be updated to a published version or removed as a support for the claim about dynamic versus static moduli.","section":"References"}],"recommendation":"major_revision","confidential_remarks":"The paper is methodologically interesting and the central scheme is likely sound, but the unproven boundary-condition extension from Ĉ<1 to Ĉ≤1 is load-bearing for the main numerical claims. The second-order accuracy assertion for α=1/2 is also missing a proof. Both issues are fixable within the manuscript's scope, hence major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Plainly: this paper earns a careful read. It is a real extension of the authors' prior staggered scheme to a dispersive rheological model, and the stability analysis for the periodic problem is done properly. The explicit Jury/Routh-Hurwitz conditions, especially the emergence of tau-hat > tau from the numerical stability requirement, are a nice touch. The dissipative and dispersive error discussion is also clear—the root plots tell a coherent story about why C-hat=1 is special.\n\nThe weaknesses are proportionate. The most important is the boundary-condition step. The strict condition C-hat<1 is derived, but the simulations run at C-hat=1, and the leap to equality relies on the stated rule-of-thumb that boundary conditions only remove bad modes. The paper identifies this as a rule-of-thumb, so the gap is out in the open. Still, for a claim of stability at the optimal setting, a referee should ask for either a proof that the pulse/free-boundary IBVP excludes the defective Nyquist mode, or numerical evidence at slightly below and at C-hat=1. The second-order accuracy of the alpha=1/2 PTZ scheme is asserted rather than demonstrated—the Taylor expansion is only shown for the Hooke case. I suspect the claim is true, but it's not verified in the text. The COMSOL comparison is limited to the elastic case, which the abstract admits; it does not test the PTZ scheme against FEM. And there's no code or data, which makes reproduction harder than it needs to be.\n\nOverall, the core of the paper—the scheme itself and its periodic stability analysis—is solid. The flaws are real but fixable, and they don't sink the central contribution. I'd send this to peer review, asking for the boundary-condition justification (or a softened claim) and the missing second-order derivation. For researchers in computational mechanics or viscoelasticity, this is worth reading; for others, it's a well-executed example of how thermodynamic structure can guide discretization.","headline":"A solid staggered finite-difference scheme for PTZ waves with rigorous periodic stability analysis, but the optimal C-hat=1 claim rests on an unproven boundary-condition assumption.","tokens_in":18658,"tokens_out":3158,"would_cite":false,"duration_ms":31088,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper introduces a staggered half-space, half-time finite-difference scheme for the Poynting–Thomson–Zener rheological solid and argues that at $\\alpha=1/2$, $\\hat C=1$ it is second-order accurate, stable, and free of dissipative and…","keywords":["Poynting–Thomson–Zener model","staggered grid","symplectic Euler","rheological waves","plane-wave stability analysis","numerical dissipation","numerical dispersion","finite difference scheme"],"falsifier":"For the actual finite sample (stress pulse at one end, free end at the other), assemble the full iteration matrix and compute its eigenvalues at $\\alpha=1/2$, $\\hat C=1$; any eigenvalue with modulus above 1, or a numerical experiment showing unbounded growth of total energy over many bounces, would falsify the claimed boundary-value stability.","tokens_in":17664,"feed_emoji":"🌊","tokens_out":23073,"duration_ms":187253,"temperature":0.7,"pith_summary":"The paper tries to establish that a finite-difference scheme built on staggering each field by half a space step and half a time step relative to the fields it couples to can solve the Poynting–Thomson–Zener (PTZ) rheological model both accurately and fast. Its central claim is that this placement, chosen from the spacetime roles of the quantities, turns the computationally identical symplectic Euler method into a second-order scheme and, at the interpolation parameter $\\alpha=1/2$ with the Courant condition based on the fastest wave speed $\\hat C\\le 1$, removes both dissipative and dispersive numerical error. The evidence is a plane-wave stability analysis, truncation-error calculations, and numerical experiments in which the scheme reproduces a stress-pulse shape with a few dozen spatial cells while a commercial finite-element package shows artificial damping, oscillations, or instability and needs run times 100 to 10,000 times longer. If the claim holds, the scheme provides a thermodynamics-motivated discretization recipe for continuum systems and shows that a widely used commercial finite-element package can fail qualitatively on the simplest elastic limit.","feed_headline":"One Courant setting frees rheological wave pulses from numerical error","feed_subtitle":"At that setting the scheme needs dozens of cells and runs 100–10,000 times faster than a commercial solver.","key_machinery":"The load-bearing mechanism is the staggered spacetime grid: velocities live at half-integer positions in both space and time relative to stress, strain shares stress's nodes, and temperature is half-shifted in time afterwards, so every discrete derivative is evaluated at the midpoint of the quantity it couples to. This arrangement makes the Hooke-case scheme computationally identical to symplectic Euler, the standard first-order geometric integrator, while raising its accuracy to second order by reflection symmetry, and it turns the PTZ constitutive equation into the explicit weighted update (33) with parameter $\\alpha$. Stability is controlled by the eigenvalues of the $3\\times 3$ iteration (transfer) matrix, which the paper tests through two classical polynomial root-location criteria; those criteria isolate the parameter-free thermodynamic condition $\\hat\\tau>\\tau$ and the fastest-wave Courant condition, and at $\\alpha=1/2$, $\\hat C=1$ the eigenvalues collapse to $1$ and $e^{\\pm ik\\Delta x}$, eliminating numerical dissipation and dispersion up to $O(\\Delta t/\\tau)$.","core_discovery":"The central discovery is that the PTZ wave system admits a staggered finite-difference realization in which stress, strain, and velocity sit at mutually half-shifted positions dictated by the equations containing them: velocity is half-shifted in both space and time from stress, strain sits with stress, and temperature is then half-shifted in time. The rheological equation, which couples each quantity to its own time derivative, is discretized by an $\\alpha$-weighted average that is explicit and second-order accurate at $\\alpha=1/2$. A plane-wave stability analysis of the resulting iteration matrix yields three conditions: the thermodynamic requirement $\\hat\\tau>\\tau$, a relation between $\\alpha$ and $\\Delta t$, and the Courant-type bound $\\hat C<1$ (extendable to $\\hat C\\le 1$ with boundary conditions), where $\\hat C=\\hat c\\,\\Delta t/\\Delta x$ uses the fast wave speed $\\hat c=\\sqrt{\\hat E/(\\tau\\varrho)}$. At the special point $\\alpha=1/2$, $\\hat C=1$, the three eigenvalues that multiply each Fourier mode per time step are exactly $1$ and $e^{\\pm ik\\Delta x}$, so all wavelengths travel at the same discrete speed and the dissipative and dispersive errors are of order $O(\\Delta t/\\tau)$. The same scheme applied to the elastic limit conserves total energy over many bounces and produces clean wave pulses where the commercial finite-element software COMSOL, with several tuned time-stepping methods, gives damped, oscillatory, or unstable results and takes 100 to 10,000 times longer.","pith_inferences":["A natural extension, listed by the authors as future work but not demonstrated, is that the same half-shift recipe transfers to other members of the Kluitenberg–Verhás family and to non-Fourier heat conduction, with parabolic limits likely requiring mixed Courant conditions analogous to the Kelvin–Voigt case.","Because the elastic limit coincides with symplectic Euler, a plausible conjecture the paper does not prove is that the full PTZ scheme inherits a discrete variational or symplectic structure, which would explain the observed total-energy conservation.","A direct testable extension would be to compare the scheme's temperature histories with an analytic PTZ solution in the force-equilibrial limit; the paper leaves analytic comparison as future work and does not yet provide a convergence study of the thermal field.","The COMSOL comparison suggests a broader benchmarking lesson: commercial finite-element packages may be unreliable for viscoelastic wave propagation unless the time stepper and tolerances are carefully chosen, so benchmark suites should include dissipative wave problems rather than only static or quasi-static cases."],"forward_implications":["Setting $\\alpha=1/2$ and $\\hat C\\le 1$ is enough for stability of the PTZ scheme, and with $\\hat C=1$ the discrete dispersion branches are linear, so the scheme has no numerical dissipation or dispersion up to $O(\\Delta t/\\tau)$.","In the Hooke limit the scheme reduces to symplectic Euler, is stable for $C\\le 1$, conserves elastic plus kinetic energy over long times, and needs $C=1$ to avoid both dissipative and dispersive artifacts.","For the Kelvin–Voigt limit ($\\tau=0$) the stability conditions become $\\alpha<1/2$ (relaxed to $\\alpha\\le 1/2$ with boundary conditions) together with a mixed parabolic-hyperbolic Courant bound combining $\\Delta t^2$ and $\\hat\\tau\\Delta t$ terms.","The stability analysis is not purely numerical: it reproduces the thermodynamic stability condition $\\hat\\tau>\\tau$ as a scheme-independent requirement, so numerical stability criteria can teach something about the underlying continuum model.","Practically, the $\\alpha=1/2$ scheme gives a reliable PTZ stress-signal shape with as few as 25 to 50 spatial cells, whereas $\\alpha=0$ needs more than 1000 cells for comparable quality."],"supporting_citations":[{"why":"It supplies the staggered half-space placement for heat conduction that the present scheme generalizes by adding half-time shifts.","marker":"[1]"},{"why":"It derives the PTZ model and the inequality $\\hat\\tau>\\tau$ from a thermodynamic internal-variable framework, providing the continuum model under study.","marker":"[18]"},{"why":"It documents the symplectic Euler method, with which the Hooke-case scheme coincides, and its first-order accuracy and long-time behavior.","marker":"[22]"},{"why":"It provides the plane-wave stability analysis technique applied to the transfer matrix of the scheme.","marker":"[23]"},{"why":"It supplies the theorem on algebraic versus geometric multiplicity used to judge neutral stability when $|\\xi|=1$.","marker":"[24]"},{"why":"It backs the stability criterion for repeated eigenvalues, supporting the extension of the stability region to closed boundaries.","marker":"[25]"},{"why":"It gives a staggered-field electrodynamics method with similar stability expectations, cited as a parallel case for staggered discretization.","marker":"[26]"},{"why":"It provides the innerwise-determinant criteria used to reduce the PTZ stability conditions to explicit inequalities.","marker":"[27]"}],"fun_headline_variants":["Exact Courant point kills error in rheological wave scheme","Staggered half-shifts make rheological waves error-free","One Courant equality makes rheological waves exact","Half-shift scheme: elastic waves stay clean, COMSOL fails","Exact transmission at one Courant point beats COMSOL"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that stability proved for waves on an infinitely long periodic medium also holds for the finite sample with a stress pulse at one end and a free end, a boundary-value extension the paper treats as a rule of thumb rather than a proof.","fun_headline_variants_meta":{"raw":{"variants":["Exact Courant point kills error in rheological wave scheme","Staggered half-shifts make rheological waves error-free","One Courant equality makes rheological waves exact","Half-shift scheme: elastic waves stay clean, COMSOL fails","Exact transmission at one Courant point beats COMSOL"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00052,"raw_usage":{"total_tokens":2565,"prompt_tokens":1039,"completion_tokens":1526,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":655,"completion_tokens_details":{"reasoning_tokens":1441}},"tokens_in":655,"tokens_out":1526,"duration_ms":11222,"temperature":1.0,"reasoning_tokens":1441,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T11:52:01.326121+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"For the actual finite sample (stress pulse at one end, free end at the other), assemble the full iteration matrix and compute its eigenvalues at $\\alpha=1/2$, $\\hat C=1$; any eigenvalue with modulus above 1, or a numerical experiment showing unbounded growth of total energy over many bounces, would falsify the claimed boundary-value stability.","supporting_citations":[{"cited_title":"Implicit numerical schemes for generalized heat conduction equations.International Journal of Heat and Mass T ransfer 2018, 126, 1177–1182","cited_arxiv_id":null,"evidence_quote":"It supplies the staggered half-space placement for heat conduction that the present scheme generalizes by adding half-time shifts."},{"cited_title":"Distinguished rheological models for solids in the framework of a thermodynamical internal variable theory","cited_arxiv_id":null,"evidence_quote":"It derives the PTZ model and the inequality $\\hat\\tau>\\tau$ from a thermodynamic internal-variable framework, providing the continuum model under study."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"It documents the symplectic Euler method, with which the Hooke-case scheme coincides, and its first-order accuracy and long-time behavior."},{"cited_title":"Numerical integration of the barotropic vorticity equation","cited_arxiv_id":null,"evidence_quote":"It provides the plane-wave stability analysis technique applied to the transfer matrix of the scheme."},{"cited_title":"An Introduction to Difference Equations, 3rd ed.; Springer, New York, USA, 2005","cited_arxiv_id":null,"evidence_quote":"It supplies the theorem on algebraic versus geometric multiplicity used to judge neutral stability when $|\\xi|=1$."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"It backs the stability criterion for repeated eigenvalues, supporting the extension of the stability region to closed boundaries."},{"cited_title":"Von Neumann stability analysis of globally divergence-free RKDG schemes for the induction equation using multidimensional Riemann solvers","cited_arxiv_id":null,"evidence_quote":"It gives a staggered-field electrodynamics method with similar stability expectations, cited as a parallel case for staggered discretization."},{"cited_title":"Inners and Stability of Dynamical Systems , John Wiley & Sons: New York, USA, 1974","cited_arxiv_id":null,"evidence_quote":"It provides the innerwise-determinant criteria used to reduce the PTZ stability conditions to explicit inequalities."}],"review_version":1}