Runge–Kutta methods at the complete-positivity boundary

In a closed classical calculation, a time-stepping method is usually judged by consistency, convergence order, stability, and computational cost. Those diagnostics remain necessary in quantum dynamics, but the object being propagated is not an arbitrary vector. A density matrix must remain Hermitian, have unit trace, and be positive semidefinite.

Positivity alone is not the final condition. A channel may act on one part of an entangled system, so its extension by the identity map must also preserve positivity. This stronger requirement is complete positivity. Together with trace preservation, it defines a CPTP map, the standard mathematical representation of a physically admissible quantum channel.

A direct Runge–Kutta update of a Lindblad equation preserves trace in the model considered here, and real coefficients preserve Hermiticity. Neither fact guarantees complete positivity. The remaining failure can be small enough to disappear on an ordinary error plot while still giving the Choi matrix a negative eigenvalue.

The Choi construction converts the tensor-level definition of complete positivity into a finite matrix test. Apply the channel to one half of a maximally entangled state. The channel is completely positive exactly when the resulting Choi matrix is positive semidefinite.

For zero-temperature qubit damping, the channel is controlled by two real multipliers. The excited-state population is multiplied by a, while the off-diagonal coherence is multiplied by c. The Choi condition reduces to

The first inequality prevents a population from becoming negative or exceeding its physical ceiling. The second couples coherence to population. A channel cannot retain more squared coherence than the surviving population can support.

Consider a qubit relaxing from its excited state at rate γ. For a step of durationh, define the dimensionless relaxation step x = γh. Optional pure dephasing contributes a second coordinate u = γ_ϕh. The exact population and coherence multipliers are

With pure amplitude damping, u = 0 and |c|² =a exactly. The physical evolution runs along the boundary of the channel set rather than through its interior. There is no structural buffer. A small numerical perturbation can land on either side.

Added dephasing changes that geometry. It suppresses coherence more rapidly while leaving the population multiplier unchanged, so |c|² < a. Near the origin, every fixed positive dephasing-to-relaxation ratio creates a nonzero local CPTP interval. Pure amplitude damping is therefore the critical case: it exposes the sign of the numerical defect without an interior margin masking it.

For the scalar test equation y′ = λy, a Runge–Kutta method replaces the exact exponential e^hλ with its stability function R(hλ). In this damping model, the Liouvillian separates into invariant population and coherence modes. The full numerical channel is therefore determined by two evaluations of the same function:

On the pure-amplitude-damping boundary, the complete-positivity test becomes

The scalar margin M_R is the determinant of the active Choi block. Its sign is the sign of the lower Choi eigenvalue. This makes the physicality test exact without repeatedly diagonalizing a matrix.

The population ceiling and the block margin are independent. Forward Euler atx = 1 keeps both scalar modes within the usual absolute-stability bound, yet its margin is −1/4. Heun’s second-order method has a positive block margin at x = 4, but its population multiplier is 5. Either condition can be the first route out of the channel set.

Suppose the first difference between the stability function and the exponential occurs at order m:

Substitute one full population step and two half coherence steps into the Choi margin. The exact exponential terms cancel because e^−x = (e^−x/2)². What remains is

The factor 1 − 2^1−m is familiar from full-step versus two-half-step error estimation. Here its signed coefficient has a geometric meaning: it tells us whether the numerical channel initially enters the CPTP region or exits it.

This is not an odd-versus-even order rule. The pattern happens to hold for minimal truncated exponential polynomials, but the Dormand–Prince principal fifth-order formula crosses inward. Backward Euler is first order and remains admissible for every nonnegative step in this model; implicit midpoint is second order and crosses outward immediately. The signed first defect, not formal order alone, controls the local direction.

The release evaluates eight representative explicit and implicit formulas using exact rational Butcher coefficients. Boundary polynomials are derived symbolically. Positive roots are isolated with Sturm counts, and exact rational samples certify the sign on every intervening interval.

The resulting admissible sets include four qualitatively different forms: the origin alone, a finite interval attached to the origin, an unbounded interval, and a disconnected set. The embedded fourth-order Dormand–Prince formula is the most instructive global example. It is non-CPTP for every sufficiently small positive step, then becomes CPTP on the isolated interval 3.0688565481… ≤ x ≤ 4.3849863208….

That remote component is a physicality statement, not an accuracy recommendation. Atx = 4, the step is CPTP but its population multiplier differs substantially from the exact exponential. Complete positivity constrains structure. It does not certify that a numerical answer is close.

Repeating a fixed phase-covariant step raises the population multiplier toa^N and the coherence multiplier to c^N. On the nonnegative-population branch, taking powers preserves the ordering betweena and c². A failed one-step map therefore remains failed after any finite number of identical compositions.

Fixed-horizon refinement is subtler. Hold a final dimensionless time X constant and take N steps of size x = X/N. For outward methods, the negative Choi eigenvalue can decay at the method’s ordinary global order. The violation becomes rapidly smaller, but it remains negative for every sufficiently fine finite refinement.

This is the central separation. Convergence describes distance to the exact channel. Complete positivity describes which side of the channel boundary the approximation occupies. A sequence can improve by every conventional error measure while approaching the boundary entirely from the nonphysical side.

The repository is not a collection of plotting scripts attached after the analysis. It is the computational record of the proof obligations. One command reconstructs each stability function from its exact tableau, identifies the first exponential defect, forms the boundary polynomials, verifies all listed positive-root brackets, rebuilds the interval sign charts, and regenerates the figures and fixed-horizon data.

A successful check certifies eight method audits, fourteen positive-root certificates, and three figures. The exact symbolic and Sturm calculations are recorded in machine-readable JSON; a compact CSV provides the method-level classification; a human-readable supplement presents the root isolation and sign charts. No external dataset, web access, or proprietary software is required.

This division matters. Decimal endpoints are useful for reading a figure, but they are not the proof. The proof resides in the polynomial identities, exact rational isolating intervals, Sturm root counts, and rational sign samples. The repository keeps those layers together.