18Simulating time evolution

Turning energy into a phase

We still want the electronic physics of neutral singlet N₂ at 1.85 Å, in the chemical model selected earlier. Its ordinary qubit Hamiltonian describes the active operator, with the nuclear/frozen-inactive core scalar retained separately. A trial-state preparation circuit supplies an initial state. What should happen next depends on the question: estimate an eigenenergy, or follow the evolution of an initial state.

For a time-independent Hermitian Hamiltonian, the evolution operator is

The reduced Planck constant $\hbar$ converts energy times time into phase. The repository’s atomic-unit convention takes $\hbar=1$, so its formulas use $e^{-iHt}$. The parameter time is simulated physical time in the model’s units, not seconds of hardware wall-clock runtime. A longer simulated time can require more gates and tighter synthesis, but those are derived costs, not the value of time itself.

If $H|E_k\rangle=E_k|E_k\rangle$, then $U(t)|E_k\rangle=e^{-iE_kt}|E_k\rangle$. This is how energy becomes a phase. In a superposition of eigenstates, the relative phases change with their energy differences; that produces real dynamics. A controlled evolution makes an eigenphase observable relative to a control-zero branch, supplying the central ingredient for later phase estimation.

Adding $E_{\mathrm{core}}I$ multiplies an uncontrolled evolution by $e^{-iE_{\mathrm{core}}t}$. That global phase is unobservable in isolated state dynamics but becomes relative under control. The ordinary mapping route deliberately leaves the core scalar outside the Pauli operator, so its measured active energies require the corresponding classical addition once. Do not remove other identity terms without tracking their phases and then assume a controlled circuit is unchanged.

This lesson covers the actual fixed-Hamiltonian builders. Explicitly time-dependent Hamiltonians and propagators belong to the later dynamics lesson. No evolution, scientific calculation, circuit compilation, or test is being executed here.

Why one exponential becomes many

Write the represented qubit Hamiltonian as

where $P_j$ is a Pauli word, $h_j$ a real coefficient, and $L$ the number of listed terms after the relevant screening. Since $P_j^2=I$, one factor has the exact form

Lesson 17 showed how basis changes, a parity network, a rotation by $2\theta$, and uncomputation implement such a factor. The remaining problem is that the exponential of a sum is generally not the product of exponentials.

The commutator $[A,B]=AB-BA$ measures whether the order of two operations matters. Commuting operators have a common compatible action and can be exponentiated separately without this splitting error. Anticommuting Pauli words produce $[h_jP_j,h_kP_k]=2h_jh_kP_jP_k$, with norm $2|h_jh_k|$. The complex phase in $P_jP_k$ matters. The Pauli commutator utility uses the native accumulator for products and reverses its canonical labels into the QubitOperator string convention.

As a toy model, take $H=X+Z$ on one qubit. Because $[X,Z]=-2iY$, X evolution followed by Z evolution is not exact combined evolution. The exact expression is

A first-order step instead executes X then Z, whose matrix is $e^{-i\delta Z}e^{-i\delta X}$. Its short-time expansion contains a leading commutator contribution absent from the target. This is an analytic toy illustration, not a molecular or simulator result.

The distinction between chronology and written multiplication remains important: a list of gates is executed first-to-last, while matrices acting on a column state multiply rightmost-first. The builders return chronological ExponentiatedPauliTerm lists. Their angle already means the coefficient in $e^{-i\theta P}$; the circuit mapper’s Q# call uses Exp(paulis,-angle,systems).

Trotter and Suzuki product formulas

Split total time into $N$ slices of size $\delta=t/N$. The first-order Lie–Trotter formula executes each term for $\delta$ per slice. For two terms $A$ and $B$, chronological A then B gives

The local error begins at order $\delta^2$; repeated unitary slices give a leading total scale proportional to $t^2/N$, with commutator-dependent constants. If all terms commute, the product is exact apart from numerical/synthesis errors and coefficient screening.

Second-order Strang splitting uses a palindrome: A for half a slice, B for a whole slice, A for half a slice. Its symmetry cancels the leading asymmetric error:

For many groups it sweeps forward with half times, applies the last group at full time, and sweeps back. The leading local error is cubic, and the leading total scaling is $|t|^3/N^2$. This does not mean “twice as accurate” for a fixed gate budget: one second-order slice contains more factors.

The native trotter builder supports order 1, order 2, and any positive even integer above 2, despite an older settings comment saying only first order. Higher orders use Suzuki recursion:

