Models and spectra

This guide constructs each Hamiltonian model in pmm_tools and computes many-body energy spectra. The Python cells execute when the documentation builds. The parameter choices illustrate the APIs; they are not fits to a device.

Model	Constructor	Fermionic modes
Isolated quadratic Kitaev chain	`kitaev_chain.chain_operators` or `symbolic_chain`	$N$ spinless sites
Interacting Kitaev chain	`interacting_kitaev_chain`	$N$ spinless sites with nearest-neighbor interactions
Spinful dot-ABS chain	`dot_abs_chain`	Two spin modes per site
Parity qubit with a middle dot	`parity_qubit_microscopic.device_operators`	Two $N$ -site chains and one spinless dot
Parity qubit with a direct link	`parity_qubit_tunnel.device_matrix_dict`	Two $N$ -site chains joined by normal hopping

Units and basis

All numerical energies below are in μeV, times in ns, and superconducting phases in radians. Hamiltonian constructors do not attach units: every chemical potential, hopping, pairing amplitude, and interaction in one calculation must use the same energy unit.

For a level separation $\Delta E$ , the angular frequency and ordinary frequency are

\omega = \frac{\Delta E}{\hbar},\qquad f = \frac{\Delta E}{2\pi\hbar}.

HBAR_UEV_NS is $\hbar$ in μeV·ns. Thus $\omega$ is in rad/ns and $f$ in ns $^{-1}$ , numerically GHz. Unitary evolution uses $\exp(-iHt/\hbar)$ . A decay rate expressed in ns $^{-1}$ converts to s $^{-1}$ by multiplying by $10^9$ ; a decay rate has no extra $2\pi$ factor.

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
from second_quantization import hilbert_space

from pmm_tools.models import dot_abs_chain, interacting_kitaev_chain
from pmm_tools.models import kitaev_chain as chain
from pmm_tools.models import parity_qubit_microscopic as dot_qubit
from pmm_tools.models import parity_qubit_tunnel as tunnel_qubit
from pmm_tools.parity import split_parity_blocks
from pmm_tools.units import HBAR_UEV_NS

np.set_printoptions(precision=5, suppress=True)
print(f"hbar = {HBAR_UEV_NS:.10f} μeV ns")
print(f"1 μeV / hbar = {1 / HBAR_UEV_NS:.6f} rad/ns")
print(f"1 μeV / (2 pi hbar) = {1 / (2 * np.pi * HBAR_UEV_NS):.6f} GHz")

hbar = 0.6582119570 μeV ns
1 μeV / hbar = 1.519267 rad/ns
1 μeV / (2 pi hbar) = 0.241799 GHz

For $M$ fermionic modes, the full occupation basis has dimension $2^M$ . The first mode is the most significant occupation bit: two modes are ordered $|00\rangle, |01\rangle, |10\rangle, |11\rangle$ . Physical parity is $P=(-1)^{\sum_j n_j}$ : even is $+1$ , odd is $-1$ . second_quantization returns the labels 0 and 1 instead, so we explicitly convert its parity matrix below. Pairing changes particle number by two and preserves parity. These examples diagonalize complete parity blocks, without projecting onto a two-level qubit.

Single Kitaev chain

The quadratic spinless chain uses

H = \sum_j \mu_j c_j^\dagger c_j + \sum_j\left(t_j c_j^\dagger c_{j+1} + \Delta_j e^{i\phi}c_j^\dagger c_{j+1}^\dagger + \mathrm{h.c.}\right).

The onsite sign is $+\mu n$ . chain_operators returns matrices multiplying real control values, with one-indexed parameter names and a fixed pairing phase. Here three sites give an $8\times8$ Hamiltonian and two $4\times4$ parity blocks.

chain_ops, chain_parity = chain.chain_operators(3, phi=0.0)
chain_parameters = {
    "mu_1": 6.0, "mu_2": 6.0, "mu_3": 6.0,
    "t_1": 20.0, "t_2": 20.0,
    "Delta_1": 20.0, "Delta_2": 20.0,
}
H_chain = sum(chain_parameters[name] * op for name, op in chain_ops.items())
chain_blocks = split_parity_blocks(H_chain, chain_parity)
chain_energies = np.linalg.eigvalsh(chain_blocks)
np.testing.assert_allclose(H_chain, H_chain.conj().T, atol=1e-12)
print("Even energies (μeV):", chain_energies[0])
print("Odd energies (μeV): ", chain_energies[1])

