Padé(0,2) cuts each directional step from 189.37 ms to 5.92 ms at 240 nodes per direction, a factor of 32. The run uses the same step size, Δt=0.05, with moderate correlation ρ(t)=0.8+0.1cos(0.7t). At every refinement level, the exponential and Padé errors match to three digits. The speed ratio therefore survives at fixed error for that accuracy.
Anyone evolving densities inside a calibration loop should take the speed and keep the diagnostic close. At ρ=0.95, the construction's promised nonnegativity is unavailable on every mesh Itkin and Kazbek test. Their abstract says as much.
One replacement, same framework
The solver advances a density under anisotropic diffusion with a cross term and jumps. It is a forward Fokker-Planck method on a tensor grid. The companion Diagonal Frog paper formed directional operators in divergence form, producing matrices that are either M-matrices or eventually exponentially positive. In the M-matrix case, the matrix exponential is entrywise nonnegative at every step. The EM case gains the same property above a computable threshold τ0.
Each directional substep previously used an Arnoldi approximation to the exponential. Krylov work dominated the bill at O(m²N+m³) per sweep. With tolerance 1e-8, the Krylov dimension m climbs from 7 at N=48 to 12 at N=240.
This paper replaces a single factor. The exponential e^{sA} becomes the rational map r(sA), evaluated through a banded solve at O(N). The spatial discretization remains, along with the zero-flux conservative closure, the implicit mixed-derivative factor handled by a factorized Picard iteration, and the symmetric Strang ordering of substeps. Mass conservation follows immediately. Any rational r with r(0)=1 preserves unit column sums when applied to a zero-column-sum generator. On the one-dimensional benchmark, measured mass defects are 1e-12 or better.
A fresh positivity argument supplies the paper's theoretical contribution. The paper uses no market data; its evidence comes from PDE benchmarks. These comprise Ornstein-Uhlenbeck in one dimension, a two-dimensional anisotropic Gaussian at ρ=0.95, and Kramers escape from a double well over 600 steps to T=30. It then runs a Péclet sweep, a stiff pure-diffusion spike, and one-dimensional jump-diffusion with Merton and Kou kernels. The code is posted, using Matlab R2025b on an Apple silicon machine.
Why does infinity control the sign?
For these generators, r(γL) approaches Π + r(∞)(I−Π) as the step grows. Here Π is the rank-one projector onto the stationary density. Large-step entrywise signs are therefore governed solely by r(∞), the damping applied to infinitely stiff modes. When 0 ≤ r(∞) < 1, positivity occupies a genuine one-sided window γ ≥ γr. When r(∞) < 0, every node with stationary mass below |r(∞)|/(1+|r(∞)|) acquires a negative diagonal entry.
Crank-Nicolson has r(∞) = −1, demanding stationary mass above one half at every node. Padé(0,2) is also second order and L-stable, while its r(∞) = 0. That single value separates the methods. Across all four advection-dominated meshes tested, from n=51 to 161, the Cayley window stays empty. Even after scanning every γ on the finest mesh, the factor's least negative entry is −7.7e-3. At γ = τ0, the exponential has already become nonnegative, yet the factor still contains an entry of −0.98. The window opens in one weakly eventually positive case, spanning [5.24e-3, 5.93e-3], only a factor of 1.13 in the substep.
Two features favor Padé(0,2) for the directional factor. Under refinement from n=51 to 161, its threshold shifts only from 3.115e-1 to 3.175e-1. The backward Euler threshold, by comparison, drops from 2.906e-2 to 5.237e-3. With sinh grading extended to hmax/hmin = 8.5, τ0 increases by 2.4x and γ0 by 4.2x, while γr moves by 0.6%. The threshold can be certified once per problem without becoming a mesh-by-mesh concern. Its direction remains unusual: positivity arrives for large steps and disappears as Δt approaches zero.
High correlation
At ρ=0.95, the composite admissible region must satisfy both the directional lower bound 2γr and the mixed factor's upper bound Θ. Across six meshes, 2γr lies between 2.067 and 2.169. Meanwhile, Θ increases from 2.121e-3 to 6.593e-2 and levels off near 5e-2. The endpoints remain about forty apart, with refinement leaving the gap open.
Only the first-order backward Euler variant creates a nonempty region. Its 2γ0 falls from 2.477e-1 to 1.643e-2, and the interval appears only from N=101 onward. The authors accurately call this the Bolley-Crouzeix trade: provably nonnegative or second order, not both.
The abstract acknowledges the empty interval while defending the criterion. It "serves to discriminate the well-behaved multiplicative factors from the stabilizing-correction schemes rather than to guarantee positivity." The evidence is the positivity table at N=120. Padé records min p of −5.91e-8, against −4.12e-8 for Crank-Nicolson. Hundsdorfer-Verwer reaches −4.83e-2. MCS reaches −2.32e+12, with an L2 error of 2.368e+12. Padé stands six orders of magnitude from HV and twenty from MCS.
The comparison needs its stated qualification. The explicit one-sided cross stencil drives the comparators along a joint-refinement path that the paper plainly says leaves their stability region. At the same mesh, halving the step brings HV back to 1.46e-2. The discrimination claim is consequently narrower than the table first suggests: multiplicative factors endure a step path outside the affine methods' tolerance. They also charge more per step. In the timing run with Δt=0.05 and ρ(t)=0.8+0.1cos(0.7t), Padé takes 5.92 ms at N=240, versus 4.04 for HV and 4.35 for MCS.
A stabilizing-correction stage takes the form (I−θΔtF)^{-1}[Y−θΔtFp]. Under the conservative closure, Fp sums to zero. Unless the datum is stationary for the block, the bracket is signed for every relevant nonnegative datum. An entrywise nonnegative factor cannot repair a signed operand. Products of nonnegative rational maps permit factor-level sign control; affine stage combinations lack it. The case for this scheme shape survives the absent composite guarantee.
Refine only time and error rises
One experiment should be read before implementation. Hold the mesh at N=64 in each direction, set ρ=0.8, integrate to T=0.1, and refine the step while advancing the mixed factor through the factorized Picard iteration. Observed temporal order is negative, ranging from −0.42 to −0.88. Error rises from 1.343e-3 to 1.525e-2. The iteration converges to a nearby system with an O(h³) defect per step and no Δt factor. Accumulation then scales as O(h³/Δt).
A direct mixed solve yields 2.00.
Symmetrizing the one-sided differences recovers second order while worsening positivity by two to four orders. On the finest mesh at ρ=0.95, min p moves from −5.9e-8 to −2.0e-5. Here, second order refers to joint refinement Δt ~ h, the path on which the two-dimensional benchmark measures a 2.11 rate. Refining both arguments together is essential because step refinement alone makes the scheme worse.
Two implementation details matter. Threshold certification required 40-digit arithmetic. The smallest projector entries have size e^{-36}, while the Padé denominator's condition number grows like γ²||L||². On fine meshes, a double-precision scan would incorrectly declare γr unbounded.
Jump blocks favor the exponential. For Merton, γr is 11.4; for Kou, it is 48.5. Those values sit one to two orders above the directional thresholds, whereas τ0 = γ0 = 0 exactly for these Metzler generators. Itkin and Kazbek call the exponential the natural choice there, "the opposite of the preference for the directional blocks that motivated this paper." A truncated Poisson series keeps that choice cheap.
The practical construction is mixed: Padé(0,2) for directional sweeps, an exponential for the jump block, joint refinement, and one high-precision threshold certification. The third paper of the series marks the empty high-correlation region for repair. Its flux-limited variant carries no lower step bound, and this paper expects that change to close the two-dimensional gap. A two-dimensional demonstration is still pending. At ρ=0.95, the present result delivers exact mass, negativity around the 1e-7 level, and the per-step cost reported above. For most calibration work, that remains a good trade.