Here $u_k$ sets the fractional time of each subformula. Some fractions are negative, meaning inverse-time substeps, not negative physical runtime. Each recursion multiplies the unmerged schedule length by roughly five. The implementation merges consecutive entries for the same group before expanding them into Pauli rotations. Higher order can save slices while increasing per-slice length and classical schedule construction.

The algorithm family is hamiltonian_unitary_builder, whose default is trotter. Relevant defaults are time=0.0, order=1, num_divisions=0, target_accuracy=0.0, error_bound="commutator", and weight_threshold=1e-12. Zero time means identity evolution; it is not a sensible phase-to-energy decoding interval. A positive explicit division count is a floor; nonpositive values fall back to one slice when automatic selection is disabled. A positive target accuracy activates the selected step estimator, and the larger of manual and automatic counts wins. Nonpositive target accuracy disables automatic selection in the builder, although calling the error helper directly with a nonpositive accuracy raises.

All these builders require a Pauli-decomposition QubitOperator, not an SOS container. Trotter checks coefficient Hermiticity to its threshold, then uses real coefficients above that cutoff. The returned PauliProductFormulaContainer contains a single step list, step_reps, num_qubits, and scale. Constructing it neither simulates the state nor synthesizes the final gates.

Grouping is a mathematical assumption, not just a label

The shared time-evolution base consumes an existing FlatPartition or LayeredPartition; without one, each listed Pauli term is its own group. It materializes subgroup operators, removes empty groups, and sorts outer groups by their number of subgroups. A flat partition therefore follows its stable group ordering, while a layered partition can be reordered by this rule.

Within each subgroup, the code emits separate Pauli exponentials through a helper named _exponentiate_commuting. That helper does not itself verify commutation. Exact subgroup exponentiation and the nominal Suzuki order require the intended compatible grouping assumptions. A manually attached partition containing noncommuting terms is not repaired or promoted into an exact group exponential. Layers are a scheduling structure, not proof that their combined group commutes.

Error helpers generally inspect the screened original term list rather than reconstructing and certifying every grouped schedule the builder may produce. Especially for order-dependent second-order constants or unusual layered groups, preserve and check the actual schedule instead of treating a partition strategy string as a proof of the requested order.

This also constrains the CSWAP controlled mapper from Lesson 17. Its computational vacuum must be preserved by the ordered finite product, not only by the unsplit chemical Hamiltonian. Randomizing or regrouping terms can destroy cancellation partners that made the shortcut valid.

What the Trotter accuracy helpers actually compute

An operator norm measures the largest vector amplification of an operator; for a unitary difference it bounds worst-case state-vector error. A dimensionless target $\epsilon$ for $\|U-\widetilde U\|$ is not an energy tolerance in Hartree.

Define the screened coefficient absolute-sum $\Lambda=\sum_j|h_j|$. The actual trotter_steps_naive implementation, for order $p$, returns at least one and otherwise the ceiling of

For first order this is $2\Lambda^2t^2/\epsilon$. Some older prose omits that factor two; source tests pin the implemented formula.

For first-order commutator selection, define

The code uses $N=\max(1,\lceil C_1t^2/(2\epsilon)\rceil)$. For the standard first-order term splitting, the corresponding commutator estimate follows the familiar norm bound. It detects exact commutation and can therefore be much smaller than a coefficient-only estimate.

For second order, the utility explicitly sums surviving nested commutators:

and chooses $N=\max(1,\lceil\sqrt{C_2|t|^3/(12\epsilon)}\rceil)$. For higher supported even order $p$, it sums all nonvanishing nested commutators of length $p+1$, bounding each by $2^p\prod|h_j|$, then uses the ceiling of $(C_p|t|^{p+1}/\epsilon)^{1/p}$.

The work needed for these estimates matters: the first-order sum is quadratic in term count, the second-order sum cubic, and the generic higher-order helper explicitly enumerates $L^{p+1}$ index tuples. Choosing an advanced bound can itself dominate classical preprocessing. These are the code’s implemented planning formulas; an asymptotic higher-order theorem or a finite test set is not a universal certificate for every possible grouping, recursion length, threshold, and time. In particular, the helper does not track a full schedule-specific remainder for arbitrary hand-supplied partitions.

