Noise and Open System Simulation

Noise and Open-System Simulation

Model discrete circuit noise, separate readout errors from quantum-state evolution, and simulate continuous-time Markovian open systems.

Scope: FlagQuantum 0.3.1. All six complete scripts below were verified on October 20, 2026 with Python 3.12.14, PyTorch 2.13.0, and macOS ARM64 using local CPU execution. Examples use synthetic parameters, not measured device calibration. No GPU, distributed runtime, provider account, or quantum hardware is required or validated here.

Install the core package using Installation and Compatibility. Familiarity with Core Concepts, Building and Simulating Circuits, and Choosing a Simulation Method is helpful.

Formula/code review — October 20, 2026: The equations are paired with code excerpts, and the local CPU examples were re-run with Python 3.12.14 and PyTorch 2.13.0. See the review record for executed commands, distributed checks, and untested hardware/provider paths. Earlier environment and output records are retained.

1. Choose a Noise Workflow

Workflow	Interface	What it computes	Main limitation
Discrete Kraus-channel evolution	fq.run(..., noise_model=..., options=ExecutionOptions(mode="density_matrix"))	Density matrix and requested exact-channel probabilities/expectations for a circuit	Dense state memory grows as 4**n.
Explicit inline channels	Four channel methods on fq.Circuit	Noise at a specific position in a circuit	Only four named scalar channel opcodes are exposed this way.
Statevector trajectories	flagquantum.runtime.run_noisy_statevector(...)	Statistics over sampled pure-state trajectories, processed in trajectory batches	Exponential statevector memory and finite-trajectory error.
MPS trajectories	flagquantum.runtime.run_noisy_mps(...)	Statistics over sampled MPS trajectories	Entanglement growth, truncation, and dense fallback for multi-qubit channels.
Continuous-time evolution	flagquantum.lindblad.plan(...) and .run(...)	Time-resolved density matrices, populations, and expectations from a Lindblad generator	CPU-only dense fixed-step integration, with a time-independent Hamiltonian.
Classical readout confusion	NoiseModel.add_readout(...), correlated readout records, or explicit probability transformation	A model of observed classical bits	Separate from the quantum density matrix and dependent on the chosen output path.

These methods answer different questions. A circuit channel acts at an explicit program location; continuous-time evolution integrates a Hamiltonian and dissipative generator over a time grid. A trajectory is one sampled noise history, while a measurement shot is one sampled readout outcome. Increasing shots does not remove finite-trajectory error.

For a Kraus channel,

$$\rho' = \sum_k K_k\rho K_k^\dagger,\qquad \sum_k K_k^\dagger K_k = I.$$

Formula / mathematical operation	Corresponding code excerpt	Check / interpretation
$\rho\mapsto\sum_k K_k\rho K_k^\dagger$	fq.run(circuit, options=options, outputs=[fq.density_matrix(), fq.probabilities()])	Set options.mode to density_matrix; inline channels are exercised in noise_channels.py.
$\sum_k K_k^\dagger K_k=I$	channel.kraus	Channel factories expose the Kraus operators and validate completeness; the script below also checks an explicit damping reference.
$\rho=\rho^\dagger,\quad\mathrm{Tr}(\rho)=1$	torch.allclose(rho, rho.mH, atol=1e-6)	The script checks Hermiticity, trace, populations, expectations, and nonnegative eigenvalues.

“Exact density-matrix evolution” means applying the full channel instead of sampling noise histories. It still has floating-point error, and its accuracy as a hardware model depends on the supplied assumptions and parameters.

2. Supported Models and Parameter Conventions

Factories and records live under flagquantum.noise. A NoiseModel attaches channel rules after matching gates; classical readout rules are stored separately.

