Fixing Divergence in Carleman Linearization via Analytical Continuation

Abstract: Nonlinear differential equations play a crucial role in modeling a wide range of phenomena, yet their solutions remain notoriously difficult to obtain. With the rapid development of quantum computing, quantum algorithms for efficiently solving such equations are actively being explored. One promising approach is based on Carleman linearization, which transforms nonlinear differential equations into linear systems. However, this method suffers from exponential divergence beyond a certain time scale. By reformulating the solutions in terms of eigenvalues and eigenvectors, we identify that this divergence originates from the Laurent expansion outside its neighborhood of convergence. To address this issue, we insert a regularized function to the divergent solution hinted by analytical continuation. We validate this divergence-correction method on both the logistic equation and some other partial differential equations like KPP-Fisher equations and Phase-Field models under periodic conditions. We implement our method for the logistic equation using the Linear Combination of Unitaries (LCU) quantum algorithm, providing a detailed complexity and error analysis.

• The paper identifies Carleman divergence as evaluating a Laurent-type spectral expansion outside its convergence radius and replaces it with Möbius-map analytic continuation implemented through regularized eigenvalue filters.
• The regularized method stabilizes logistic, KPP–Fisher, and phase-field simulations for long times, with explicit map constraints such as c < (1−x₀)/x₀ and exponentially decreasing truncation error in the logistic case.
• The authors show that quantum LCU implementations can realize the eigenvalue transformation with total Toffoli cost O(ΛK² log K/εₘₑₐₛ), while highlighting unresolved issues in eigenbasis conditioning, conformal-map selection, and general block encodings.

Carleman linearization embeds a polynomial nonlinear ODE into an infinite-dimensional linear system, and has become the standard route toward quantum algorithms for nonlinear dynamics. Its well-known defect is that truncated Carleman solutions diverge exponentially beyond a finite time horizon, regardless of truncation order. This paper by Zhu et al. (2607.05873) traces this divergence to a concrete analytic cause—the evaluation of a Laurent-type spectral expansion outside its radius of convergence—and repairs it via analytic continuation implemented as an eigenvalue regularization. The result is a "regularized Carleman method" validated on the logistic equation, KPP–Fisher equations, and phase-field models, together with an LCU-based quantum implementation and resource analysis.

Origin of the divergence

For the logistic equation $\dot{x} = x(1-x)$, the Carleman lift uses variables $y_k = x^k$ satisfying $\dot{y}_k = k(y_k - y_{k+1})$. The infinite Carleman matrix $\mathbf{A}$ has point spectrum $\{k\}$ with binomial eigenvectors. Expanding the initial condition in this eigenbasis yields coefficients $a_k = (x_0/(1-x_0))^k$ (for $|x_0|<1$), so the first lifted component reproduces the exact solution

$$x(t) = -\sum_{k=1}^{\infty}\left(\frac{x_0 e^t}{x_0-1}\right)^k,$$

a geometric series in $\zeta = e^t$ converging only for $|\zeta| < (1-x_0)/x_0$. Since $y_k = x^k$0 grows without bound along the positive real axis, divergence at late times is unavoidable even though the underlying solution is globally analytic. This identifies the failure mode precisely: not an artifact of finite truncation, but a convergence-domain problem of the spectral representation itself—a diagnosis consistent with earlier empirical observations of divergence in lattice Boltzmann and related settings.

Analytic continuation via a Möbius map

The remedy is a conformal map $y_k = x^k$1 mapping the positive real ray to $y_k = x^k$2. Substituting into the divergent series and using the binomial expansion produces a series convergent on the entire unit disk—hence for all $y_k = x^k$3—provided

$y_k = x^k$4

an explicit constraint coupling the map scale parameter to the initial condition. Numerically, the correction takes a particularly simple form: inserting a regularized function $y_k = x^k$5 into each eigenmode of the divergent solution,

$y_k = x^k$6

where $y_k = x^k$7 reduces to the regularized incomplete beta function $y_k = x^k$8 for reaction–diffusion systems, equals 1 when the original series converges, and decays to 0 as $y_k = x^k$9, suppressing the divergent modes. Nonpositive eigenvalues (from diffusion terms) are left untouched since they are already convergent.

Numerical experiments confirm both the mechanism and the bound: for $\dot{y}_k = k(y_k - y_{k+1})$0 (critical $\dot{y}_k = k(y_k - y_{k+1})$1), $\dot{y}_k = k(y_k - y_{k+1})$2 requires large $\dot{y}_k = k(y_k - y_{k+1})$3; $\dot{y}_k = k(y_k - y_{k+1})$4 matches the exact solution already at $\dot{y}_k = k(y_k - y_{k+1})$5; $\dot{y}_k = k(y_k - y_{k+1})$6 shows oscillatory convergence; and $\dot{y}_k = k(y_k - y_{k+1})$7 diverges, exactly as predicted. Stability diagrams in the $\dot{y}_k = k(y_k - y_{k+1})$8 plane delineate monotonic, oscillatory, and divergent regimes, matching the analytical boundary. A practical caveat is that the mapped-series cutoff $\dot{y}_k = k(y_k - y_{k+1})$9 must be coupled to the Carleman order $\mathbf{A}$0: taking $\mathbf{A}$1 introduces additional error or instability, so $\mathbf{A}$2 is recommended rather than treating $\mathbf{A}$3 as freely optimizable.