Even energies (μeV): [-31.6293    6.13197  12.       49.49734]
Odd energies (μeV):  [-31.49734   6.       11.86803  49.6293 ]

For symbolic calculations, symbolic_chain returns the expression, ordered modes, and a symbol dictionary. This interface uses zero-indexed names and allows an independent phase on each bond. For example, a two-site chain is:

H_symbolic, modes, symbols = chain.symbolic_chain(2)
print("Symbols:", tuple(symbols))
display(H_symbolic)

Symbols: ('mu_0', 'mu_1', 't_0', 'Delta_0', 'phi_0')

\displaystyle \Delta_{0} e^{i \phi_{0}} {{c_0}^\dagger} {{c_1}^\dagger} + \Delta_{0} e^{- i \phi_{0}} {c_1} {c_0} + \mu_{0} {{c_0}^\dagger} {c_0} + \mu_{1} {{c_1}^\dagger} {c_1} + t_{0} {{c_0}^\dagger} {c_1} + t_{0} {{c_1}^\dagger} {c_0}

For a real, noninteracting chain, the singular values of majorana_block are positive quasiparticle excitation energies. They are not the full many-body spectrum. parity_gap returns the signed ground-state difference $E_{\mathrm{even},0}-E_{\mathrm{odd},0}$ ; its sign identifies the lower sector. The reduced representation should agree with the Fock-space calculation:

B = chain.majorana_block([6.0, 6.0, 6.0], 20.0, 20.0)
quasiparticle_energies = np.sort(np.linalg.svd(B, compute_uv=False))
fock_gap = chain_energies[0, 0] - chain_energies[1, 0]
np.testing.assert_allclose(chain.parity_gap(B), fock_gap, atol=1e-12)
print("Quasiparticle energies (μeV):", quasiparticle_energies)
print(f"Signed parity gap: {fock_gap:.6f} μeV")

Quasiparticle energies (μeV): [ 0.13197 37.6293  43.49734]
Signed parity gap: -0.131967 μeV

We next vary the common onsite energy. Every plotted level uses the same energy reference at each point: the lowest energy across both parity sectors. At $\mu=0$ and $t=\Delta$ , the open chain has a zero mode and matching even/odd spectra.

mu_values = np.linspace(-40.0, 40.0, 161)
chain_spectra = []
for mu in mu_values:
    parameters = chain_parameters | {f"mu_{j}": mu for j in range(1, 4)}
    H = sum(parameters[name] * op for name, op in chain_ops.items())
    chain_spectra.append(np.linalg.eigvalsh(split_parity_blocks(H, chain_parity)))
chain_spectra = np.asarray(chain_spectra)
zero_index = np.flatnonzero(mu_values == 0.0)[0]
np.testing.assert_allclose(
    chain_spectra[zero_index, 0], chain_spectra[zero_index, 1], atol=1e-12
)
chain_excitations = chain_spectra - chain_spectra.min(axis=(1, 2))[:, None, None]

fig, ax = plt.subplots(figsize=(6.4, 3.8), layout="constrained")
for sector, color, style, label in [(0, "C0", "-", "Even"), (1, "C1", "--", "Odd")]:
    lines = ax.plot(mu_values, chain_excitations[:, sector], color=color, ls=style)
    lines[0].set_label(label)
ax.set(xlabel="Uniform onsite energy μ (μeV)", ylabel="E − E₀ (μeV)",
       title="Three-site Kitaev chain: t = Δ = 20 μeV")
ax.legend()
plt.show()

Figure generated by the preceding code cell

Interacting Kitaev chain

interacting_kitaev_chain adds $\sum_j U_j n_j n_{j+1}$ . Its parameter mapping accepts scalars, site/bond sequences, "constant" for a uniform symbol, or "space-dependent" for indexed symbols. It returns (Hamiltonian, fermions, parameters), with parameter values normalized to tuples. An explicit mapping supplies all four keys: mu, t, Delta, and U.

Here we compare zero interaction with $U=8$ μeV. We construct the symbolic Hamiltonian, convert its terms to Fock matrices, and diagonalize the parity blocks. The $U=0$ result must reproduce the previous quadratic chain.

