{"id":"29f36dee-f57f-4537-ae86-071cc1b594c7","arxiv_id":"2607.28309","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"Large-time-step marching-on-in-time volume integral equation analysis for dielectrics is demonstrated, with an FFT-based matrix-free present-time operator reaching 15.6 million unknowns in single-threaded execution.","lead":"A time-domain volume integral equation solver for dielectric scattering is shown to run with time steps up to 16 times larger than the Courant reference for band-limited pulses, provided the Nyquist sampling condition is met. To handle the denser present-time interaction matrix at large time steps, the paper applies an FFT-based matrix-free multiplication, enabling a 15.6 million-unknown simulation on a single CPU thread.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Large-time-step accuracy claim is validated only at single observation points; the promised global relative ℓ2-norm error never appears, so the 16× accuracy claim is not yet established.","rationale":"The reader's weakest assumption correctly identifies the load-bearing premise: the large-time-step advantage depends on temporal sampling and on the validity of pointwise accuracy checks. My stress-test agrees and sharpens the concern: the paper promises a global relative ℓ2-norm error in Section V-A but never delivers it, and the existing comparisons are single-point. The benchmark in Section III-C uses a single current component at one point with an absolute error that may be a large fraction of the local signal; the homogeneous-cube check uses the same solver's α=1 solution as reference, so it cannot detect systematic errors introduced by the large time step itself. The FFT-based MVP implementation is separately verified at α=1 with 1e-4 pointwise error, but that verifies the matrix-free operator, not the accuracy of α=16 time stepping. My recommendation is not to reject the work—the computational-shift analysis and FFT acceleration are plausible and internally consistent—but to keep the verdict CONDITIONAL until a global accuracy metric is supplied. Hence UNCHANGED relative to the reader's verdict.","tokens_in":16314,"tokens_out":8443,"duration_ms":86468,"concrete_test":"Reproduce the homogeneous-cube simulation (fine 80³ voxelization, α=1 and α=16) and compute the global temporal relative ℓ2-norm error E = sqrt( Σ_t ||J_α(t) − J_1(t)||²_2 / Σ_t ||J_1(t)||²_2 ) over all 80³×3 unknowns. Also compute the same global error for the benchmark cube of Section III-C against the MOD-JVIE reference over the full domain, not just at the single observation point (0.025,0.075,0.025). If either error at α=16 exceeds a stated threshold (e.g., 5%), the 'accurate up to 16×' claim is not supported; if both are small, the concern is resolved.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim—accurate MOT-JVIE at Δt=16Δt_CFL—rests on the assumption that the temporal sampling condition (5) with κ=2 is sufficient for the full transient solution. Equation (5) bounds only the incident-field bandwidth; the actual current-density error is governed by the interaction between the quadratic-spline temporal basis and the retarded Green function, which is not guaranteed by Nyquist sampling alone. The validation is pointwise rather than global: Fig. 3 shows one Jx voxel against the α=1 self-reference; Fig. 6 shows one azimuthal current component against the MOD-JVIE reference, with only an absolute-error curve whose peak (≈0.06 A/m²) is comparable to the signal amplitude (≈0.2 A/m²); Fig. 10 checks one voxel at α=1 only. Section V-A explicitly promises 'relative ℓ2-norm errors,' but no such global metric is reported. Moreover, Section V states that the underlying direct MOT-JVIE was validated previously at small time steps, so the large-time-step regime has no independent ground truth. If the global relative error at α=16 is large—especially at late times or near material discontinuities—the 'up to 16×' assertion overstates what the evidence supports.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript studies large-time-step operation of a marching-on-in-time contrast-current-density volume integral equation (MOT-JVIE) for transient dielectric scattering. The time step is chosen by a Nyquist sampling condition (Eq. (5)) rather than by the CFL limit of the spatial voxel grid, and the paper demonstrates, for band-limited excitations, accurate solutions for time steps up to 16 times the reference CFL-limited step. A central computational observation is that as the time step grows, the present-time interaction matrix Z0 densifies, so the dominant cost shifts from history-term evaluations to repeated Z0 matrix-vector products inside the iterative solver. To address this, the authors apply a matrix-free FFT-based evaluation of Z0 exploiting its multilevel Toeplitz structure, and demonstrate more than an order-of-magnitude speedup on a cubic test case plus a single-threaded simulation of an 8x8 bilayer nanopillar metasurface with about 15.6 million unknowns.","tokens_in":16459,"tokens_out":3223,"duration_ms":30326,"significance":"If the accuracy and efficiency claims hold, the paper would be a substantial practical advance: it identifies a regime change in MOT-JVIE solvers, provides a concrete FFT-based remedy, and gives a credible large-scale demonstration without parallel hardware. The manuscript has notable strengths: the mathematical core (Toeplitz/circulant FFT MVP) is standard and the pointwise comparison against the direct sparse implementation at the 1e-4 level is a useful consistency check; there are no fitted parameters; and the use of an independent MOD-JVIE reference for the benchmark cube provides some external grounding. The main weakness is that the central accuracy claim is supported only by pointwise or single-observation-point comparisons, despite the paper explicitly promising a global relative l2-norm error in Section V-A.","major_comments":[{"comment":"The text at the start of Section V-A states that the comparison is performed 'in terms of induced contrast current density, relative ℓ2-norm error, runtime, and TFQMR convergence,' but no global relative ℓ2-norm error is ever reported. The only accuracy plot is Fig. 10(b), a normalized pointwise error at a single voxel and only for α=1. Since the direct implementation is available for the cube, please add the global relative ℓ2-norm error of the current-density vector over all voxels and all time steps for α=2, 4, 8, and 16, and indicate where in the domain and at which times the error is largest.","section":"Section V-A and Fig. 10"},{"comment":"The large-time-step accuracy claim is validated by one current-density component at voxel (41,41,41) for the homogeneous cube and by one azimuthal component at one observation point for the benchmark cube. These single-point checks do not by themselves support the abstract's global statement that 'accurate transient analysis is demonstrated for time steps up to 16 times larger.' Please provide a multi-voxel or global error assessment, ideally a convergence table in α for both the coarse and fine voxelizations, and explicitly check late-time behavior and regions near material discontinuities.","section":"Section III-B and Fig. 3"},{"comment":"Equation (5) is a Nyquist criterion on the incident-field bandwidth; it does not by itself guarantee that the full MOT-JVIE solution is temporally accurate, because the temporal discretization also involves the quadratic-spline basis and the retarded Green-function interaction. The sufficiency of κ=2 should be demonstrated by a time-step convergence study, not inferred from the sampling formula alone. The observed degradation at α=32 in Fig. 3(b) is consistent with the Nyquist bound, but it does not quantify how close α=16 is to the accuracy threshold; a global error-versus-α curve would settle this.","section":"Section III-B, Eq. (5)"},{"comment":"The 8x8 nanopillar-array demonstration includes no quantitative accuracy check: Fig. 12 illustrates a spatial current pattern, but there is no comparison against a reference solution or a converged solution on a coarser or smaller problem. Because the array result is the flagship scalability claim, please add a validation step (for example, a smaller 2x2 or 4x4 array computed with the direct sparse implementation, or a coarse-grid comparison) to show that the FFT-based large-time-step solver retains the accuracy established for the cube.","section":"Section V-B"}],"minor_comments":[{"comment":"The title and running header contain a typo: 'V olume' should be 'Volume.'","section":"Title/header"},{"comment":"The manuscript uses 'MOT-JVIE,' 'MOT-TDJVIE,' and 'TDJVIE' interchangeably; please standardize the terminology.","section":"Throughout"},{"comment":"Figure 4 would benefit from an explicit statement that the vertical axis is on a logarithmic scale and from a caption defining what 'LHS' and 'RHS' mean for readers who skip Eq. (6).","section":"Fig. 4"},{"comment":"The sentence referencing 'Eq. (35) in [25]' and the definition of the δ-regularized formulation is not self-contained; a brief explanation of the regularization parameter δn would improve readability.","section":"Section II"},{"comment":"In Fig. 10(b), the normalized pointwise error is normalized by the temporal maximum of the direct solution at that voxel; stating this in the main text (rather than only in the figure) would avoid ambiguity about what 'pointwise' means.","section":"Section V-A"}],"recommendation":"major_revision","confidential_remarks":"The paper is built directly on the authors' prior work [23] and [25], which is legitimate, but the novelty relative to those papers should be made crisper: the new contribution is the identification of the Z0 bottleneck and the FFT-based Z0 MVP, not the MOT-JVIE formulation itself. The missing global ℓ2 error is a concrete, fixable issue, but it is load-bearing for the advertised 'up to 16x' accuracy claim, so I recommend major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things you should know before reading: the paper is a competent incremental advance in transient volume integral equations, and the headline claim of accurate operation at 16× the CFL time step is less well supported than the abstract suggests.\n\nWhat is genuinely new is the bottleneck analysis. The authors show clearly, with sparsity tables and runtime breakdowns, that as Δt grows, the present-time interaction matrix Z0 densifies and the dominant cost shifts from history-term accumulation to repeated matrix–vector products with Z0. The FFT-based multilevel Toeplitz acceleration is a standard trick from frequency-domain integral equations, but applying it to the present-time matrix rather than only to history terms as in their prior work is a sensible and effective adaptation. The implementation check is credible: a pointwise relative error around 1e-4 between the direct sparse MVP and the FFT-based MVP, and the 15.6-million-unknown single-threaded run is a real data point even if the 'first beyond 15 million' claim is unverifiable and not worth much. The citations to their own earlier work are legitimate; they are building on published formulation and acceleration methods.\n\nThe soft spot is the validation of the large-time-step accuracy claim. Section V-A explicitly promises 'relative ℓ2-norm errors,' but no such global metric is reported. Fig. 10 compares direct versus FFT only at α=1 and only at one voxel. The benchmark against MOD-JVIE in Fig. 6 is also single-point, and the peak absolute error around 0.06 A/m² is not negligible against the signal amplitude around 0.2 A/m². The paper leans on the Nyquist sampling condition (5) with κ=2 to justify the large steps, but that condition only constrains the incident-field bandwidth; it does not by itself guarantee global accuracy of the MOT solution, particularly near material discontinuities or at late times. The paper honestly shows α=32 deviating, so the operating point is bandwidth-limited rather than a stability margin, but then 'up to 16×' should be framed as a demonstrated operating point, not a general capability.\n\nNone of this invalidates the method. The math is standard, the implementation appears sound, and the bottleneck-shift observation is a useful contribution to the MOT-TDVIE literature. What the paper needs before it is a definitive demonstration is straightforward: report global relative ℓ2-norm errors for the inhomogeneous cube at α=8 and α=16, show the direct-versus-FFT comparison at a large time step, and temper the 'first beyond 15 million' phrasing.\n\nFor peer review: yes, send it. The core idea is sound and the missing items are fixable in revision. I would accept after those additions.","headline":"Competent, useful incremental advance in large-time-step MOT-JVIE, but the 16× accuracy claim is under-supported by pointwise checks and a missing global error norm.","tokens_in":17059,"tokens_out":4972,"would_cite":true,"duration_ms":41632,"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":"A marching-on-in-time volume integral solver for dielectric scattering takes time steps up to 16 times the CFL reference and shifts its dominant cost to an FFT-accelerated present-time matrix.","keywords":["marching-on-in-time","volume integral equation","dielectric scattering","CFL condition","large time step","matrix-free FFT","Toeplitz structure","metasurface"],"falsifier":"Run the 200 nm homogeneous cube with $\\varepsilon_r=12$ at $\\Delta t = 16\\,\\Delta t_{\\mathrm{CFL}}$ and compute the global relative $\\ell^2$ error between the $\\alpha=16$ and $\\alpha=1$ solutions over all voxels and all time steps; if that global error is large even though single-point traces agree, the claim that 16$\\times$ steps are accurate for band-limited excitations is falsified.","tokens_in":16025,"feed_emoji":"⚡","tokens_out":6993,"duration_ms":53392,"temperature":0.7,"pith_summary":"The paper argues that a marching-on-in-time volume integral equation for dielectric scattering can be run with time steps up to 16 times larger than the CFL-limited step tied to the voxel grid, provided the excitation is band-limited enough to satisfy the Nyquist sampling condition $\\Delta t \\le 1/(2\\kappa f_{\\max})$ with oversampling factor $\\kappa=2$. In this large-time-step regime the computational bottleneck is not the accumulated history terms but the present-time interaction matrix $\\mathbf{Z}_0$, whose nonzero population grows with the time step. The paper makes $\\mathbf{Z}_0$ matrix-free by applying its Green-function part through multidimensional FFTs, exploiting the multilevel Toeplitz structure of the volume integral operator. On an inhomogeneous dielectric cube the accelerated solver reproduces the direct solution's current density and iteration counts while cutting runtime by more than an order of magnitude, and on an $8\\times8$ bilayer nanopillar metasurface it handles about 15.6 million unknowns on a single CPU thread. The result is that large-time-step MOT-JVIE analysis is computationally viable and the field's acceleration target shifts from history updates to the present-time operator.","feed_headline":"Time-domain scattering solver steps 16x past the CFL limit","feed_subtitle":"An FFT-based matrix-free trick cuts the main cost, enabling 15.6M-unknown metasurface runs on one CPU thread.","key_machinery":"The load-bearing object is the present-time interaction matrix $\\mathbf{Z}_0$ and its decomposition $\\mathbf{Z}_0 = \\mathbf{D} - \\mathbf{X}\\mathbf{G}$, where $\\mathbf{D}$ and $\\mathbf{X}$ are diagonal matrices for local constitutive and material-contrast terms and $\\mathbf{G}$ is the volume-integral Green-function operator. Because the voxel grid is uniform and the background is homogeneous, the entries of $\\mathbf{G}$ depend only on relative source--observer voxel offsets, giving it a multilevel block Toeplitz structure; embedding it in a circulant operator lets each matrix--vector product be computed by zero-padding, forward three-dimensional FFT, spectral multiplication with precomputed kernels, inverse FFT, and truncation back to the physical grid. This matrix-free application of $\\mathbf{Z}_0$ is what removes the large-time-step bottleneck. The temporal admissibility condition $\\Delta t \\le 1/(2\\kappa f_{\\max})$ with $\\kappa=2$ is what licenses the 16$\\times$ step in the first place, together with quadratic-spline temporal basis functions and TFQMR iteration for the implicit marching system.","core_discovery":"The central discovery is that in a marching-on-in-time current-density volume integral equation (MOT-JVIE) solver, the time step is limited by the Nyquist sampling of the band-limited incident and scattered fields rather than by a CFL stability condition, and that choosing $\\Delta t = 16\\,\\Delta t_{\\mathrm{CFL}}$ preserves accurate transient responses for the pulses considered. As the time step grows, the spherical causal wavefront that fills the present-time matrix in one step encloses more of the voxelized scatterer, so $\\mathbf{Z}_0$ densifies rapidly. The paper shows the cost of the iterative linear solve then dominates, while TFQMR iteration counts stay low, pinpointing the matrix--vector product with $\\mathbf{Z}_0$ as the true cost driver. That product is evaluated without forming $\\mathbf{Z}_0$ by writing $\\mathbf{Z}_0 = \\mathbf{D} - \\mathbf{X}\\mathbf{G}$ and applying $\\mathbf{G}$ as a multilevel Toeplitz/circulant convolution through three-dimensional FFTs, cutting the operator application to $\\mathcal{O}(N\\log N)$ and removing the memory ceiling of explicit sparse storage. Evidence includes matching current-density traces, normalized pointwise errors near $10^{-4}$, up to 15$\\times$ speedups, and a full $8\\times8$ GST/a-Si nanopillar array simulation with roughly $1.56\\times10^{7}$ vector unknowns in single-threaded execution.","pith_inferences":["The same bandwidth-limited large-time-step picture should carry over to other volume integral formulations on uniform grids — magnetic, dispersive, or nonlinear materials — whenever the Green-function operator retains Toeplitz structure; the paper's scalar nondispersive permittivity is a simplifying choice, not the mechanism.","The accuracy ceiling is the Nyquist condition with $\\kappa=2$: the paper's own $\\alpha=32$ result shows the method is bandwidth-limited rather than unconditionally stable, so a user must know the effective $f_{\\max}$ of the excitation before choosing $\\Delta t$.","FFT-based $\\mathbf{Z}_0$ application should combine naturally with plane-wave time-domain or hierarchical FFT history acceleration, potentially also speeding up moderate-time-step regimes where neither history terms nor $\\mathbf{Z}_0$ alone dominate.","A global relative $\\ell^2$ error over all voxels and all time steps, rather than pointwise traces, would settle how much conservatism is hidden in the $\\kappa=2$ oversampling choice."],"forward_implications":["Band-limited transient dielectric scattering can be marched with up to 16 times fewer time steps than the voxel-CFL reference, directly reducing time-marching, storage, and history-accumulation cost.","In the large-time-step regime, further acceleration should target the present-time operator $\\mathbf{Z}_0$ rather than history terms; the FFT-based Toeplitz/circulant application scales as $\\mathcal{O}(N\\log N)$ and never forms $\\mathbf{Z}_0$ explicitly.","The matrix-free present-time product lifts the memory wall that makes direct MOT-JVIE impractical for large metasurface arrays, demonstrated at about 15.6 million vector unknowns on a single CPU thread.","Because TFQMR iteration counts remain low and stable as the time step grows, the FFT speedup is not eroded by slower convergence in the tested cases."],"supporting_citations":[{"why":"Supplies the MOT-JVIE space-time discretization, quadratic-spline temporal basis, and the alpha=1 reference solution used throughout.","marker":"[23]"},{"why":"Provides the FFT-based delayed-history acceleration and the delta-regularization notation; the present work reuses this framework for the RHS and motivates the Z0 bottleneck.","marker":"[25]"},{"why":"Supplies the MOD-JVIE benchmark solution against which the large-time-step MOT-JVIE responses are validated.","marker":"[31]"},{"why":"Defines TFQMR, the iterative solver whose per-iteration Z0 matrix-vector products become the dominant cost.","marker":"[24]"},{"why":"Supplies the three-dimensional weak-form conjugate-gradient FFT method for scattering, the established Toeplitz/circulant FFT treatment applied to G.","marker":"[35]"},{"why":"Establishes the CFL condition for explicit differential-equation solvers, used only as the reference Delta-t-CFL scale.","marker":"[4]"}],"fun_headline_variants":["Dielectric solver beats CFL by 16x with FFT trick","MOT-JVIE scales to 15.6M unknowns on one thread","Breaking CFL: Volume integral solver takes 16x time steps","Matrix-free FFT accelerates dielectric scattering beyond CFL","Large time steps cut dielectric scattering cost 15x"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the incident and scattered fields are effectively band-limited and that the oversampling factor $\\kappa=2$ in the sampling condition is sufficient to guarantee the claimed transient accuracy; if that sampling assumption fails, the 16$\\times$ time-step claim no longer applies.","fun_headline_variants_meta":{"raw":{"variants":["Dielectric solver beats CFL by 16x with FFT trick","MOT-JVIE scales to 15.6M unknowns on one thread","Breaking CFL: Volume integral solver takes 16x time steps","Matrix-free FFT accelerates dielectric scattering beyond CFL","Large time steps cut dielectric scattering cost 15x"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001126,"raw_usage":{"total_tokens":4784,"prompt_tokens":1148,"completion_tokens":3636,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":764,"completion_tokens_details":{"reasoning_tokens":3547}},"tokens_in":764,"tokens_out":3636,"duration_ms":24180,"temperature":1.0,"reasoning_tokens":3547,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T15:22:07.426089+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the 200 nm homogeneous cube with $\\varepsilon_r=12$ at $\\Delta t = 16\\,\\Delta t_{\\mathrm{CFL}}$ and compute the global relative $\\ell^2$ error between the $\\alpha=16$ and $\\alpha=1$ solutions over all voxels and all time steps; if that global error is large even though single-point traces agree, the claim that 16$\\times$ steps are accurate for band-limited excitations is falsified.","supporting_citations":[{"cited_title":"The influence of contrast and temporal expansion on the marching-on-in-t ime contrast current density volume integral equation,","cited_arxiv_id":null,"evidence_quote":"Supplies the MOT-JVIE space-time discretization, quadratic-spline temporal basis, and the alpha=1 reference solution used throughout."},{"cited_title":"FFT - acceleration and stabilization of the 3D marching-on-in-t ime contrast current density volume integral equation for scattering fr om high contrast dielectrics,","cited_arxiv_id":null,"evidence_quote":"Provides the FFT-based delayed-history acceleration and the delta-regularization notation; the present work reuses this framework for the RHS and motivates the Z0 bottleneck."},{"cited_title":"A time-domain volume integral equ ation and its marching-on-in-degree solution for analysis of disper sive dielectric objects,","cited_arxiv_id":null,"evidence_quote":"Supplies the MOD-JVIE benchmark solution against which the large-time-step MOT-JVIE responses are validated."},{"cited_title":"A transpose-free quasi-minimal residua l algorithm for non-hermitian linear systems,","cited_arxiv_id":null,"evidence_quote":"Defines TFQMR, the iterative solver whose per-iteration Z0 matrix-vector products become the dominant cost."},{"cited_title":"The three-dimensional weak form of the conjugate gradient FFT method for solving scattering problems,","cited_arxiv_id":null,"evidence_quote":"Supplies the three-dimensional weak-form conjugate-gradient FFT method for scattering, the established Toeplitz/circulant FFT treatment applied to G."},{"cited_title":"Taflove, Computational Electrodynamics: The Finite-Difference Time- Domain Method","cited_arxiv_id":null,"evidence_quote":"Establishes the CFL condition for explicit differential-equation solvers, used only as the reference Delta-t-CFL scale."}],"review_version":2}