For the illustrative $H=X+Z$, $t=1$, and $\epsilon=0.1$, the code formulas give $\Lambda=2$, $C_1=2$, and $C_2=6$. Thus the first-order commutator and naive choices are 10 and 80 slices, and the second-order choices are 3 and 18. These are hand substitutions into the documented implementation formulas, not executed accuracy measurements. The error tests and evolution tests contain formula checks and small-matrix comparisons. A step count alone does not establish chemical accuracy.

Standalone bound helpers call real-coefficient accessors and should be given already validated Hermitian input. They are not substitutes for checking an arbitrary complex operator. Coefficient screening also means their target concerns a screened model, not every term originally present.

Repeating a unitary versus rescaling its time

Shared settings power=1 and power_strategy="repeat" determine how a requested power is represented. Under repeat, the builder constructs its approximation for time $t$ once and multiplies its repetition count by the power $m$. Under rescale, it constructs an approximation at time $mt$ instead. These agree for exact evolution but need not agree for a finite approximation or random sampled schedule.

Automatic step selection under repeat is based on the base time, not an automatically tightened total-error budget for all $m$ uses. If each base unitary differs by at most $\epsilon_0$, telescoping gives a bound $m\epsilon_0$ for repeated unitary uses. Long phase-estimation powers therefore matter. Under rescale, a positive target accuracy makes the helper see the longer time, potentially selecting many more divisions/samples.

The output scale is the effective time supplied to construction: base time for repeat, multiplied time for rescale. A raw phase from a powered circuit cannot be decoded by ignoring the power and modular phase branches. Likewise reorder_terms currently resets scale metadata, as noted in Lesson 17. Validate phase conversion alongside circuit schedules, not after an unexplained energy discrepancy. Use valid positive powers for these intended workflows; shared power resolution is not a comprehensive independent input validator.

qDRIFT: importance sampling the terms

qdrift replaces a deterministic sweep by random draws. It samples term $j$ with probability

and applies, for every draw, the factor

Large coefficients are selected more often rather than assigned proportionally larger angles in each draw. Every sampled angle has the same magnitude, with the coefficient sign retained. Expanding one random step to first order shows its mean generator equals the desired Hamiltonian divided among the $N$ steps. Finite random fluctuations and noncommuting higher-order terms remain.

Defaults are num_samples=100, target_accuracy=0.0, error_bound="campbell", weight_threshold=1e-12, seed=-1, merge_duplicate_terms=True, and commutation_type="general", plus the common time/power settings. A nonnegative seed initializes NumPy’s generator reproducibly; −1 uses nondeterministic initialization. The generator is recreated when building, not embedded as a quantum random sampler in the output. Executing the same constructed circuit repeatedly reuses that compiled sequence.

The qdrift_samples_campbell helper uses

A positive target activates it; num_samples remains a manual floor and must be positive. The default floor of 100 can therefore dominate a smaller automatic result. Sampling and the norm calculation screen the same coefficient list. The sampling helper additionally returns an empty list when the remaining norm is below 1e-14, a separate very-small-operator cutoff. Identity terms can participate in sampling too; their phases should not be dropped without considering control.

Duplicate fusion is exact within the sampled unitary: it collects consecutive mutually commuting runs and combines identical Pauli words by adding angles. Distinct commuting words remain separate, and it does not cross noncommuting boundaries. commutation_type="qubit_wise" is a stricter alternative to general Pauli commutation. This optimization changes the emitted rotation count, not the number of random draws used in the sample-count prescription.

The crucial qualification is the error metric. The Campbell prescription concerns the averaged quantum channel: average the density-matrix transformations $\rho\mapsto U_\omega\rho U_\omega^\dagger$ over independently sampled random circuits $\omega$. A channel describes how states and their probabilities evolve; its diamond-distance metric also considers inputs entangled with a reference system. A small error of that averaged channel is not the same claim as small coherent operator error for one fixed sampled circuit.

The qDRIFT tests explicitly average channel representations over seeds and allow a finite-sampling tolerance. They do not assert every sampled circuit has operator error at most the requested epsilon. A single sample can fluctuate on a scale resembling $\Lambda|t|/\sqrt N$, rather than the averaged-channel $\Lambda^2t^2/N$ scale. Controlling and powering one fixed sample for QPE does not magically invoke the averaged-channel theorem. Energy phases, repeat/rescale choices, and randomness across circuit builds need a separately justified algorithmic analysis.

