{"id":"5a4fb0d5-b141-4da6-a56e-fd9e3a4375e9","arxiv_id":"2607.18103","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"A time-domain harmonic balance method with a spatiotemporal splitting solves radio-frequency capacitively coupled plasma fluid simulations about 10x faster than conventional time stepping, with <0.3% error on a 1D argon benchmark.","lead":"This paper brings a method called harmonic balance to the simulation of radio-frequency plasmas, letting computers skip the long warm-up phase and jump directly to the steady oscillating state. It runs more than ten times faster than standard time-marching simulations on a test case, with errors below 0.3% for the main plasma quantities.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"NH=8 is selected from the DTS reference against which accuracy is measured; the '<0.3%' claim is an in-sample fit and is contradicted by Table 4's E∞ values.","rationale":"We focused on the strongest pillar of the central claim: the sufficiency of NH=8. Because the harmonic count is chosen from the same reference solution used for validation, the reported accuracy is not an independent test of the method's general validity. This is a circularity in the validation protocol, not an error in the derivation. The reader's stated weakest assumption (T-periodicity) is less compelling: for the modeled dissipative system with a periodic drive, the asymptotic state is expected to be T-periodic, and HB is designed to compute it directly; slow species affect transient duration, not the existence of the periodic orbit. The abstract's <0.3% wording is also imprecise (Table 4's E∞ values exceed it), but the authors do disclose the Te caveat in §3.5. The proposed fixed-NH re-run on a new operating point would settle whether the method's accuracy claim generalizes. Since this supports the reader's CONDITIONAL verdict rather than overturning it, we recommend UNCHANGED.","tokens_in":18538,"tokens_out":11205,"duration_ms":117297,"concrete_test":"Run the HB solver with NH=8 fixed (no spectral preprocessing) on a second 1D argon CCP case, e.g., p=0.5 Torr or V0=150 V, keeping all other algorithms identical. Generate a fresh fully converged DTS baseline using the paper's own criteria (T/Δt=200, pseudo-CFL=10000, ninner=300). Compute the Table 4 metrics E2, EMA, E∞ for ne, εe, Te, ϕ. If ne/εe/ϕ E2 remain below 0.3% and Te E2 below ~1%, the circularity concern is answered. Also report the HB residual norm at convergence and the amplitude of harmonic NH+1 from the HB solution, to show the spectral tail is negligible without using DTS.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's headline accuracy claim depends on choosing NH=8. That choice comes from an FFT of the DTS baseline itself (§3.3: 'This preliminary analysis relies on the high-fidelity transient signals extracted directly from the reference DTS baseline'). The HB solution with NH=8 is then compared to that same baseline (§3.5, Table 4), so the reported <0.3% error is in-sample: it demonstrates that a truncation tailored to the reference can represent the reference, not that the HB method reliably predicts the periodic state when the harmonic content is unknown. The convergence sweep NH=6,8,10,12 does not resolve this: the tested range is derived from the reference spectrum, and no error-vs-NH table is given, so one cannot tell whether NH=8 is a true plateau or just an interpolated sweet spot. As a result, the broader conclusion that 'retaining eight harmonics perfectly resolves' RF CCP dynamics is not established for operating points outside this tuned 1D benchmark. Separately, the abstract's 'strictly below 0.3%' is not supported even for ne: Table 4 gives E∞(ne)=0.392%, and Te has E2=1.18%, E∞=6.385%. The periodicity assumption (Eq. 15) is standard and less concerning for this dissipative periodically forced model.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents a time-domain harmonic balance (HB) method for the periodic steady state of a fluid model of radio-frequency capacitively coupled plasmas, specifically a drift-diffusion-Poisson system with electron-energy transport (LMEA). The physical time derivative is replaced by a dense spectral operator at NT = 2NH + 1 temporal collocation points, and the resulting quasi-steady system is solved by an implicit pseudo-time relaxation that combines a spatial block-implicit sweep (decoupled across phases) with a cell-local dense temporal inversion; a semi-implicit Poisson update handles the dielectric-relaxation limit. The method is validated against a 1D argon CCP benchmark, using a dual-time-stepping (DTS) baseline with T/Δt = 200 as reference. The authors select NH = 8 from an FFT of that baseline and report E2 errors below 0.3% for ne, εe, and ϕ, and a 10.26x speedup over the DTS baseline, with the conclusion that eight harmonics 'perfectly resolves' the discharge dynamics.","tokens_in":18942,"tokens_out":7194,"duration_ms":82370,"significance":"If the central claims hold, this is a meaningful contribution: prior HB plasma work (ref. [25]) was limited to the local-field approximation, whereas this paper includes full electron-energy transport and proposes a memory-efficient, cell-local temporal inversion that avoids global Jacobian assembly. The paper contains a careful time-step refinement study showing second-order convergence of the DTS reference and a detailed spatiotemporal comparison against a classical benchmark. However, the headline accuracy claim is weakened by the in-sample selection of NH from the DTS reference, the absence of an error-versus-NH table, and the restriction of all validation to a single 1D operating point. The speedup result, while plausible, is also only demonstrated in 1D on a single processor.","major_comments":[{"comment":"The harmonic count NH = 8 is chosen from an FFT of the DTS baseline (§3.3: 'This preliminary analysis relies on the high-fidelity transient signals extracted directly from the reference DTS baseline'), and the error metrics in Table 4 compare the NH = 8 HB solution to that same baseline. The reported <0.3% errors are therefore in-sample: they show that a truncation tailored to the reference can represent that reference, not that the method reliably predicts the periodic state when the harmonic content is unknown. The sweep NH = 6, 8, 10, 12 in §3.4 is not accompanied by an error table; only convergence histories and selected profiles are shown, so the reader cannot tell whether NH = 8 is a true plateau or a tuned sweet spot. Please provide E2/EMA/E∞ for all four truncation levels for all variables, and ideally an out-of-sample test (e.g., a different V0 or pressure) to support the genera","section":"§3.3, §3.4, and Table 4"},{"comment":"The abstract and conclusion state that 'macroscopic relative errors [are] strictly below 0.3%' compared to DTS solutions. Table 4 contradicts this as stated: E∞(ne) = 0.392%, E2(Te) = 1.18%, and E∞(Te) = 6.385%. The body text carefully qualifies the 0.3% statement to the L2 and mean-absolute errors of the directly integrated variables (ne, εe, ϕ), but the abstract and conclusion are unqualified. Revise the headline claim to match the data, e.g., 'L2 and mean-absolute errors below 0.3% for directly integrated variables, with larger localized pointwise errors for derived quantities such as Te near the walls.'","section":"Abstract and §3.5, Table 4"},{"comment":"The entire validation is performed in one dimension on a 90-cell mesh. The abstract promises a 'physically rigorous, memory-efficient, and highly accelerated paradigm for practical RF plasma simulations', but no multidimensional test is presented. The proposed cell-local temporal inversion is dimension-independent in principle, but solver robustness, convergence behavior, and the claimed speedup have not been demonstrated in 2D or 3D. Please either add a multidimensional proof-of-concept or explicitly restrict the conclusions and title claims to 1D benchmarks.","section":"§3.1, §3.5, Fig. 14"},{"comment":"The DTS reference itself carries finite temporal discretization error: Table 3 shows about 0.024% error in spatial averages for T/Δt = 200 versus T/Δt = 400, and local pointwise errors could be larger. The HB errors in Table 4 are measured relative to this reference. In particular, E∞(Te) = 6.385% is attributed to the sensitivity of the ratio εe/ne in depleted sheaths, but that attribution is not verified; part of this deviation could be temporal error of the DTS baseline. To separate HB truncation/aliasing error from DTS discretization error, compare the HB solution against a T/Δt = 400 (or Richardson-extrapolated) DTS reference, at least for Te.","section":"§3.5, Table 4, Eq. (71)"}],"minor_comments":[{"comment":"The notation in Eq. (66) is ambiguous: 'ne,me' appears to mean n_{e,m} times the elementary charge e, but should be written explicitly as such. Also, the dimensionless form of the equations is said to be omitted; providing it, or at least the values of cD, α, and the entries of D_damp, would improve reproducibility.","section":"Eq. (66) and §3.1"},{"comment":"The numerical diffusion constant c_D, the damping factor α, and the diagonal damping matrix D_damp are user-set parameters but their values are never specified. The paper should list these values and briefly discuss sensitivity to them, since the convergence and stability claims depend on them.","section":"§2.4, Eq. (43)–(44), §2.4, Eq. (49)"},{"comment":"The error definition in Eq. (71) samples Ns = 200 phase instances per RF cycle, while the HB solution has only NT = 17 collocation points for NH = 8. The text should state explicitly how the HB solution is evaluated at those 200 samples (e.g., inverse DFT reconstruction from the Fourier coefficients) so that the error metric is unambiguous.","section":"§3.5, Eq. (71)"},{"comment":"Reference [25] is cited as 'Journal of Computational Physics (2026) 115027' without volume/page details; add the full citation or DOI if available. Also, the caption of Fig. 3 should identify the pseudo-CFL values by line style or color, as the reader cannot distinguish them from the text description alone.","section":"References"}],"recommendation":"major_revision","confidential_remarks":"The novelty claim relative to [25] should be verified by the editor or a separate check; the manuscript does not compare directly with [25] or discuss its algorithmic differences in detail. The in-sample harmonic selection is the main correctness risk; a blind or out-of-sample test would substantially increase confidence. The 1D-only scope should be acknowledged in the abstract and conclusions even after revision."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The short version: solid incremental work, probably correct, but the headline accuracy claim is in-sample and overstated. The paper is worth a careful review, not a desk reject.\n\nWhat is actually new: previous harmonic balance work on CCPs [25] used the local-field approximation and skipped the electron-energy equation. This paper extends the method to a fully coupled drift-diffusion-Poisson system with electron-energy transport (LMEA). That is a real step, because the energy coupling adds stiffness. The proposed spatiotemporal operator splitting—spatial implicit relaxation plus cell-local temporal inversion—is a sensible way to avoid global Jacobians, and the implementation appears careful. The authors verify their DTS baseline with a time-step refinement study showing clean second-order convergence, which is more than many papers do. For the 1D benchmark, the HB solution with NH=8 gives relative L2 errors around 0.2-0.3% for ne, εe, and ϕ, and delivers a 10x speedup over the fully converged DTS baseline. Those numbers are plausible.\n\nSoft spots, in order of importance. First, the abstract's \"strictly below 0.3%\" is not supported by Table 4. For ne, E∞ is 0.392%, and for Te, E2 is 1.18% and E∞ is 6.385%. The paper explains the Te peak as a division artifact near walls, and that explanation is reasonable, but the blanket claim should be qualified. Second, NH=8 is selected from an FFT of the same DTS reference used to measure the error. That makes the reported <0.3% partly in-sample. The NH=6,8,10,12 sweep is described qualitatively but no error-vs-NH table is given, so the reader cannot tell whether NH=8 is a true plateau or a tuned sweet spot. I would ask for a table of E2/E∞ as a function of NH for the conservative variables. Third, it is 1D only; the method is aimed at multidimensions, and no 2D evidence is provided. That does not invalidate the 1D claim, but it limits the strength of the \"paradigm\" language.\n\nThe periodicity assumption in Eq. (15) is standard for dissipative, periodically forced systems and I do not consider it a serious flaw here; the reader's concern about slow species is legitimate in general but not a reason to doubt this benchmark.\n\nBottom line: the math, the data, and the citation pattern look solid. The claims overreach modestly. With a qualified abstract and an NH-convergence table, this is a good contribution. I would bring it to the reading group and cite it as a useful benchmark for time-spectral plasma solvers.","headline":"Solid incremental work that is likely correct, but the headline <0.3% claim is in-sample and overstated; the paper deserves a careful revision, not a desk reject.","tokens_in":19372,"tokens_out":2660,"would_cite":true,"duration_ms":379628,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Harmonic balance solves RF plasma periodic state 10x faster than time-marching","keywords":["harmonic balance","capacitively coupled plasma","drift-diffusion-Poisson","electron-energy transport","time-spectral method","operator splitting","RF plasma simulation","local-mean-energy approximation"],"falsifier":"Run the same 1D argon benchmark at lower pressure (e.g., 0.1 Torr) or with a metastable population whose lifetime is many RF periods; if the HB solution with NH=8 differs from a fully converged DTS solution by more than 0.3% in the bulk plasma, the strict periodicity assumption is not generally valid.","tokens_in":18454,"feed_emoji":"⚡","tokens_out":1092,"duration_ms":16384,"temperature":0.7,"pith_summary":"This paper shows that a time-domain harmonic balance method can directly compute the periodic steady state of radio-frequency capacitively coupled plasmas, skipping the long transient phase that conventional time-marching solvers must integrate. The authors extend harmonic balance beyond simpler plasma models to include full electron-energy transport, using an operator-splitting scheme that keeps all implicit inversions cell-local. On a 1D argon benchmark, retaining eight harmonics resolves the quasi-steady bulk and the sharply nonlinear sheath dynamics, with macroscopic errors under 0.3% versus a highly converged dual-time-stepping reference. The same accuracy is reached in about a tenth of the wall-clock time of a fully converged baseline, even on a single CPU core.","feed_headline":"Plasma simulation skips transients with 10x speedup","feed_subtitle":"Time-domain harmonic balance nails the periodic state of RF CCPs with 0.3% error and a tenth of the wall-clock cost.","key_machinery":"The central mechanism is the Fourier-collocation time-spectral operator E, a real, dense, skew-symmetric matrix mapping time derivatives at the discrete phase points to a coupling across all collocation points (Eq. 23). Applying E to the drift-diffusion and electron-energy equations converts the unsteady system into a quasi-steady one, and the proposed implicit relaxation factorizes the local Jacobian (Eq. 39) into a spatial factor and a temporal factor, so the dense temporal inversion is only per-cell and per-variable.","core_discovery":"The paper establishes that the time-domain harmonic balance method can be extended to a fully coupled drift-diffusion-Poisson system with electron-energy transport (the local-mean-energy approximation) for RF capacitively coupled plasmas, and that with eight harmonics the periodic solution is strictly consistent with an established time-marching reference. The central technical contribution is converting the periodic-in-time plasma flow problem into a pseudo-steady system sampled at NT=2NH+1 collocation points per RF cycle, then solving it with a spatiotemporal operator splitting: a spatial implicit relaxation sweep followed by a cell-local dense temporal inversion that treats the time-spect","pith_inferences":["If the method holds in 2D/3D, where the same operator-splitting remains cell-local, the memory advantage over monolithic Jacobian approaches should become even more pronounced.","The 10x speedup for a single condition suggests that for optimization studies or design loops that require many successive RF conditions, the effective speedup could compound; a designed experiment varying pressure or voltage would test this.","Because the electron temperature is a derived ratio, the near-wall Te error (6.4%) hints that derived quantities in depleted sheaths will always be less accurate than conserved quantities; a post-processing smoothing or a different temperature definition might be worth testing.","A natural next test is a discharge with a slowly relaxing neutral metastable population whose period exceeds the RF period, since the method assumes an exactly T-periodic quasi-steady state and could fail or require treating slow manifolds separately."],"forward_implications":["For RF CCP simulations, the periodic state can be obtained without simulating hundreds of RF cycles, making routine parameter sweeps and reactor design iterations much cheaper.","Because the method evaluates nonlinear kinetics directly at each phase point, it extends naturally to more complex chemistries and multidimensional geometries without changing the core formulation.","The reported errors for ne, εe, and ϕ are all below 0.3% at NH=8, confirming that spectral truncation is not a bottleneck for this problem class.","The speedup is achieved on a single core, so the method can be combined with spatial parallelism for further gains.","The approach inherits harmonic balance’s clean error control: truncating at NH only changes the spectral error, not the physical time-step error, making it a suitable verification tool for time-marching solvers."],"fun_headline_variants":["RF plasma simulation achieves 10x speedup","Harmonic balance skips transients for 10x plasma sim speed","Plasma simulation: 10x faster by skipping transients","Time-domain harmonic balance: 10x faster RF plasma sims","Simulating RF plasmas 10x faster without physical transients"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The method assumes the plasma reaches a strictly periodic steady state with period T of the RF drive; if any species or coupled dynamics relax on a longer timescale, the computed periodic solution would not match the time-asymptotic state of a time-marching simulation.","fun_headline_variants_meta":{"raw":{"variants":["RF plasma simulation achieves 10x speedup","Harmonic balance skips transients for 10x plasma sim speed","Plasma simulation: 10x faster by skipping transients","Time-domain harmonic balance: 10x faster RF plasma sims","Simulating RF plasmas 10x faster without physical transients"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000237,"raw_usage":{"total_tokens":1377,"prompt_tokens":810,"completion_tokens":567,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":554,"completion_tokens_details":{"reasoning_tokens":480}},"tokens_in":554,"tokens_out":567,"duration_ms":5822,"temperature":1.0,"reasoning_tokens":480,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-01T15:58:42.538935+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the same 1D argon benchmark at lower pressure (e.g., 0.1 Torr) or with a metastable population whose lifetime is many RF periods; if the HB solution with NH=8 differs from a fully converged DTS solution by more than 0.3% in the bulk plasma, the strict periodicity assumption is not generally valid.","supporting_citations":[],"review_version":1}