{"id":"5ce80660-a55a-404c-a686-738862adc1e2","arxiv_id":"2507.13997","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Two methods, an asymptotic expansion and a predictor-corrector scheme, compute slow manifolds defined by zero fast isostable coordinates and support reduced-order models of nonlinear dynamics.","lead":"This paper presents two new numerical methods to compute slow manifolds of nonlinear dynamical systems using isostable coordinates, which track the slowest decay modes near a stable fixed point. The methods allow backward-time integration along the manifold to model nonlinear behaviors such as period-doubling bifurcations that linear approximations miss.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The predictor-corrector's key residual assumption in Eq. (42) already fails in the paper's planar example; the method's success there appears to rely on forward contraction rather than accurate backward propagation along W^s.","rationale":"The reader's weakest_assumption identifies the same mechanism I find most load-bearing: the predictor-corrector hinges on Eq. (42), whose residual bound r_j=O(ε) is not established and is acknowledged by the authors as difficult to check. My independent analysis sharpens this. In the derivation, an O(ε) error in g_j(t2) that has any component in the fast spectral subspace is amplified by the inverse of an O(ε) eigenvalue, so it produces O(1) residual at t1. Biorthogonality alone does not prevent this; an additional normality or invariant-subspace condition is needed. The paper's own planar example allows this to be checked in closed form, and the check shows the residual is O(1) at the far end of the claimed backward integration. The method still produces useful reduced-order models in the examples because the correction step and exponential attraction of W^s in forward time mask the tangent error. This does not invalidate the numerical results, but it means the examples do not verify the stated backward-propagation mechanism, and the method has no quantitative error control in systems with weaker contraction. This supports the reader's CONDITIONAL verdict; I would not change it, but the condition should include an explicit check of Eq. (42)'s residual along the computed trajectory.","tokens_in":19049,"tokens_out":15935,"duration_ms":195547,"concrete_test":"Re-run the planar example (47) with β=1 and compute the residual exactly. Use the exact slow manifold x2=(5/4)x1^4−(20/9)x1^2 and exact g1=(1,h'(x1)), h'=5x1^3−(40/9)x1, along the predictor-corrector backward trajectory; plot ||(Id−v1v1†)g1||_2 as a function of x1. If this quantity is O(1) near x1≈1.8 while the initial ψ1 is 0.001, Eq. (42)'s validity condition fails. Then check whether the forward-in-time orbit from the computed endpoint stays within a prescribed tolerance of the exact manifold; if it does, the observed accuracy is due to contraction, not to the assumed residual bound. The same residual check can be run for the Goodwin example using the eighth-order expansion of g1 as an independent proxy.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 3.7's approximation (42) is load-bearing. It replaces g_j(t1) by a vector in Vβ after asserting r_j(t1)=(Id−PβPβ†)g_j(t1) has O(ε) norm and σ_min of [v̂_{β+1}...v̂_N] is O(1). The derivation (39)–(41) does not justify that residual bound. Writing g_j(t2)=v_j+O(ε), any O(ε) component of this correction lying in the fast spectral subspace is multiplied by 1/(λ̂_k exp(−λ_j(t2−t1))) ∼ O(1/ε) during inversion of Φ_gj, and therefore contributes O(1) to r_j(t1). Biorthogonality annihilates v_j itself but not the O(ε) correction. So (42) tacitly assumes the O(ε) correction has no fast component—an additional spectral/nonnormality condition that is neither stated nor checked. Section 4.1 gives a direct counterexample. For (47) with β=1, ψ1=x1 and the exact slow manifold is x2=(5/4)x1^4−(20/9)x1^2, so g1=(1,h'(x1)) with h'=5x1^3−(40/9)x1. Starting near x1=0.001 and integrating backward to τ=150 reaches x1≈1.8, where ||r_1||_2=|h'(1.8)|≈21.2, not O(ε) for ε=0.001. Thus Eq. (42) is violated exactly in the nonlinear regime the method claims to enter. The favorable Figure 3 result appears to be produced by the correction step (46) plus exponential contraction in forward time, not by accurate backward propagation on W^s.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper defines the slow manifold of a stable fixed point as the set where the fast principal isostable coordinates vanish, W^s = {x | ψ_k(x)=0 for k>β}, and proposes two computational strategies for following this manifold backwards in time from the fixed point into the nonlinear regime. The first strategy (§3.6) uses an asymptotic expansion of the dual vectors g_j in powers of the isostable coordinates; the second (§3.7) is a predictor-corrector method that approximates g_j by its projection onto the slow eigenspace V_β (Eq. (42)). The backward-time evolution equation (19) is exact given the gradients I_k, and the paper derives evolution equations for I_k and g_k along trajectories. The methods are applied to a planar model, a damped pendulum, the Goodwin oscillator, and a population of coupled planar oscillators, and reduced-order models are constructed that reproduce forced responses and a period-doubling bifurcation missed by linearization.","tokens_in":19464,"tokens_out":13906,"duration_ms":153494,"significance":"The idea of defining slow manifolds through principal isostable coordinates is attractive, and the exact formulation in Eq. (19) is a useful contribution that connects Koopman/isostable theory with geometric slow-manifold reduction. The numerical examples span nontrivial systems, and the reduced-order models show behavior, such as a period-doubling bifurcation and frequency-dependent resonance, that a local linearization misses. However, the central approximation (42) of the predictor-corrector method is an ansatz whose validity the paper itself, in Section 5, says cannot be checked explicitly; moreover, Section 4.1 provides a direct instance where the residual assumption behind (42) fails badly. The advertised claim of accurate backward propagation along W^s is therefore not established, and the numerical successes appear to rely on the forward contraction property rather than on accurate computation of the slow manifold. The paper would be substantially strengthened by reframing what the predictor-corrector algorithm actually computes and by supplying quantitative error measures.","major_comments":[{"comment":"The approximation g_j(t1) ≈ ∑_{k=1}^β (1/(\\hat{λ}_k exp(-λ_j(t2-t1)))) \\hat{v}_k \\hat{w}_k^T v_j is load-bearing for the predictor-corrector strategy, but the residual condition ||r_j(t1)||=O(ε) is asserted rather than derived. From (39), r_j(t1) contains the terms (1/(\\hat{λ}_k exp(-λ_j(t2-t1)))) \\hat{v}_k \\hat{w}_k^T(v_j+O(ε)) for k=β+1,...,N. Even when \\hat{w}_k^T v_j=0 by biorthogonality, the O(ε) correction at t2 generally has a fast component, and that component is amplified by the reciprocal of \\hat{λ}_k exp(-λ_j(t2-t1)), which is exponentially large for k>β. Thus (42) tacitly imposes an additional spectral/nonnormality condition on the O(ε) correction. Section 5 concedes that there is no obvious metric to check this condition, so the central approximation is unverified.","section":"3.7.1, Eq. (42)"},{"comment":"For β=1, the planar system (47) has exact slow manifold x2=(5/4)x1^4-(20/9)x1^2 and hence g1=(1,h'(x1))^T with h'=5x1^3-(40/9)x1. The predictor-corrector integration starts near x1=0.001 and, after 150 time units of backward integration, reaches x1≈0.001 exp(0.05·150)=1.8. At x1=1.8, ||r_1||_2=|h'(1.8)|≈21.2, which is not O(ε) for ε=0.001. Equation (42) is therefore violated precisely in the nonlinear regime that the method claims to enter. The favorable result in Figure 3 appears to be produced by the correction step (46) followed by exponential contraction in forward time, not by accurate backward propagation on W^s. The manuscript should either demonstrate backward propagation with quantitative error measures or redefine what the algorithm actually computes.","section":"4.1, Eq. (47)"},{"comment":"The error analysis for the computation of I_1,...,I_β is heuristic. Equation (27) bounds the error under the assumptions that per-step solver errors are O(ε), that \\bar{s}_k||I_k||<O(1/ε), and that the relevant exponential factors remain O(1). The second assumption is state-dependent and is not verified along the computed trajectories, and the analysis does not account for the fact that the trajectory itself is approximate when (19) is integrated. Since (19) requires I_1,...,I_β, this gap is load-bearing. A numerical estimate of the error in I_j along the manifolds of Section 4 would help.","section":"3.4, Eq. (27)"},{"comment":"The numerical validation is qualitative throughout. For the planar example, the exact slow manifold is known and the errors of both methods can be tabulated as functions of ψ_1 and of the expansion order, but no error norms are reported. For the pendulum, Goodwin, and coupled-oscillator examples, agreement is assessed visually by forward integration, and the period-doubling bifurcation in Section 4.3 is located at a=0.020 versus the full model's a=0.023 without uncertainty quantification. Given that the checkable hypotheses behind (42) are admitted in Section 5 to be unverifiable, quantitative error measures are needed to support the paper's accuracy claims.","section":"4, Figures 3-9"}],"minor_comments":[{"comment":"The phrase 'in terms a linear' should read 'in terms of a linear'.","section":"Abstract"},{"comment":"In the linear solution formula, the sum is over j but the exponential uses λ_k; it should be exp(λ_j t).","section":"Eq. (14)"},{"comment":"The third-order term in the expansion should involve h_{ijk}, not h_{jk}.","section":"Eq. (34)"},{"comment":"The sentence referring to 'initial conditions evolving under the flow of (48)' should refer to the Goodwin model (17), not the pendulum (48).","section":"4.3, Figure 7 text"},{"comment":"The set {\\hat{v}_{β+1},...,\\hat{v}_{β+N}} should be {\\hat{v}_{β+1},...,\\hat{v}_N}.","section":"3.7.3, Eq. (46)"}],"recommendation":"major_revision","confidential_remarks":"The planar example in §4.1 is presented as a success, but the same computation demonstrates that the residual assumption behind Eq. (42) fails in the regime the method claims to handle. This should be addressed head-on in revision, either by proving additional conditions for the examples or by explicitly reframing the algorithm as a forward-contraction-based approximation rather than accurate backward propagation on W^s."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this paper deserves peer review, but the main new method has a load-bearing unsupported assumption, and the examples are not as clean as the text makes them sound.\n\nThe genuinely new thing is the combination of isostable coordinates with backward-time integration along the slow manifold W^s = {ψ_k=0, k>β}. Equation (19) is a clean and correct formulation. The examples are well chosen: the period-doubling capture in the Goodwin model and the resonance shift in the coupled oscillator population are exactly the nonlinear effects that linearization misses. The author is admirably explicit in Section 5 that the validity of the approximation cannot be checked.\n\nThe problem is in Section 3.7. Equation (42) replaces g_j(t1) by its projection onto the slow eigenspace after asserting that the residual r_j(t1) is O(ε). That assertion is not derived, and it is actually false in the paper's own planar example. For (47), the exact slow manifold is x2 = (5/4)x1^4 - (20/9)x1^2, so g1=(1,h'(x1)) with h'=5x1^3-(40/9)x1. Starting at x1=0.001 and integrating backward to t1=150 gives x1≈1.8, where |h'|≈21. So the residual is O(1), not O(ε). The favorable Figure 3 result is doing what the author admits in Section 5: forward-time contraction masks the error. The backward trajectory is not on the slow manifold; it is a nearby curve that converges to it forward in time. For the stated goal of 'approximate backward-time solutions on the slow manifold,' that is a mismatch.\n\nThis is not fatal to the model-reduction use case. The reduced-order models are validated against full simulation, and the examples show real improvement over linearization. But the paper should either justify (42) with a spectral/nonnormality condition, or reframe what the predictor-corrector actually computes, and add quantitative error measures against known slow manifolds or high-order expansions. There is also no code or convergence study, which would help a lot.\n\nThe asymptotic expansion method is less problematic but is a direct application of prior Taylor-expansion machinery. The novelty is mainly in the predictor-corrector.\n\nVerdict: send to serious review, but expect major revisions.","headline":"Genuinely useful isostable-based computation of slow manifolds, but the predictor-corrector's key approximation is not valid in the nonlinear regime and needs serious revision.","tokens_in":19923,"tokens_out":4503,"would_cite":true,"duration_ms":49505,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["34C45","37M99","37C10"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that slow manifolds of stable fixed-point systems can be computed backward in time with isostable coordinates, well beyond the linear regime.","keywords":["slow manifolds","isostable coordinates","Koopman operator","model order reduction","backward-time integration","predictor-corrector methods","asymptotic expansion","stable fixed points"],"falsifier":"Take a system whose slow manifold is known exactly, such as the planar model $\\dot{x}_1=-0.05x_1$, $\\dot{x}_2=-(x_2-x_1^4+2x_1^2)$, whose slow manifold is the nullcline $x_2=x_1^4-2x_1^2$. Integrate the predictor-corrector scheme (19), (20), (42), (46) backwards far enough into the nonlinear regime and compare the resulting trajectory to the exact curve: if the computed trajectory departs from the known curve before the stated integration range is reached, the projection assumption (42) has failed, and the central claim is refuted for that regime.","tokens_in":18814,"feed_emoji":"⏪","tokens_out":6771,"duration_ms":76671,"temperature":0.7,"pith_summary":"This paper tries to establish that the slow manifold of a dynamical system with a stable fixed point\\u2014the invariant surface on which the fast-decaying isostable coordinates vanish\\u2014can be computed explicitly, not just near the fixed point but far into the nonlinear regime. The obstacle is that backward-time integration amplifies any error at the fast decay rates, so naive integration leaves the manifold instantly. The paper develops two ways around this: an asymptotic expansion of the manifold in a basis of isostable coordinates, and a predictor-corrector scheme that repeatedly projects the solution back onto the slow eigenspace. If the methods work, they give reduced-order models that reproduce nonlinear behaviors such as period-doubling bifurcations in forced oscillations.","feed_headline":"Two schemes trace slow manifolds past the linear regime","feed_subtitle":"Isostable-coordinate reduced models reproduce nonlinear forced responses, including period-doubling, that linearization misses.","key_machinery":"The central objects are the principal isostable coordinates $\\psi_1,\\ldots,\\psi_N$\\u2014level-set coordinates of the slowest Koopman eigenfunctions, ordered so that each obeys $\\dot{\\psi}_k=\\lambda_k\\psi_k$\\u2014and their spatial gradients $I_k=\\partial\\psi_k/\\partial x$. The slow manifold is defined as the zero level set of the fast coordinates, $W^s=\\{x:\\psi_k(x)=0,\\ k>\\beta\\}$. The argument runs through equation (19), an ODE whose right-hand side is the inverse of the matrix of gradients $I_k$ applied to the vector $(-\\lambda_1\\psi_1,\\ldots,-\\lambda_\\beta\\psi_\\beta,0,\\ldots,0)$; this integrates trajectories exactly along $W^s$ in backward time. Because the slow gradients can be integrated reliably through (20), the only missing ingredient is the span of the fast gradients, supplied either by the asymptotic expansion (35) for $g_1,\\ldots,g_\\beta$ or by the predictor-corrector projection (42) combined with correction step (46).","core_discovery":"The central claim is that trajectories lying on the slow manifold $W^s = \\{x \\mid \\psi_k(x)=0 \\text{ for } k>\\beta\\}$ can be propagated backward in time from a small neighborhood of the fixed point into the strongly nonlinear region by integrating equation (19), which expresses $dx/d\\tilde{t}$ using the gradients $I_k$ of the slow isostable coordinates and zeros out the fast coordinates. The slow gradients $I_1,\\ldots,I_\\beta$ are computable accurately because adjoint equation (20) suppresses numerical errors over the integration window, while the fast gradient directions are never integrated directly: they are replaced by the span of $g_1,\\ldots,g_\\beta$, obtained either from a Taylor expansion of the state in isostable coordinates or from the predictor-corrector approximation $g_j(t_1) \\approx \\sum_{k=1}^\\beta (1/(\\hat{\\lambda}_k \\exp(-\\lambda_j(t_2-t_1)))) \\hat{v}_k \\hat{w}_k^T v_j$, which projects onto the slow linear eigenspace. The paper demonstrates on a planar example, a damped pendulum, the Goodwin circadian model, and a ten-oscillator network that the computed surface is invariant under forward flow and yields accurate reduced-order forced models, including reproducing a period-doubling bifurcation that a linearized model misses.","pith_inferences":["Beyond the paper itself, one could map where the projection assumption (42) breaks down by computing the residual $r_j(t_1)$ along the same trajectory, giving a practical validity horizon for the manifold that could be compared across systems with different spectral gaps.","Beyond the paper itself, a data-driven variant could estimate $I_1,\\ldots,I_\\beta$ and the span of $g_1,\\ldots,g_\\beta$ from short forward simulations, extending the approach to systems whose governing equations are unknown.","Beyond the paper itself, the correction-step subspace may select one branch of the slow manifold when resonance or near-resonance among eigenvalues makes the expansion (33) non-unique, so uniqueness questions could be studied by varying that projection."],"forward_implications":["If the methods hold, a slow manifold defined by $\\psi_k=0$ for $k>\\beta$ can be traced from the fixed point into the nonlinear regime, giving a concrete geometric object for model reduction rather than a local linear approximation.","The one- or two-dimensional isostable reduced models built from the computed manifold reproduce nonlinear forced responses, such as the period-doubling bifurcation of the Goodwin oscillator, that local linearization misses.","The predictor-corrector strategy only needs local Jacobian evaluations along a trajectory, so it scales to higher-dimensional systems where the asymptotic expansion becomes computationally prohibitive, as demonstrated on a 20-dimensional oscillator network.","Because forward trajectories converge exponentially fast to the slow manifold, moderate errors in the computed manifold still give useful reduced-order models."],"supporting_citations":[{"why":"Defines the principal isostable coordinates and their decay dynamics, which ground the slow-manifold definition.","marker":"[21]"},{"why":"Supplies the asymptotic expansion of the state in a basis of isostable coordinates and the strategy for computing its coefficients.","marker":"[35]"},{"why":"Provides the gradient dynamics $\\dot{I}_k=-(J^T-\\lambda_k I)I_k$ used to compute slow isostable gradients.","marker":"[39]"},{"why":"Establishes existence and uniqueness of global Koopman eigenfunctions for stable fixed points, justifying the principal isostable coordinates.","marker":"[18]"},{"why":"Gives the numerical solver error form used in the error analysis of Section 3.4.","marker":"[2]"},{"why":"Supplies the Goodwin oscillator model used as the main circadian testbed for slow-manifold computation and reduced-order modeling.","marker":"[10]"}],"fun_headline_variants":["Isostable coordinates expose slow manifolds past linear regime","Backward-time integration finds slow manifolds beyond linearity","Reduced models via isostable coordinates handle nonlinear dynamics","New computational routes to slow manifolds with isostable coordinates","Slow manifold identification using isostable coordinate system"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing assumption is that, well away from the fixed point, each fast-direction vector stays close enough to the linear slow eigenspace that it can be replaced by its projection; the paper itself notes there is no obvious metric for how far into the nonlinear regime this remains true.","fun_headline_variants_meta":{"raw":{"variants":["Isostable coordinates expose slow manifolds past linear regime","Backward-time integration finds slow manifolds beyond linearity","Reduced models via isostable coordinates handle nonlinear dynamics","New computational routes to slow manifolds with isostable coordinates","Slow manifold identification using isostable coordinate system"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000158,"raw_usage":{"total_tokens":1238,"prompt_tokens":968,"completion_tokens":270,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":584,"completion_tokens_details":{"reasoning_tokens":190}},"tokens_in":584,"tokens_out":270,"duration_ms":3743,"temperature":1.0,"reasoning_tokens":190,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T16:11:22.444623+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a system whose slow manifold is known exactly, such as the planar model $\\dot{x}_1=-0.05x_1$, $\\dot{x}_2=-(x_2-x_1^4+2x_1^2)$, whose slow manifold is the nullcline $x_2=x_1^4-2x_1^2$. Integrate the predictor-corrector scheme (19), (20), (42), (46) backwards far enough into the nonlinear regime and compare the resulting trajectory to the exact curve: if the computed trajectory departs from the known curve before the stated integration range is reached, the projection assumption (42) has failed, and the central claim is refuted for that regime.","supporting_citations":[{"cited_title":"Mauroy, I","cited_arxiv_id":null,"evidence_quote":"Defines the principal isostable coordinates and their decay dynamics, which ground the slow-manifold definition."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the numerical solver error form used in the error analysis of Section 3.4."},{"cited_title":"Gonze, S","cited_arxiv_id":null,"evidence_quote":"Supplies the Goodwin oscillator model used as the main circadian testbed for slow-manifold computation and reduced-order modeling."}],"review_version":1}