Extension to PDEs

KPP–Fisher equations: after spatial discretization with periodic boundaries ($\mathbf{A}$4), the truncated spectrum consists of integers from $\mathbf{A}$5 to $\mathbf{A}$6 with high degeneracy, giving a Laurent-series structure; degeneracy is broken by a small perturbation $\mathbf{A}$7 to stabilize eigenvector computation. The regularized solution agrees with Euler-method reference integration up to $\mathbf{A}$8. The $\mathbf{A}$9 case has non-integer eigenvalues, but the gamma-function definition of $\{k\}$0 generalizes without modification—an important robustness check.

Cubic nonlinearities: for the cubic logistic equation $\{k\}$1, the standard map fails because the branch points $\{k\}$2 are mapped inside the unit disk for all $\{k\}$3. Introducing $\{k\}$4 moves the branch point onto the negative real axis, where the modified map keeps it outside the disk whenever $\{k\}$5. The regularized function becomes $\{k\}$6. This distinction is essential and carries over to the discrete phase-field model, where stability diagrams separate regimes according to the sign of the mean initial value.

Quantum implementation

Because the regularized solution is an eigenvalue transformation $\{k\}$7 of the Carleman matrix, it can be realized quantumly: block-encode the sparse Carleman matrix (sparsity $\{k\}$8 for the logistic case, built from diffusion, value, and position oracles following Camps et al.), approximate the transformed function by a degree-$\{k\}$9 polynomial, and implement it via LCU. For the logistic equation, interpolation on the integer spectrum $a_k = (x_0/(1-x_0))^k$0 gives $a_k = (x_0/(1-x_0))^k$1 with zero polynomial approximation error in exact arithmetic.

The error budget decomposes as mapped-series truncation error (bounded uniformly by $a_k = (x_0/(1-x_0))^k$2 over all time, hence requiring $a_k = (x_0/(1-x_0))^k$3), eigenvector-conditioning-amplified polynomial error $a_k = (x_0/(1-x_0))^k$4, block-encoding error $a_k = (x_0/(1-x_0))^k$5, plus measurement and state-preparation errors. With $a_k = (x_0/(1-x_0))^k$6, $a_k = (x_0/(1-x_0))^k$7, the total Toffoli cost scales as

$a_k = (x_0/(1-x_0))^k$8

A notable tension is flagged explicitly: while increasing $a_k = (x_0/(1-x_0))^k$9 suppresses the regularization error exponentially as $|x_0|<1$0, the LCU normalization factor $|x_0|<1$1 may grow with $|x_0|<1$2 and the interpolation coefficients, so parameter choices must balance these competing effects. A classical simulation of the full algorithm (Quri SDK, $|x_0|<1$3, $|x_0|<1$4, $|x_0|<1$5) reproduces the exact logistic solution without divergence.

Limitations and open questions

The authors are candid about scope. The explicit exponential error bound $|x_0|<1$6 holds only for the logistic equation; for reaction–diffusion systems, a complete analysis must additionally account for spatial discretization error, degeneracy-breaking perturbation error, and conditioning of the Carleman eigenbasis, none of which receive sharp a priori bounds here. Systematic selection of the optimal conformal map for a given singularity structure remains unsolved—the cubic case shows that map choice is nontrivial and solution-dependent. On the quantum side, efficient block encodings of Carleman matrices for general polynomial systems, phase-field models, and fluid equations are not yet available; the logistic construction exploits an unusually simple two-nonzero-per-row sparsity pattern. Application to computational fluid dynamics, where Carleman-based lattice Boltzmann simulations exhibit the same qualitative divergence, is proposed but untested.

Conclusion

This work converts a long-standing empirical pathology of Carleman linearization into a solved analytic problem: divergence is the Laurent expansion evaluated outside its convergence radius, and a Möbius-map analytic continuation—realizable as insertion of a regularized incomplete-beta eigenvalue filter—restores convergence for all time under explicit constraints relating the map scale to initial data. The method covers quadratic and cubic nonlinearities, extends to discretized reaction–diffusion and phase-field systems, and admits a concrete LCU-based quantum implementation with explicit complexity scaling $|x_0|<1$7. It thereby establishes a foundation for long-time stable Carleman-based simulation of nonlinear differential equations on both classical and quantum hardware, contingent on future progress in general block-encoding construction and systematic conformal-map selection.