For $H=X+Z$, the sampling probabilities are one-half each. With $N$ draws, each angle is $2t/N$; the particular count and order of X and Z draws vary. This toy example shows both why the mean generator is right and why an individual schedule is not exact. A fixed seed makes experiments reproducible, not scientifically exact.

Partially randomized evolution: a heavy head and a light tail

partially_randomized divides the screened, magnitude-sorted terms into a deterministic head $H_D$ and a random tail $H_R$. Its weight_threshold has a different meaning from Trotter/qDRIFT: it is the head/tail split, not the removal tolerance. Terms at or above a nonnegative split threshold are deterministic. The separate tolerance=1e-12 screens small terms away.

With default split value −1, the implementation chooses the top ten percent by integer term count, at least one. For a single surviving term, it remains deterministic. Any negative split value enters this automatic branch in the implementation. The default is not a chemical orbital threshold or a fraction of cumulative coefficient weight.

Other defaults are trotter_order=2, num_random_samples=100, seed=-1, merge_duplicate_terms=True, and commutation_type="general". Only deterministic order 1 or 2 is supported. Unlike Trotter and qDRIFT, this settings class offers no target_accuracy auto-budgeting helper and no num_divisions time-slicing control.

Order 1 emits a full-time deterministic sweep followed by the tail samples. Order 2 emits a forward half-time deterministic sweep, the full-time randomized tail, and the reverse half-time deterministic sweep. Tail sampling uses $\Lambda_R=\sum_{j\in R}|h_j|$, not the full norm. Existing term partitions are not consumed by this method; magnitude sorting and its own sequence construction determine the order.

For a toy $H=X+0.1Z+0.05Y$, split threshold 0.5, time 0.2, and five tail draws, the second-order sequence begins and ends with an X factor of angle 0.1. The central random angles have magnitude $(0.1+0.05)0.2/5=0.006$. There are seven factors before any duplicate fusion. The tests check this head/tail structure and related norm/seed behavior; the numbers here are explanatory arithmetic.

This can be attractive when a few large terms dominate and a long tail has small total weight. But increasing tail samples does not remove the fixed deterministic splitting error of the outer sweep. Nor does a second-order outer palindrome turn one sampled tail into a deterministic second-order approximation with a universal operator-norm bound. The cited partially randomized phase-estimation research motivates the method; the current adapter does not implement a complete theorem-driven choice of split, time subdivisions, and full QPE error budget. Do not infer such a guarantee from a complexity formula in its docstring.

Zassenhaus: correcting the noncommuting product

Zassenhaus takes a different route: begin with exponentials of the pieces, then add exponentials of commutators that correct the mismatch. For generators $A$ and $B$, a formal product is

“Nesting” means taking the commutator of an operator with an already formed commutator. These higher terms account for order-sensitive interactions between the pieces. For quantum evolution the primitive generators are $-iH_j$; their factors of $-i$ must be retained when converting a correction to a real Pauli rotation angle.

For the toy choice $A=-iX$, $B=-iZ$, the second correction is $C_2=-iY$. The written product through that correction is $e^{-itX}e^{-itZ}e^{-it^2Y}$. Its chronological circuit applies the Y correction first, then Z, then X. That is not the same chronological order as the first-order example earlier. The current zassenhaus builder explicitly reverses its assembled factor list to preserve the required multiplication convention.

Default order is 2. Other defaults mirror Trotter’s time, power, division, accuracy, error-bound, and coefficient-threshold settings, with an additional term_grouper=AlgorithmRef("term_grouper","commuting") used when decomposing noncommuting correction exponents. Order 1 delegates to first-order Trotter. The current generic code accepts integer orders at least 2 through symbolic generation; it is not hard-capped at order 4 despite the older module summary. The inspected numerical convergence tests principally cover orders 2–4, while symbolic tests include degree 5. Higher acceptance is not a practical scaling or comprehensive validation guarantee. There is also no live automatic optimal-order search in the inspected execution path despite broad introductory prose suggesting one.

utils/zassenhaus_generation.py is an in-memory symbolic generator, not a file-generation step we need to run. It uses exact rational Fraction coefficients, antisymmetry to canonicalize commutators, and a dependency-ordered directed acyclic graph of shared subexpressions. A DAG is a graph whose arrows have no cycles, so children can be evaluated before parents. An LRU cache of up to 128 symbolic plans reuses structures for the same leaves/order; coefficient- and time-dependent Pauli calculations are separate.

