23Following systems that change in time

Changing the Hamiltonian is a different scientific question

So far our nitrogen story has primarily asked about a fixed electronic model at a fixed geometry: neutral singlet N₂ at 1.85 Å, with a selected active space drawn from the 28-spatial-orbital cc-pVDZ model. A geometry scan would solve a sequence of different static problems. A classical geometry optimizer would use energy derivatives to move nuclei. Neither is automatically a real-time quantum simulation of an electronically driven molecule.

Here the question changes: given an initial state, how does it respond when the Hamiltonian itself varies during evolution? An external driving field is one possible motivation. But a scalar time-dependent coefficient is not, by itself, a validated laser–molecule model, solvent model, thermal bath, or nuclear trajectory. The field-coupling operator, gauge, units, approximations, and neglected degrees of freedom have to be specified scientifically.

The repository currently offers a concrete driven form, an order-1 Magnus time-averaging propagator, a stepwise circuit builder, and an evolution-plus-measurement orchestrator. We will trace those exact capabilities. The supplied demonstration is a small driven qubit model, not a newly computed nitrogen reaction or spectroscopic rate. No simulations, quadratures, imports, or tests are executed here.

A driven operator and the meaning of its parameters

The concrete supported form is

$H_0$ is the fixed base operator. $H_1$ is another fixed operator whose strength is modulated by a scalar function $f$. The Pauli support and register size remain fixed; coefficients change with time. “Driven” refers to this explicit modulation, not to a state being automatically kept in the instantaneous ground state.

For example, $f(t)=A\sin(\omega t+\phi)$ has amplitude A, angular frequency $\omega$, and phase $\phi$. These are values the user puts in the supplied function, not built-in DrivenQubitHamiltonian settings. If both operators carry energy units, f is usually dimensionless; then $\omega$ has inverse-time units and $\phi$ is an angle in radians. Alternatively units can be split between f and the coupling operator, but the product must still have energy units. The library cannot infer field strength units from a Python lambda.

A real f and Hermitian $H_0,H_1$ make $H(t)$ Hermitian. A complex modulation is not generally a valid Hermitian drive unless paired with the appropriate conjugate operator structure, which this simple two-operator form does not invent automatically.

DrivenQubitHamiltonian(h0,h1,drive) is a convenience wrapper around DrivenContainer. Construction checks that the two operators have equal qubit counts. evaluate(t) returns h0 + drive(t) * h1; properties expose h0, h1, and the function. It does not solve the Schrödinger equation or prepare a state.

The underlying QubitOperator addition checks matching widths, encoding strings, mode-order metadata, and tapering specifications. It concatenates Pauli terms rather than combining every duplicate. Matching partition types can be combined with shifted term indices; missing partition data removes that grouping, and incompatible partition types can raise. Zero drive leaves zero-weight H1 entries in the list rather than necessarily deleting them. The driven-data tests check these evaluation and partition behaviors.

The general TimeDependentQubitHamiltonian wrapper delegates to a container with evaluate(t) and num_qubits. These current wrappers are not the serialized DataClass containers used for many static chemistry objects. They provide no general JSON/HDF5 round-trip for an arbitrary Python callable. Keep the drive’s definition, parameters, units, and operator provenance separately; a closure is not a portable scientific specification merely because it can be evaluated locally.

Why the order of times matters

For a changing Hamiltonian, the evolution from $t_0$ to T is a time-ordered exponential:

The symbol $\mathcal T$ means later-time operations stand to the left in matrix multiplication, so they act after earlier ones. In general, $H(s_1)$ and $H(s_2)$ need not commute, and replacing the entire process by the exponential of the average Hamiltonian loses information about chronology.

For the driven form,

If the base and drive operators commute, time ordering is unnecessary and the scalar time integral is enough. If the drive is constant, the Hamiltonian is also constant, even if its separate H0 and H1 pieces do not commute. Their subsequent product-formula splitting can still introduce a different approximation.

A full-period sine drive can average to zero while producing a nontrivial evolution when $[H_0,H_1]\ne0$. Taking one large average step and concluding that the drive “did nothing” would then be wrong. That is a concrete reason to choose a time resolution related to the changing dynamics rather than just the final simulation duration.