Model / factory	Parameters and meaning
bit_flip_channel(probability)	Apply X with probability p, otherwise identity.
phase_flip_channel(probability)	Apply Z with probability p, otherwise identity.
depolarizing_channel(probability)	Identity with probability 1-p; X, Y, and Z each with probability p/3.
two_qubit_depolarizing_channel(probability)	Identity with probability 1-p; each of the other 15 two-qubit Paulis with probability p/15.
amplitude_damping_channel(gamma)	Excited-state decay toward `
phase_damping_channel(gamma)	Preserve populations and multiply off-diagonal coherence by sqrt(1-gamma).
reset_error_channel(probability_zero, probability_one=0)	Reset to zero with probability p0, reset to one with probability p1, otherwise leave the state; require p0+p1 <= 1.
coherent_overrotation_channel(angle, axis="x")	A deterministic extra rotation in radians around X, Y, or Z, represented by one Kraus operator.
thermal_relaxation_channel(t1, t2, duration, excited_population=0)	Relaxation/dephasing over a duration, including an optional thermal equilibrium excited population.
ReadoutError / CorrelatedReadoutError	Row-stochastic true-to-observed classical confusion matrices for one or multiple qubits.
DeviceNoiseProfile	Timestamped T1/T2/readout calibration and gate durations used to lower timing and idle relaxation.

Channel probabilities are real scalars in [0, 1]. Channel construction validates finite, compatible Kraus matrices and the completeness condition. A valid channel definition does not establish that every execution representation handles its arity efficiently.

The depolarizing parameter above is the total nonidentity Pauli probability. It is not automatically the same parameter as another library's convention rho -> (1-p)rho + p I/2. Compare the equations before transferring a noise strength between frameworks.

Only bit_flip, phase_flip, depolarizing, and amplitude_damping are direct circuit methods in this release. Use model rules for the other built-in channels. Channel probabilities are bound model values; trainable noisy-channel parameters and general noisy gradients are outside the documented support contract.

3. Run Four Basic Channels with Density Matrices

Save as noise_channels.py:

import torch
import flagquantum as fq

options = fq.ExecutionOptions(mode='density_matrix', backend='pytorch', device='cpu', precision='complex128')
cases = [
    ('bit_flip', fq.Circuit(1).bit_flip(0, 0.25), [0.75, 0.25], 0.0, 0.5),
    ('phase_flip', fq.Circuit(1).h(0).phase_flip(0, 0.25), [0.5, 0.5], 0.5, 0.0),
    ('depolarizing', fq.Circuit(1).depolarizing(0, 0.3), [0.8, 0.2], 0.0, 0.6),
    ('amplitude_damping', fq.Circuit(1).x(0).amplitude_damping(0, 0.25), [0.25, 0.75], 0.0, -0.5),
]
for name, circuit, populations, expected_x, expected_z in cases:
    result = fq.run(circuit, options=options, outputs=[
        fq.density_matrix(), fq.probabilities(),
        fq.expectation(fq.X(0), name='x'), fq.expectation(fq.Z(0), name='z'),
    ])
    rho = result.density_matrix[0]
    p = result.probabilities[0]
    x = result.expectation('x').item()
    z = result.expectation('z').item()
    assert torch.allclose(p, torch.tensor(populations, dtype=p.dtype), atol=1e-6)
    assert abs(x - expected_x) < 1e-6 and abs(z - expected_z) < 1e-6
    assert torch.allclose(rho, rho.mH, atol=1e-6)
    assert abs(torch.trace(rho).real.item() - 1) < 1e-6
    assert torch.linalg.eigvalsh(rho).min().item() >= -1e-6
    print(f'{name}: P={[round(v, 6) for v in p.tolist()]}, <X>={x:.6f}, <Z>={z:.6f}')
# Independently apply amplitude-damping Kraus matrices to |1><1|.
from flagquantum.noise import amplitude_damping_channel
channel = amplitude_damping_channel(0.25, dtype=torch.complex128)
rho_in = torch.diag(torch.tensor([0.0, 1.0], dtype=torch.complex128))
identity = torch.eye(2, dtype=torch.complex128)
completeness = sum(k.mH @ k for k in channel.kraus)
rho_reference = sum(k @ rho_in @ k.mH for k in channel.kraus)
assert torch.allclose(completeness, identity, atol=1e-12)
assert torch.allclose(rho, rho_reference, atol=1e-6)
print('Analytic channel and density-matrix checks passed.')

Run python noise_channels.py. Expected output:

bit_flip: P=[0.75, 0.25], <X>=0.000000, <Z>=0.500000
phase_flip: P=[0.5, 0.5], <X>=0.500000, <Z>=0.000000
depolarizing: P=[0.8, 0.2], <X>=0.000000, <Z>=0.600000
amplitude_damping: P=[0.25, 0.75], <X>=0.000000, <Z>=-0.500000
Analytic channel and density-matrix checks passed.

The phase-flip example starts in |+>. Its computational-basis populations stay [0.5, 0.5], but <X> decreases, exposing coherence loss that a population-only check would miss. The amplitude-damping example starts in |1> and transfers population to |0>.

result.density_matrix and result.probabilities are properties; result.expectation(name) is a method. A single circuit retains a leading batch dimension in this path. The script checks trace, Hermiticity, and density-matrix eigenvalues as well as the analytical observables.

The explicit density-matrix route applies inline channels without a separate NoiseModel. Pure-amplitude modes statevector, mps, and tensor_network do not become noisy-channel simulators merely by selecting that mode; use the exact-channel route or the dedicated trajectory APIs. Inline channels do not silently turn into noiseless gates.

4. Attach Noise Rules and Model Readout Separately

Save as noise_rules_readout.py:

import torch
import flagquantum as fq
from flagquantum.noise import NoiseModel, ReadoutError, amplitude_damping_channel

options = fq.ExecutionOptions(mode='density_matrix', backend='pytorch', device='cpu', precision='complex128')
circuit = fq.Circuit(1).x(0)
model = NoiseModel().add('x', amplitude_damping_channel(0.25, dtype=torch.complex128), qubits=(0,))
rule_result = fq.run(circuit, noise_model=model, options=options, outputs=fq.density_matrix())
inline_result = fq.run(fq.Circuit(1).x(0).amplitude_damping(0, 0.25), options=options, outputs=fq.density_matrix())
assert torch.allclose(rule_result.density_matrix, inline_result.density_matrix, atol=1e-6)
true_p = torch.diagonal(rule_result.density_matrix[0]).real
# Rows: true 0/1. Columns: observed 0/1.
readout_model = NoiseModel().add_readout(0, ReadoutError(((0.9, 0.1), (0.2, 0.8))))
observed_p = readout_model.apply_readout_probabilities(true_p, n_qubits=1)
assert torch.allclose(observed_p, torch.tensor([0.375, 0.625], dtype=true_p.dtype), atol=1e-6)
restored = NoiseModel.from_dict(model.to_dict())
assert restored.identity == model.identity
print('Gate rule matches inline channel: True')
print('Quantum-state probabilities:', [round(v, 6) for v in true_p.tolist()])
print('Observed probabilities:', [round(v, 6) for v in observed_p.tolist()])
print(f'State <Z>: {(true_p[0] - true_p[1]).item():.6f}')
print(f'Observed-bit <Z>: {(observed_p[0] - observed_p[1]).item():.6f}')
print('Noise model serialization identity preserved: True')

Run python noise_rules_readout.py. Expected output:

Gate rule matches inline channel: True
Quantum-state probabilities: [0.25, 0.75]
Observed probabilities: [0.375, 0.625]
State <Z>: -0.500000
Observed-bit <Z>: -0.250000
Noise model serialization identity preserved: True

The gate-bound rule and inline channel agree because both place the same damping channel once, after the same X gate. A rule applies after every matching occurrence: a rule on cx throughout a ladder is a different model from one channel inserted only at the end. Matching gate names and target qubits therefore matters.

A one-qubit channel attached to a multi-qubit gate can be placed on its selected wires; a genuine two-qubit channel is a distinct operation. Inspect the lowered instruction sequence with flagquantum.compiler.lower_noise_model when placement is unclear.

Readout matrix rows are true states and columns are observed states. Here the observed probability of zero is 0.25*0.9 + 0.75*0.2 = 0.375. The underlying quantum state still has probabilities [0.25, 0.75]; classical reporting errors do not change that density matrix.

This script deliberately transforms probabilities explicitly, making the boundary visible. Runtime measurement paths can also consume readout rules; the trajectory APIs apply configured readout rules to their per-qubit Z statistics. Do not apply the same confusion matrix again to an already transformed output, or assume every arbitrary observable includes a classical readout model automatically.

Models support versioned to_dict() / from_dict() serialization and a content identity. Record that identity along with circuit, versions, precision, and execution settings. NoiseModel itself is mutable; adding rules changes the model being simulated.

5. Use T1/T2 Relaxation and Timing Profiles

T1, T2, and duration must use the same time unit. Require positive T1/T2, nonnegative duration, and T2 <= 2*T1. At zero equilibrium excited population,

$$P_1(t)=P_1(0)e^{-t/T_1},\qquad \rho_{01}(t)=\rho_{01}(0)e^{-t/T_2}.$$

Formula / mathematical operation	Corresponding code excerpt	Check / interpretation
$P_1(t)=e^{-t/T_1}$	math.exp(-duration / t1)	noise_thermal.py starts in the excited state, so P1(0)=1.
$\langle X\rangle(t)=e^{-t/T_2}$	math.exp(-duration / t2)	For the initial plus state, real coherence gives this X expectation.
$\mathcal E_{T_1,T_2,t}$	thermal_relaxation_channel(t1, t2, duration, dtype=torch.complex128)	Apply after X or H, then compare both measured quantities with the exponentials, with tolerance 1e-6.

Save as noise_thermal.py:

import math
import torch
import flagquantum as fq
from flagquantum.noise import NoiseModel, thermal_relaxation_channel

options = fq.ExecutionOptions(mode='density_matrix', backend='pytorch', device='cpu', precision='complex128')
# All three times use the same unit, for example microseconds.
t1, t2, duration = 50.0, 70.0, 10.0
channel = thermal_relaxation_channel(t1, t2, duration, dtype=torch.complex128)
excited = fq.run(fq.Circuit(1).x(0), noise_model=NoiseModel().add('x', channel), options=options, outputs=fq.probabilities()).probabilities[0, 1].item()
x = fq.run(fq.Circuit(1).h(0), noise_model=NoiseModel().add('h', channel), options=options, outputs=fq.expectation(fq.X(0))).expectation().item()
assert abs(excited - math.exp(-duration / t1)) < 1e-6
assert abs(x - math.exp(-duration / t2)) < 1e-6
print(f'Excited population after relaxation: {excited:.6f}')
print(f'<X> after dephasing and relaxation: {x:.6f}')
try:
    thermal_relaxation_channel(t1=50.0, t2=110.0, duration=10.0)
except ValueError:
    print('Unphysical T2 > 2*T1 rejected: True')
else:
    raise AssertionError('Expected invalid relaxation times to be rejected')
print('Thermal analytic checks passed.')

Run python noise_thermal.py. Expected output:

Excited population after relaxation: 0.818731
<X> after dephasing and relaxation: 0.866878
Unphysical T2 > 2*T1 rejected: True
Thermal analytic checks passed.

This is a finite-duration channel applied at a circuit location. It does not jointly integrate an arbitrary drive Hamiltonian during that gate. Use the Lindblad workflow below when the intended model is simultaneous coherent driving and dissipative evolution.

Calibration-Driven Gate and Idle Noise

Save as noise_profile.py:

import torch
import flagquantum as fq
from flagquantum.compiler import lower_noise_model
from flagquantum.noise import DeviceNoiseProfile, GateDuration, NoiseModel, QubitNoiseCalibration

# Illustrative values, not a real device calibration.
profile = DeviceNoiseProfile(
    qubits=(
        QubitNoiseCalibration(0, t1=50_000.0, t2=70_000.0),
        QubitNoiseCalibration(1, t1=48_000.0, t2=65_000.0),
    ),
    gate_durations=(GateDuration('h', 35.0), GateDuration('cx', 280.0)),
    source='wiki-synthetic-profile', captured_at='2026-10-20T00:00:00+00:00', time_unit='ns',
)
noise = NoiseModel.from_device_profile(profile)
circuit = fq.Circuit(2).h(0).cx(0, 1)
lowered = lower_noise_model(circuit, noise)
channels = [gate for gate in lowered.instructions if gate.metadata.get('is_channel')]
assert len(channels) == 4
for gate in channels:
    assert gate.metadata['device_profile_identity'] == profile.identity
    print(f'{gate.metadata["placement"]}: qubit={gate.wires[0]}, duration={gate.metadata["duration"]:.1f} ns')
result = fq.run(circuit, noise_model=noise, outputs=fq.density_matrix(), options=fq.ExecutionOptions(mode='density_matrix', device='cpu'))
rho = result.density_matrix[0]
assert abs(torch.trace(rho).real.item() - 1) < 1e-5
assert torch.linalg.eigvalsh(rho).min().item() >= -1e-5
print('Profile channels and final density-matrix checks passed.')

Run python noise_profile.py. Expected output:

during_gate_approximation: qubit=0, duration=35.0 ns
idle_before_gate: qubit=1, duration=35.0 ns
during_gate_approximation: qubit=0, duration=280.0 ns
during_gate_approximation: qubit=1, duration=280.0 ns
Profile channels and final density-matrix checks passed.

The ASAP qubit clocks let qubit 1 idle while H acts on qubit 0, then synchronize them for CX. Lowering can insert relaxation during gate durations, before synchronization, and for terminal idle time. Generated instructions record placement, duration, time unit, source gate, and profile identity. Missing required gate durations or qubit calibration is rejected.

This example uses floating-point calibration times and durations consistently. In the tested 0.3.1 release, supplying these values as integers can produce a noise-model identity mismatch after serialization reconstructs them as floats. Use floating-point values when constructing the profile and verify a serialization round trip when saving it.

Profiles implement an explicit Markovian timing approximation. They do not model pulse overlap, crosstalk, leakage, or arbitrary correlated environment dynamics, and a timestamped profile alone is not evidence of agreement with hardware.

6. Compare Trajectories with an Exact Reference

Each trajectory samples Kraus branches and evolves a pure state. Averaging their observables estimates the full-channel result. The expert statevector and MPS APIs expose trajectory statistics rather than the same result envelope as fq.run.

Save as noise_trajectories.py:

import torch
import flagquantum as fq
from flagquantum.noise import NoiseModel, amplitude_damping_channel
from flagquantum.runtime import run_noisy_mps, run_noisy_statevector

circuit = fq.Circuit(1).x(0)
noise = NoiseModel().add('x', amplitude_damping_channel(0.25, dtype=torch.complex128))
exact = fq.run(circuit, noise_model=noise, outputs=fq.expectation(fq.Z(0)), options=fq.ExecutionOptions(mode='density_matrix', device='cpu', precision='complex128')).expectation().item()
mps = run_noisy_mps(circuit, noise, trajectories=512, seed=42, device='cpu', dtype=torch.complex128, retain_trajectories=False, max_bond=None, cutoff=0.0)
sv = run_noisy_statevector(circuit, noise, trajectories=512, trajectory_batch_size=32, seed=42, device='cpu', dtype=torch.complex128, retain_trajectories=False)
print(f'Exact <Z>: {exact:.6f}')
for name, mean, se, count in [
    ('MPS', mps.expectation_z_mean, mps.statistics.standard_error, mps.n_trajectories),
    ('Statevector', sv.expectation_z, sv.standard_error, sv.statistics.count),
]:
    value, error = mean.item(), se.item()
    assert count == 512 and error > 0
    # A smoke check, not a guaranteed confidence bound for arbitrary runs.
    assert abs(value - exact) < 6 * error
    print(f'{name}: trajectories={count}, <Z>={value:.6f}, standard_error={error:.6f}')
print('Both trajectory estimates agree with the exact reference within six reported standard errors.')

Run python noise_trajectories.py. Output observed in the verified environment:

Exact <Z>: -0.500000
MPS: trajectories=512, <Z>=-0.566406, standard_error=0.036422
Statevector: trajectories=512, <Z>=-0.566406, standard_error=0.036422
Both trajectory estimates agree with the exact reference within six reported standard errors.

The seed reproduces this test configuration. Exact stochastic values can vary across dependency versions or runtime implementations; agreement of both engines on this simple case is not a promise that every circuit produces identical trajectory samples. The assertions use reported uncertainty rather than requiring the mean to equal the reference exactly. Six standard errors is a smoke-test threshold, not a guaranteed confidence interval for arbitrary data.

Important distinctions:

• trajectories controls noise-history sampling. More trajectories generally reduce Monte Carlo standard error; this is separate from final measurement-shot noise.
• trajectory_batch_size controls the statevector trajectory work block, not a batch of independent input circuits. A block still stores exponentially many amplitudes per trajectory.
• MPS bond dimensions and cutoff control representation error. This one-qubit test needs no entanglement truncation and requests no bond cap. A low standard error cannot certify negligible MPS truncation bias.
• Single-qubit MPS channels sample branches in MPS form. Multi-qubit channels currently use a dense statevector correctness fallback before rebuilding the MPS, which can eliminate the expected memory advantage.
• retain_trajectories=False keeps aggregate statistics rather than retaining states for later shot sampling. Use the appropriate retention setting if you need that additional operation.

The MPS expert API supports single-rank adaptive stopping with min_trajectories and target_standard_error. Inspect completed counts, converged, stopped_early, and failures; reaching the requested ceiling does not prove a target was met. Distributed adaptive stopping and production multi-GPU scheduling are outside the documented noise support boundary.

The stable mode="auto" route can select noisy MPS when a density matrix exceeds an explicit memory budget and approximation is authorized with allow_approximate=True. Its documented policy uses 32 trajectories; increasing measurement shots does not increase that number. Use the expert APIs when trajectory count, stopping policy, bond dimension, or cutoff must be controlled explicitly. This page validates the small CPU expert examples, without making a general capacity or accelerator claim.

7. Simulate Continuous-Time Lindblad Evolution

For a time-independent Hamiltonian and Markovian collapse operators, the solver integrates

$$\frac{d\rho}{dt}=-i[H,\rho]+\sum_j\left( L_j\rho L_j^\dagger-\frac12\{L_j^\dagger L_j,\rho\}\right).$$

Formula / mathematical operation	Corresponding code excerpt	Check / interpretation
$\dot\rho=-i[H,\rho]+\mathcal D(\rho)$	result = fql.run(plan)	noise_lindblad.py builds a plan with a fixed Hamiltonian, time grid, and collapse operators.
$L=\sqrt\kappa\lvert0\rangle\langle1\rvert$	fql.amplitude_damping(rate=rate, qubit=0)	The collapse operator uses inverse-time rate kappa, not a per-gate probability.
$P_1(t)=e^{-\kappa t}$	expected_p1 = torch.exp(-rate * times)	The zero-drive case is checked against result.populations[:, 1], with atol=1e-8.
$\langle Z\rangle(t)=1-2P_1(t)$	1 - 2 * expected_p1	Compare against result.expectation('z'), with atol=1e-8. Hermiticity, positivity, and trace drift are checked too.

The helper fql.amplitude_damping(rate, qubit) represents L = sqrt(rate)|0><1|. Its rate has units of inverse time and is not a probability. With no drive and an initially excited qubit, P1(t)=exp(-rate*t). The equivalent finite-interval damping probability is gamma=1-exp(-rate*duration) for this decay-only model. Do not substitute a per-gate probability directly for a physical rate.

Save as noise_lindblad.py:

import torch
import flagquantum as fq
import flagquantum.lindblad as fql

rate = 0.7  # Inverse time; not a channel probability.
times = torch.linspace(0.0, 2.0, 201, dtype=torch.float64)
options = fq.ExecutionOptions(mode='density_matrix', device='cpu', precision='complex128')
plan = fql.plan(
    torch.zeros((2, 2), dtype=torch.complex128),
    '1', times,
    collapse_operators=[fql.amplitude_damping(rate=rate, qubit=0)],
    outputs=fq.expectation(fq.Z(0), name='z'),
    options=options,
    return_density_matrices=True,
)
result = fql.run(plan)
expected_p1 = torch.exp(-rate * times)
assert torch.allclose(result.populations[:, 1], expected_p1, atol=1e-8)
assert torch.allclose(result.expectation('z'), 1 - 2 * expected_p1, atol=1e-8)
assert result.density_matrices.shape == (201, 2, 2)
assert torch.allclose(result.density_matrices, result.density_matrices.mH, atol=1e-10)
assert torch.linalg.eigvalsh(result.density_matrices).min().item() >= -1e-8
assert result.maximum_trace_drift < 1e-10 and result.population_bounded
print('Time points:', plan.n_times)
print('Population shape:', tuple(result.populations.shape))
for i in (0, 50, 100, 200):
    print(f't={times[i].item():.1f}: P(1)={result.populations[i, 1].item():.6f}, <Z>={result.expectation("z")[i].item():.6f}')
print(f'Integrator: {result.method}, order={result.order}, step_size={result.step_size:.6f}')
print('Trace drift below 1e-10:', result.maximum_trace_drift < 1e-10)
print('Analytic decay, Hermiticity, and positivity checks passed.')

# The same facade supports a time-independent coherent drive plus damping.
driven = fql.run(
    0.5 * fq.X(0) + 0.1 * fq.Z(0), '1', times,
    collapse_operators=[fql.amplitude_damping(rate=rate, qubit=0)],
    outputs=fq.expectation(fq.Z(0), name='z'), options=options,
)
assert driven.population_bounded and driven.maximum_trace_drift < 1e-10
print('Driven final populations:', [round(p, 6) for p in driven.populations[-1].tolist()])
print(f'Driven final <Z>: {driven.expectation("z")[-1].item():.6f}')

Run python noise_lindblad.py. Expected output:

Time points: 201
Population shape: (201, 2)
t=0.0: P(1)=1.000000, <Z>=-1.000000
t=0.5: P(1)=0.704688, <Z>=-0.409376
t=1.0: P(1)=0.496585, <Z>=0.006829
t=2.0: P(1)=0.246597, <Z>=0.506806
Integrator: runge-kutta, order=4, step_size=0.005000
Trace drift below 1e-10: True
Analytic decay, Hermiticity, and positivity checks passed.
Driven final populations: [0.80357, 0.19643]
Driven final <Z>: 0.607140

The first run checks the entire population and expectation trajectory against the decay formula. The second uses a constant coherent drive H=0.5 X + 0.1 Z together with damping. Hamiltonian coefficients, rates, and times must use a consistent unit convention with hbar=1.

Inputs, Results, and Numerical Accuracy

• The Hamiltonian can be a supported Observable, an operator-term description, or a dense Hermitian matrix. It is built once, so an arbitrary callable or time-dependent generator is unsupported.
• One request evolves one state on a finite, strictly increasing time grid. Supported initial-state forms include a bitstring, a valid statevector, or a valid density matrix. This interface does not batch independent experiments.
• fql.plan validates and seals the request; fql.run(plan) executes it without semantic overrides. Build a new plan to change the generator, time grid, or initial state.
• The returned EvolutionResult has populations of shape (time_points, 2**n) and named expectation trajectories of shape (time_points,). There is no leading circuit-batch dimension here. Density matrices, when returned, have shape (time_points, 2**n, 2**n).
• Populations are always available as populations / probabilities; explicit output requests currently support expectations. This facade does not accept circuit options such as shots, seed, backend selection, or gradient requests. Copy the restricted ExecutionOptions used in the example rather than a circuit-execution configuration.
• Execution supports CPU complex64 and complex128. It uses fixed-step fourth-order Runge–Kutta with two internal steps per output interval. The example's interval is 0.01, so its integration step is 0.005. There is no adaptive error tolerance.

Refine the time grid and compare the observables to assess convergence. Trace preservation alone does not establish accurate integration, and bounded diagonal populations do not prove that every density-matrix eigenvalue is nonnegative. The script checks eigenvalues for the analytical decay case; the driven run additionally checks returned trace/population diagnostics. Fixed-step RK4 does not provide an unconditional positivity guarantee for arbitrary coarse grids.

This is dense integration: density storage scales as 4**n, and the implementation assembles the time trajectory internally. Requesting return_density_matrices=False omits those matrices from the returned result but does not remove all internal trajectory storage. A finer grid increases both accuracy work and memory use.

8. Limits and Common Problems

Problem or assumption	What to check
A noise rule appears to have no effect	Match the gate name and qubits, then inspect lowered instructions. Out-of-range wires and a model whose gate rules match nothing are rejected; some unused rules in a broader model can remain legal.
A phase-noise result has unchanged probabilities	Inspect coherence or an X/Y expectation; computational-basis populations alone can miss dephasing.
Two libraries disagree at the same noise parameter	Compare channel equations, parameter units, wire order, and readout conventions.
A profile raises noise model identity mismatch	In 0.3.1, construct timing/calibration values as floats consistently and check serialized reconstruction.
Density simulation runs out of memory	One complex64 density matrix needs about 8*4**n bytes before batches and temporaries. Use a reduced exact reference and evaluate trajectory accuracy and representation costs.
A trajectory estimate changes between runs	Record seed, versions, trajectory IDs/count, precision, representation, and truncation settings; distinguish histories from measurement shots.
MPS reaches its bond cap	Inspect truncation/discarded-weight diagnostics and rerun with a larger bond or tighter cutoff. Standard error only describes sampling variation.
Lindblad evolution is inaccurate or nonphysical	Refine the grid, verify units, Hamiltonian Hermiticity, nonnegative rates, and the initial state. Check trace, populations, and eigenvalues where feasible.
A Lindblad GPU, shot, or gradient request is rejected	This facade is CPU-only deterministic dense evolution; those circuit/runtime capabilities do not transfer to it.
Backward appears possible through part of a noisy calculation	General noisy gradients and trainable channel parameters remain outside the documented contract. An accepted flag or a tensor with grad_fn is insufficient evidence of supported differentiation.

The documented boundaries also exclude non-Markovian environment dynamics, leakage, pulse-overlap/crosstalk modeling, and general calibration-provider adapters. The Lindblad solver additionally excludes stochastic trajectories, sparse solvers, GPU/distributed execution, hardware submission, and arbitrary time-dependent generators. The circuit trajectory APIs are separate workflows and do not widen those Lindblad limits.

A synthetic noise simulation is a prediction under its model. It does not establish measured hardware fidelity, error correction, fault tolerance, error mitigation, or general large-system scalability. Match the release capability evidence to the exact workload before making such claims.

Sources and Related Guides

Technical sources below are pinned to v0.3.1:

• Noisy simulation guide: execution, trajectories, readout, and boundaries
• Noise models and calibration · Channel factory definitions
• Noise lowering and timing placement
• Continuous-time Lindblad facade · Integrator and result contract
• Release noisy-circuit example · Release driven Lindblad example
• Capability catalog
• Choosing a Simulation Method · Circuit Compilation and Interoperability · Troubleshooting and FAQ