Separate the initial condition rho(0) from the driving term

[Jirnyak]

rho(0) was hardcoded in run_simulation, and the lead bias did double duty: it set both the initial fill and the target the dissipator drives towards. Those are different objects in the Cauchy problem — the initial data $y_0$ and a coefficient $b$ of the equation $\dot y = Ay + b$ — so varying one always dragged the other along, and "same equation, different initial condition" was not expressible.

Changes

• build_initial_rho replaces the five hardcoded lines, with three recipes:
  • from_bias — each lead at its own bias, molecule empty. Previous behaviour, bit-identical, still the default.
  • block — every block including the extended molecule filled to one common level. Removes the charging transient.
  • connected — $f(H,\mu,T)$ over the eigenstates of the whole Hamiltonian, rotated into the state representation. Lead–EM coherences appear on their own and $0 \preceq \rho \preceq 1$ holds by construction.
• Reservoir T moves from Simulation Parameters up next to the bias fields, where it belongs — it is part of the driving term, not of the numerics.
• Sampling order — the loop recorded a point after stepping, shifting the whole plot by one dt and omitting the true $t=0$ point. Now it records first.
• Fermi occupation factored into fermi() instead of being inlined.

Verification

Built clean on macOS, app launches. Numerics checked against a 106-site system:

recipe	trace	EM pop.	eig range	$\lVert\dot\rho\rVert$ at $\gamma=0$
from_bias	40	0	[0, 1]	0.4297
block, $\mu=0$	53	13	[0, 1]	0.4297
connected, $\mu=0$	53	13	[0, 1]	0
connected, $T=300$K	53	13	[2e-7, 1−2e-7]	0

connected is an exact stationary state of the closed system, as it must be since $f(H)$ commutes with $H$. The legacy path reproduces the old hardcoded matrix exactly (zero difference). The from_bias row shows the artifact this addresses: 13 of 26 extended-molecule levels start empty and flood on connection, a charging transient unrelated to the applied bias.

Equal derivative norms in rows 1–2 are a particle-hole coincidence of the half-filled symmetric band, not an identity — off half-filling ($\mu = 0.08$) they separate to 0.4297 vs 0.4113.

Note

The new recipes need $\gamma > 0$. The bias enters the equation only through $b = \Gamma\rho^0$, so at $\gamma = 0$ it has no effect at all once it is no longer also setting $\rho(0)$.

🤖 Generated with Claude Code

[github-actions]

🚀 Thank you for opening your first Pull Request! Our maintainers will review your changes shortly.