Magnus expansion: one exponential with commutator corrections

The Magnus expansion writes the time-ordered evolution as $U=\exp\Omega$, with $\Omega=\Omega_1+\Omega_2+\cdots$. Using $A(s)=-iH(s)$, the first terms are

The minus sign comes from the two factors of $-i$. The commutator of two Hermitian operators is anti-Hermitian, so these terms can form the anti-Hermitian exponent required for unitary evolution. Higher terms contain further nested commutators and time integrals.

As a hand-worked toy, choose consistent dimensionless units and $H(s)=Z+sX$ on the interval [0,δ]. Then

The second term follows from $[Z,X]=2iY$ and the ordered integral of $s_2-s_1$. Dividing $i\Omega$ by δ would give an effective Hamiltonian $Z+\delta X/2-\delta^2Y/6+\cdots$. The induced Y correction illustrates information that a scalar average alone misses. It is a theoretical calculation, not a supported order-2 Magnus result produced by the current implementation.

For sufficiently smooth H, $H(s_1)-H(s_2)$ is small within a short interval. The leading omitted commutator integral after $\Omega_1$ scales locally as $\delta^3$; over a fixed total duration that leads to second-order global scaling in the step size when exact averaged-step exponentials are used. Discontinuous drives or uncontrolled quadrature/splitting errors do not automatically satisfy that conclusion.

What the native Magnus propagator actually returns

The registered/default method is propagator/magnus. Its Propagator interface takes a time-dependent wrapper and interval endpoints and returns a QubitOperator, not an evolved state, circuit, or executable unitary.

The concrete MagnusPropagator implements only order 1 and only DrivenContainer. It computes

A later builder must apply $e^{-i\Delta tH_{\mathrm{eff}}}$. The interval factor has not already been multiplied into the returned operator, so multiplying time twice would be an error.

The sole setting is order=1. Requesting order greater than one raises NotImplementedError; a settings object accepting order=2 is not evidence that it can run that method. The current guard does not reject every nonpositive order value, which would follow the same leading-order path rather than name a new algorithm. Use the implemented value 1. A different time-dependent container is rejected despite satisfying the abstract evaluate interface.

The drive integral is evaluated by external scipy.integrate.quad. No custom quadrature tolerances or sampling grid are exposed by this adapter. It takes the first return value and discards the quadrature error estimate. Therefore order=1 is not a certificate that the scalar integral was resolved to a requested chemistry accuracy. Rapid oscillation, discontinuities, singularities, or non-finite function values require scientific/numerical care beyond the wrapper’s settings.

The function is called adaptively at quadrature points, not merely at the beginning or midpoint of the interval. It should be a deterministic, finite function of its argument; a stateful function that changes after each call would not define the intended H(t). For a linear drive the average equals the midpoint value, but that is a property of that function, not a midpoint-sampling implementation.

The propagator tests check constant, zero, sine, linear, and step-function averages, subintervals, interval validation, and rejection of higher order. For example, the average of sin(s) over [0,π] is 2/π. These are source-defined reference checks, not quadrature runs performed here.

The Euler-named builder is a staged exponential integrator

The registered/default circuit builder is evolution_circuit_builder/euler. Despite the name, its default propagation is interval averaging, not the literal nonunitary numerical update $I-i\Delta tH$, nor a fixed left-endpoint sampler. It builds a sequence of unitary approximations to averaged intervals.

The settings default to total_time=1.0, dt=0.0, nested propagator/magnus, nested hamiltonian_unitary_builder/trotter, and circuit_mapper/pauli_sequence. The zero dt default must be replaced before running. Nested Trotter retains its own defaults, including first order and one effective division when no positive division/accuracy target is supplied.

The implementation starts at time zero, splits [0,T] into intervals of length dt, obtains one effective Hamiltonian per interval, sets the child time-evolution builder’s time to that interval’s duration, and concatenates its product-formula containers. The parent therefore overrides the child’s configured time. It requires an actual Propagator and TimeEvolutionBuilder, not an arbitrary algorithm reference of another family.

Only after constructing the combined mathematical evolution does it call the circuit mapper and prepend state preparation. This differs from a description that implies each interval is immediately compiled and executed. The ordinary mapping is lazy Q# circuit construction; execution is another layer.

