Time evolution

pmm_tools.evolution.evolve integrates the pure-state Schrödinger equation for

$$H(t)=H_0+\sum_j u_j(t)O_j,\qquad \frac{d|\psi\rangle}{dt}=-\frac{i}{\hbar}H(t)|\psi\rangle.$$

The controls $u_j$ interpolate linearly between supplied knots. Each knot is an integration boundary, even when output times omit it. SciPy's DOP853 solver adapts its steps within each interval. No state projection or renormalization is applied; the integration tolerances determine the numerical error.

Use μeV for $H_0$ and the products $u_jO_j$, and ns for knots and output times. For occupation operators, $O_j$ is dimensionless and $u_j$ is a chemical potential in μeV. All operators and the initial state must use the same basis. This routine evolves one state vector; it does not average noise realizations or evolve a dissipative density matrix.

Pulse a parity qubit's middle dot

We use the chain–dot–chain model with two sites per chain. Its full Fock ordering is (L1, L2, M, R1, R2), with the first occupation as the most significant bit. All terms conserve total parity, so we retain the complete 16-state even sector. There is no two-level approximation. The initial state is the product of the two isolated even ground states and an empty middle dot.

import numpy as np
import matplotlib.pyplot as plt

from pmm_tools.evolution import evolve
from pmm_tools.models import parity_qubit_microscopic as qubit

operators, parity = qubit.device_operators(2, phi_R=0.4)
parameters = dict.fromkeys(operators, 0.0)
parameters.update(t_L_1=20.0, Delta_L_1=20.0,
                  t_R_1=20.0, Delta_R_1=20.0, t_LM=2.0, t_MR=2.0)
indices = qubit.basis_indices(2, "even")
full_static = sum(parameters[name] * matrix for name, matrix in operators.items())
H_static = full_static[np.ix_(indices, indices)]
number_middle = operators["mu_M"][np.ix_(indices, indices)]
psi0 = qubit.clean_product_state(2, parameters | {"phi_R": 0.4}, label="ee")

/builds/qt/pmm_tools/.pixi/envs/docs/lib/python3.13/site-packages/second_quantization/hilbert_space.py:74: FutureWarning: Input has data type uint8, but the output has been cast to float64.  In the future, the output data type will match the input. To avoid this warning, set the `dtype` parameter to `None` to have the output dtype match the input, or set it to the desired output data type.
  return scipy.sparse.diags(parity_diag, format="csr")

The dot energy ramps from 60 μeV to zero in 2 ns, stays there for 4 ns, and returns in 2 ns. The links remain on throughout. mu_M is zero in H_static, so the control values specify its absolute energy without double-counting a baseline. values has shape (number of knots, number of controls).

knots = np.array([0.0, 2.0, 6.0, 8.0])
values = np.array([[60.0], [0.0], [0.0], [60.0]])
times = np.linspace(0.0, 8.0, 401)
states = evolve(H_static, [number_middle], knots, values, psi0, times,
                rtol=1e-10, atol=1e-12)
print("State-array shape:", states.shape)

State-array shape: (401, 16)

We compute the dot occupation and the return probability to the prepared state. The latter includes coherent transfer and excitation within the retained parity sector; a decrease does not by itself imply relaxation or logical leakage.

norms = np.sum(abs(states)**2, axis=1)
occupation = np.einsum("ti,ij,tj->t", states.conj(), number_middle, states).real
return_probability = abs(states @ psi0.conj())**2
np.testing.assert_allclose(norms, 1.0, atol=1e-8, rtol=0)
print(f"Maximum norm error: {np.max(abs(norms - 1)):.2e}")
print(f"Final middle-dot occupation: {occupation[-1]:.6f}")
print(f"Final return probability: {return_probability[-1]:.6f}")

Maximum norm error: 7.83e-10
Final middle-dot occupation: 0.260191
Final return probability: 0.354384

fig, axes = plt.subplots(2, 1, figsize=(6.4, 5.0), sharex=True, layout="constrained")
axes[0].plot(knots, values[:, 0], "o-")
axes[0].set(ylabel="Middle-dot energy (μeV)", title="Coherent chain–dot–chain evolution")
axes[1].plot(times, occupation, label="Dot occupation")
axes[1].plot(times, return_probability, label="Return probability")
axes[1].set(xlabel="Time (ns)", ylabel="Occupation / probability", ylim=(-0.03, 1.03))
axes[1].legend()
plt.show()

Output times must increase strictly and include both endpoints of the control interval. Interior control knots need not be output times. For a static Hamiltonian, pass an empty operator list and values=np.empty((2, 0)) with two endpoint knots. Pulse construction, preparation choices, noise sampling, and readout remain explicit in the calling code.