The builder evaluates nested commutators with Pauli arithmetic, skips ones whose constituent operators all commute, and caches repeated evaluations. For a degree-$n$ correction it multiplies the symbolic coefficient by $i(-i)^n$ to express the exponent as $-i$ times a Hermitian Pauli operator, then uses time $t^n$. If a correction itself contains noncommuting terms, it selects an internal order $\lceil(p+1)/n\rceil-1$ for requested outer order $p$, recursively decomposing where necessary. Merely listing its terms once would not always preserve the desired overall order.

A final optimization merges equal Pauli rotations only by moving past commuting neighbors, with a final angle cutoff 1e-15. Correction screening at weight_threshold and this angle cutoff are additional numerical approximations. Nested commutators can create many new words and grow rapidly with order and group count; cached symbolic dependencies do not make the underlying Pauli expansion or subgroup combinations free.

Zassenhaus error estimates are leading-remainder estimates

The zassenhaus_error helpers inspect the first omitted exponent $C_{p+1}$. The commutator-aware route combines its actual Pauli coefficients after cancellation and uses their absolute-sum as an upper bound $\beta$ on that exponent’s norm. It selects

The naive route replaces the nested operator norms by coefficient-norm estimates and a symbolic coefficient absolute-sum $\kappa_{p+1}$, giving the analogous expression with $\beta=\kappa_{p+1}2^p\Lambda^{p+1}$. A provided symbolic-plan dictionary must contain the omitted degree; otherwise it raises. The public helpers reject nonpositive accuracy or order and incompatible containers.

These formulas capture the leading local omitted exponent and its expected small-step scaling. They do not sum a rigorously bounded infinite remainder at arbitrary time, fully certify each internal correction split, or automatically account for a different grouped leaf decomposition in every use. A zero leading estimate can mean that leading term vanishes, not a generic proof that all higher terms and implementation truncations vanish. Treat the result as a source-defined planning estimate requiring convergence/error validation, not a universal all-orders chemical guarantee.

The Zassenhaus tests verify symbolic coefficients, small-matrix convergence slopes near the expected order, optimizations, and particular model/molecular fixtures. Those tests support the implemented mechanism. We have read them, not rerun them or promoted their finite cases into an unrestricted theorem.

Choosing a method and an error budget

For a compact unexecuted use pattern:

from qdk_chemistry.algorithms import create

builder = create(
    "hamiltonian_unitary_builder", "trotter",
    time=0.2, order=2,
    num_divisions=4, target_accuracy=1e-3,
    error_bound="commutator", weight_threshold=1e-12,
    power=1, power_strategy="repeat",
)
plan = builder.run(qubit_hamiltonian)
circuit = create("circuit_mapper", "pauli_sequence").run(plan)

Here qubit_hamiltonian must be the intended Hermitian Pauli operator in the same encoding as the trial state. Four is a manual minimum, not necessarily the final slice count. The numerical choices illustrate API semantics; they are not a recommended nitrogen production budget or an executed result.

Keep errors in compatible units. Chemical basis/active-space/DFT-versus-wavefunction choices concern the physical model. Integral factorization and Pauli screening change the represented Hamiltonian. If that change is bounded by $\|\Delta H\|$, exact evolutions can differ by up to $|t|\|\Delta H\|$. Product-formula and rotation-synthesis errors instead concern the implemented unitary. Trial-state infidelity and finite-shot statistical error are other contributions, not automatically covered by target_accuracy.

For normalized states, a unitary-error bound $\epsilon$ bounds the resulting state-vector error; an observable expectation can then change by up to roughly $2\|O\|\epsilon$. Phase-to-energy conversion additionally involves time, modular branches, and eigenstate identification. Do not add a dimensionless circuit epsilon directly to a Hartree chemical error or label a simulation norm bound an energy confidence interval.

A sensible engineering comparison holds the represented Hamiltonian, target state/sector, total algorithmic objective, and error metric fixed, then counts complete controlled operations, repeated state loading, scratch, and synthesis precision. The longest QPE powers can dominate that accounting. For random methods, specify whether randomness is frozen in one compiled unitary, redrawn across jobs, or part of a justified partially randomized protocol.

Time evolution is therefore a method with explicit contracts, not a single black-box “simulate” button. The next lesson studies another route to exposing Hamiltonian information: place a normalized operator inside a larger unitary and use quantum-walk phases rather than a direct approximate time exponential.