There are two independent discretization controls. dt controls how well the changing Hamiltonian is represented by successive effective operators. The child Trotter num_divisions, order, and optional accuracy setting control exponentiating each fixed effective operator. Increasing Trotter divisions does not recover omitted Magnus time-ordering terms across a coarse dt interval. Decreasing dt does not by itself guarantee a chosen rotation-synthesis precision or eliminate shot noise.

Even though smooth Magnus order 1 has second-order global behavior when averaged exponentials are exact, the default first-order Trotter split can introduce a leading first-order global splitting contribution at fixed inner subdivision count. A blanket “the entire default workflow is second order” claim would therefore overstate the implementation. Choose and test the nested approximation as part of the method.

Interval validation and the final partial-step limit

The Euler builder rejects zero total time, zero dt, a step larger in magnitude than the same-sign total duration, and opposite signs. Its local checks allow matched negative time and dt. But the default Magnus propagator requires t_end > t_start, so the default complete pipeline does not support those backward intervals. Backward unitary evolution exists mathematically as an adjoint; that does not establish a working negative-time drive pipeline here.

The builder computes num_full_steps=int(total_time/dt) and a residual duration, and attempts to append a shorter cleanup step whenever the residual is nonzero. There is a directly coupled current limit: PauliProductFormulaContainer.combine requires matching scale values. The time-evolution builder puts each interval duration into that scale. A genuinely shorter residual generally has a different scale and the combination raises rather than completing the advertised cleanup path.

For example, requested T=0.25 and dt=0.1 would create two 0.1 slices and a 0.05 remainder; those scale values do not match. This is a source-derived interface consequence, not a failure we executed. The inspected tests use compatible equal-step or single-step cases and do not establish general shorter-final-step support. Floating-point residuals also deserve care because the branch checks exact nonzero residual before the approximate scale comparison.

For source-backed examples at this revision, choose a uniform grid that exactly reaches the desired endpoint, preferably with exactly representable values when practical, such as T=1.0 and dt=0.125. Do not silently claim arbitrary endpoint cleanup works. Nor should you relabel scales without understanding whether later consumers interpret them as physical time.

The combined container flattens step repetitions and merges adjacent identical Pauli rotations, preserving chronological order. It retains a compatible step scale, not a complete history of all time intervals. The final time-dependent product need not equal the exponential of one static Hamiltonian with that scale; using its ordinary phase-to-energy decoder to read a molecular static eigenenergy would therefore be inappropriate.

Composition and output capabilities

Prepending state preparation requires Q# callables on both the preparation and evolution circuits. The helper rejects an explicit encoding mismatch only when both encoding strings are present. The ordinary Pauli mapper often returns no encoding label, so a successful composition is not proof of matched orbitals, mode ordering, or a tapered representation.

It uses the Hamiltonian’s system width as its target register. It has no general machinery to provide a caller-owned shared phase gradient or transform an arbitrary wider preparation. The composed Circuit carries a sequential Q# factory/callable and a chosen encoding string, but not an explicit num_qubits declaration or copied shared-gradient metadata in this implementation.

The composition tests check ordering, missing-callable errors, endpoint/dt validation, and small executable/countable cases. The usage example demonstrates building a circuit and passing it to a resource-estimation adapter separately. Those example operations were read, not invoked here.

This builder is appropriate when one wants a final-time evolution circuit without automatically collecting statistics. It does not maintain a live quantum state between Python calls, report a time series by itself, or guarantee every possible representation conversion.

EulerIntegrator measures a final state, one observable at a time

The registered/default orchestration method is hamiltonian_simulation/euler_integrator. Its base settings contain three AlgorithmRefs: evolution_circuit_builder/euler, circuit_executor/qdk_sparse_state_simulator, and observable_estimator/expectation_estimator qdk by the corresponding family/name convention.

The actual run signature is run(hamiltonian, observables, state_prep, shots=1000, noise=None), with noise keyword-only. The EulerIntegrator implementation checks that every observable has the Hamiltonian’s qubit width, builds one complete state-preparation-plus-evolution circuit, then invokes the configured expectation estimator separately for each observable.