interaction_spectra = []
for U in (0.0, 8.0):
    expression, modes, parameters = interacting_kitaev_chain(
        3, {"mu": [6.0, 6.0, 6.0], "t": 20.0, "Delta": 20.0, "U": U}
    )
    terms = hilbert_space.to_matrix(expression, modes, sparse=False)
    H = sum(complex(coefficient) * matrix for coefficient, matrix in terms.items())
    parity = np.eye(H.shape[0]) - 2 * hilbert_space.parity_operator(modes, sparse=False)
    energies = np.linalg.eigvalsh(split_parity_blocks(H, parity))
    interaction_spectra.append(energies)
    print(f"U = {U:g} μeV; even/odd energies (μeV):\n{energies}")
np.testing.assert_allclose(interaction_spectra[0], chain_energies, atol=1e-12)

U = 0 μeV; even/odd energies (μeV):
[[-31.6293    6.13197  12.       49.49734]
 [-31.49734   6.       11.86803  49.6293 ]]
U = 8 μeV; even/odd energies (μeV):
[[-28.22772   6.30728  20.       53.92044]
 [-29.30061   6.       18.48806  56.81255]]

The Majorana singular-value reduction does not apply when $U\ne0$ .

Spinful dot-ABS chain

dot_abs_chain orders modes as (site0 up, site0 down, site1 up, site1 down, ...). On each site, the spin energies are $\mu+E_z$ and $\mu-E_z$ , so the Zeeman splitting is $2E_z$ . Delta creates an onsite spin pair, U couples the opposite-spin occupations, and t_n and t_so describe spin-preserving and spin-flip hopping. Delta may be complex; do not apply a phase twice.

Two spinful sites give 16 Fock states, with eight states per parity sector.

expression, spin_modes, spin_parameters = dot_abs_chain(2, {
    "mu": [0.0, -15.0], "E_z": [30.0, 0.0],
    "Delta": [0.0, 20.0 * np.exp(0.3j)], "U": [40.0, 0.0],
    "t_n": 5.0, "t_so": 2.0,
})
terms = hilbert_space.to_matrix(expression, spin_modes, sparse=False)
H_spinful = sum(complex(coefficient) * matrix for coefficient, matrix in terms.items())
spin_parity = np.eye(16) - 2 * hilbert_space.parity_operator(spin_modes, sparse=False)
np.testing.assert_allclose(H_spinful, H_spinful.conj().T, atol=1e-12)
spin_spectra = np.linalg.eigvalsh(split_parity_blocks(H_spinful, spin_parity))
print("Modes:", spin_modes)
print("Even/odd energies (μeV):\n", spin_spectra)

Modes: [c_{0\uparrow}, c_{0\downarrow}, c_{1\uparrow}, c_{1\downarrow}]
Even/odd energies (μeV):
 [[-48.39183 -45.      -37.50233  -0.25546   9.17824  15.       16.06563
   50.90575]
 [-70.34826 -22.99303 -16.38821 -12.03024  -9.69356  25.29293  25.37154
   40.78884]]

Parity qubit with a middle dot

The microscopic device consists of two Kitaev chains connected through one spinless dot. The modes are (L1, ..., LN, M, R1, ..., RN). t_LM and t_MR are normal hoppings to the middle dot; mu_M is its onsite energy. The left pairing phase is zero and phi_R is the right phase. For $N=2$ , this gives 32 states, or 16 in each total-parity sector.

The logical even-sector states are usually prepared from the isolated-chain products ee and oo, with the middle dot empty. At finite coupling, the two lowest eigenstates need not remain in that logical subspace. The spectrum below retains the dot and all excited chain states.

dot_ops, dot_parity = dot_qubit.device_operators(2, phi_R=0.4)
dot_parameters = dict.fromkeys(dot_ops, 0.0)
dot_parameters.update(
    t_L_1=20.0, Delta_L_1=20.0, t_R_1=20.0, Delta_R_1=20.0,
    mu_L_1=2.0, mu_R_2=-3.0, mu_M=60.0, t_LM=5.0, t_MR=5.0,
)
dot_values = np.linspace(-40.0, 100.0, 141)
dot_spectra = []
for mu_M in dot_values:
    parameters = dot_parameters | {"mu_M": mu_M}
    H = sum(parameters[name] * op for name, op in dot_ops.items())
    dot_spectra.append(np.linalg.eigvalsh(split_parity_blocks(H, dot_parity)))