The parent’s circuit-executor reference is explicitly written into the nested estimator, overriding that child’s prior executor selection. shots is passed as that estimator’s total_shots for each observable, not divided across the entire observable list. Within each observable, the QDK estimator still divides its shot request by nonidentity measurement groups with floor rounding. Unsupported flat groups and ignored layered partitions behave exactly as taught in Lesson 20.

The return is a list of (EnergyExpectationResult, MeasurementData) pairs, in the observable list’s order. It is not a list indexed by simulation time. An observable called Z returns a dimensionless Z expectation despite the data field name energy_expectation_value. Passing H0 measures the base part of energy; passing hamiltonian.evaluate(T) measures the instantaneous represented Hamiltonian at T. The orchestrator does not choose that observable automatically.

A fixed scalar energy offset does not affect closed-system expectation dynamics, but belongs in a total-energy observable when comparing physical energies. A time-dependent scalar shift also produces only an overall state phase in unconditioned dynamics; its integral matters under controlled interference. The wrapper has no molecular-core bookkeeping field, so retain and apply constants consistently instead of reusing SOS or ordinary-Pauli offset rules without checking the representation.

A time series requires repeated prefixes

To observe multiple times $T_1,T_2,\ldots$, the present interface must be called with appropriately configured final-time circuits for those endpoints. Each circuit begins at zero and includes the full prefix evolution from the initial preparation. There is no hidden state object advanced in place from one measured time to the next.

This is scientifically sensible for ordinary measurements: measuring at an intermediate time can disturb the state whose unmeasured evolution one wanted to study. Separate runs at each observation time estimate the unmeasured trajectory’s observable values without that intervention, assuming reproducible preparation and controls.

If the j-th endpoint needs j equal time steps, M endpoints require a total of $1+2+\cdots+M=M(M+1)/2$ step uses per corresponding shot schedule, before within-step Pauli expansion. With S actual shots at each endpoint for one observable, this becomes $S M(M+1)/2$ step executions rather than SM. Multiple separately measured observables add their own shot costs. These are illustrative workload relations, not measured timings.

The same prefixes are also reconstructed classically by separate builder calls unless a caller implements its own reuse strategy; the current integrator does not expose a general time-series cache. Because algorithm settings lock after ordinary run, configure fresh builders/references for different endpoint settings rather than assuming a completed algorithm can always be mutated. A T=0 baseline is obtained from the initial-state measurement path; the Euler circuit builder itself rejects zero total time.

Returned result pairs do not automatically carry their time, drive definition, all approximation settings, or the original molecular provenance. A caller constructing a series must retain those labels alongside each pair. The energy_variance field remains the QDK estimator’s diagonal plug-in variance of term means, omitting within-group covariance; it is neither intrinsic energy variance nor a certified error bar for the time trajectory. MeasurementData has the same allocated/usable-shot and serialization-provenance limitations from Lesson 20.

The notebook’s concrete model and reference calculation

The source notebook time_evolve_and_measure.ipynb uses a two-site Ising-style qubit model. Its base operator is a ZZ coupling, its driven operator a transverse X-field sum, and the drive is sin(2πt). Declared parameters are total time 2.0, dt=0.1, two inner first-order Trotter subdivisions, and an initial all-zero state from the identity helper. It measures ZZ at the final time, comparing reference, noiseless, and depolarizing-noise simulation paths.

This is an endpoint demonstration, not a stored multi-time trajectory and not a laser-driven N₂ calculation. Its initial all-zero state is appropriate to the toy spin model; it would not automatically represent our six-electron chemical trial state. Lattice-model semantics will be explored in the next lesson.

The notebook’s reference uses SciPy quadrature for each coarse interval, then expm_multiply repeatedly on the same averaged H_eff inside that interval. The fine internal substeps do not resample the drive or add Magnus commutator terms. In exact arithmetic, repeated exponentials of that same constant H_eff equal one exponential for the coarse interval. Therefore this reference is useful for testing product-formula and sampling error relative to the averaged-step model, but it is not independently exact continuous-time evolution merely because it has many internal substeps.