dot_spectra = np.asarray(dot_spectra)
dot_excitations = dot_spectra - dot_spectra.min(axis=(1, 2))[:, None, None]
print("Full dimension:", H.shape[0], "; states per parity:", dot_spectra.shape[-1])

Full dimension: 32 ; states per parity: 16

fig, ax = plt.subplots(figsize=(6.4, 3.8), layout="constrained")
for sector, color, style, label in [(0, "C0", "-", "Even"), (1, "C1", "--", "Odd")]:
    lines = ax.plot(dot_values, dot_excitations[:, sector, :4], color=color, ls=style)
    lines[0].set_label(label)
ax.set(xlabel="Middle-dot onsite energy μM (μeV)", ylabel="E − E₀ (μeV)",
       title="Chain–dot–chain: lowest four levels per parity")
ax.legend()
plt.show()

Parity qubit with a direct link

parity_qubit_tunnel removes the middle dot and connects the chains through $t_c(c_{LN}^\dagger c_{R1}+\mathrm{h.c.})$ , without pairing on that link. Both chain phases are zero. The modes are (L1, ..., LN, R1, ..., RN). For $N=2$ , each total-parity block has eight states.

device_matrix_dict returns sparse coefficient matrices. We compile them once with second_quantization and reuse them across a coupling sweep. clean_parameters(..., gap=40) sets each chain to $t=\Delta=20$ μeV and zero onsite energies: gap denotes the clean bulk excitation energy, not the zero-mode splitting or the coupled qubit splitting.

tunnel_evaluate = {
    sector: hilbert_space.make_dict_callable(tunnel_qubit.device_matrix_dict(2, sector))
    for sector in ("even", "odd")
}
couplings = np.linspace(0.0, 20.0, 101)
tunnel_spectra = []
for coupling in couplings:
    parameters = tunnel_qubit.clean_parameters(2, coupling=coupling, gap=40.0)
    blocks = [tunnel_evaluate[sector](**parameters).toarray() for sector in ("even", "odd")]
    tunnel_spectra.append(np.linalg.eigvalsh(blocks))
tunnel_spectra = np.asarray(tunnel_spectra)
tunnel_excitations = tunnel_spectra - tunnel_spectra.min(axis=(1, 2))[:, None, None]
np.testing.assert_allclose(tunnel_excitations[0, :, :2], 0.0, atol=1e-12)
print("Even/odd energies at t_c = 20 μeV:\n", tunnel_spectra[-1])

Even/odd energies at t_c = 20 μeV:
 [[-51.23106 -31.23106 -20.       -0.        0.       20.       31.23106
   51.23106]
 [-51.23106 -31.23106 -20.       -0.        0.       20.       31.23106
   51.23106]]

/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")

fig, ax = plt.subplots(figsize=(6.4, 3.8), layout="constrained")
for sector, color, style, label in [(0, "C0", "-", "Even"), (1, "C1", "--", "Odd")]:
    lines = ax.plot(couplings, tunnel_excitations[:, sector, :4], color=color, ls=style)
    lines[0].set_label(label)
ax.set(xlabel="Direct hopping tc (μeV)", ylabel="E − E₀ (μeV)",
       title="Directly coupled chains: lowest four levels per parity")
ax.legend()
plt.show()

At zero coupling, the products ee and oo are degenerate in the even sector; eo and oe are degenerate in the odd sector. Exact sweet-spot chains also retain matching even/odd spectra at finite coupling because the outer zero modes remain uncoupled. The plotted lines therefore overlap between sectors. The finite Fock calculation remains valid for strong links; a two-level interpretation requires a separate separation-of-scales and leakage check.

Noise and coherence

models.quasistatic supplies Gaussian Sobol samples and model-specific noisy Hamiltonian coefficients. Each realization is static during a shot; averaging phases across shots describes inhomogeneous dephasing. coherence.first_envelope_crossing locates the first $1/e$ envelope crossing of the supplied frequency ensemble. This is distinct from energy relaxation: pmm_tools does not yet provide thermal Lindblad jumps or dissipative time evolution.