The source prints numeric results only when executed; none are being asserted here. A damped noisy ZZ signal in such a model would not, on its own, establish a molecular relaxation rate, heating law, solvent response, or chemical reaction rate. Those interpretations need an appropriate physical model and calibrated noise/environment assumptions.

Consistency checks and an unexecuted pattern

A helpful analytic limit is $H(t)=f(t)X$. All Hamiltonians at different times commute, so

for initial $|0\rangle$. A constant f=1 gives cos(2T), and a linear f(s)=s gives cos(T²). The integrator tests use such cases to catch missing or doubled duration factors. Although Z does not commute with X, the time-dependent Hamiltonians in this family commute with one another; those tests do not validate arbitrary noncommuting time ordering.

For a zero Hamiltonian, the correct conceptual evolution is identity, with unchanged observables. A constant Hamiltonian should recover the fixed-H methods of Lesson 18, subject to the chosen product-formula error. For a genuinely driven noncommuting model, convergence should be checked by refining the outer intervals and inner approximation separately, not by observing that one easy commuting limit passed.

Here is a small unexecuted pattern with an exactly representable uniform endpoint grid:

from math import pi, sin
import numpy as np
from qdk_chemistry.algorithms import create
from qdk_chemistry.algorithms.state_preparation import identity_state_prep
from qdk_chemistry.data import AlgorithmRef, DrivenQubitHamiltonian, QubitOperator

h0 = QubitOperator(["Z"], np.array([1.0]))
h1 = QubitOperator(["X"], np.array([0.5]))
driven = DrivenQubitHamiltonian(h0, h1, drive=lambda t: sin(2 * pi * t))
initial = identity_state_prep(1)
simulation = create(
    "hamiltonian_simulation", "euler_integrator",
    evolution_circuit_builder=AlgorithmRef(
        "evolution_circuit_builder", "euler", total_time=1.0, dt=0.125,
        propagator=AlgorithmRef("propagator", "magnus", order=1),
        evolution_builder=AlgorithmRef(
            "hamiltonian_unitary_builder", "trotter", order=2, num_divisions=2,
            power=1, power_strategy="repeat",
        ),
    ),
    circuit_executor=AlgorithmRef("circuit_executor", "qdk_sparse_state_simulator", seed=42),
)
# Authorized execution would call:
# results = simulation.run(driven, [QubitOperator(["Z"], np.array([1.0]))], initial, shots=1000)

The energies/times here are consistent illustrative model units, not nitrogen parameters or a converged accuracy recommendation. H0 and H1 do not commute, so this example genuinely exercises the purpose of time slicing. It requires the QDK runtime when eventually executed; no object construction or imports from this example were run by the course writer.

Separating the error sources

There are at least four numerical layers before statistics: approximating the physical driven model, quadrature of f, truncation of Magnus/time ordering, and finite product-formula/rotation synthesis. Pauli coefficient cutoffs and container-combination cancellation thresholds add further choices. They are not all controlled by dt, order, or one nested target_accuracy.

A smooth Magnus-order-1 method may improve quadratically with outer dt when its assumptions hold; the complete circuit may converge differently if first-order splitting, quadrature error, or a coarse approximation dominates. Highly oscillatory and discontinuous drives require resolution and interval choices appropriate to their structure. The adapter does not expose adaptive step-size control, error-triggered refinement, or a general higher-order Magnus method.

Finite-shot error sits on top of that evolution approximation. Noise changes the measured distribution rather than only enlarging a clean statistical error bar. Loss postselection and ignored grouped covariance can make a small reported variance misleading. Increasing shots cannot correct a wrongly modeled drive, a faulty representation pairing, or a too-coarse time-ordering approximation.

Finally, simulated time is not hardware runtime. Hardware must implement every interval’s gates, every prefix, every shot and observable group, all ancillas, and eventually fault-tolerant synthesis/error correction if required. A statevector reference stores exponentially many amplitudes and is useful only for limited validation sizes, not an assumed scalable production comparison.

The achievement of this workflow is a traceable chain from an explicitly specified time-dependent qubit operator, through an interval-effective generator and circuit, to final observable statistics. It does not automatically model real nuclear motion or environmental response. That precise boundary lets us extend the chemistry story responsibly and prepares us to examine the repository’s concrete lattice/model Hamiltonians next.