Quantum Thermodynamics Toy Models

Work measured in a driven quantum system is a distribution of energy differences, not generally the outcome of one fixed “work operator.” This example lists every branch of a two-energy-measurement experiment, verifies the Jarzynski and Crooks relations by exact sums, and compares those sums with simulated records.

Helpful background. The two-point measurement scheme defines the protocol, and energy, heat, and work fixes the thermodynamic interpretation. The calculation uses Python 3 and NumPy.

Two energy measurements around an isolated drive

Use initial and final Hamiltonians $H_0=\epsilon_0|e\rangle\langle e|$ and $H_f=\epsilon_f|e\rangle\langle e|$, with ground energy zero. Start in the Gibbs state of $H_0$ at inverse temperature $\beta>0$, measure $H_0$, apply a unitary $U$, then measure $H_f$. The system is isolated during the drive; its measured energy change is work done on the system:

$$W_{mn}=E_m^f-E_n^0, \qquad p_F(m,n)=\frac{e^{-\beta E_n^0}}{Z_0}|U_{mn}|^2.$$

The equilibrium free-energy difference is $\Delta F=-\beta^{-1}\log(Z_f/Z_0)$. The actual final state need not be thermal: $Z_f$ describes the equilibrium comparison state, not a replacement for the driven state.

Choose energy units with $\epsilon_0=1$, $\epsilon_f=2$, $\beta=\log2$, and $U=e^{-i\theta X/2}$ with $\theta=\pi/3$. Thus the transition probability is $r=\sin^2(\theta/2)=1/4$, the initial populations are $(2/3,1/3)$, and the four branches are

Initial → final	Probability	Work
$g\to g$	$1/2$	$0$
$g\to e$	$1/6$	$2$
$e\to g$	$1/12$	$-1$
$e\to e$	$1/4$	$1$

Their mean work is $1/2$. The negative-work branch is allowed even though the mean satisfies $\langle W\rangle\ge\Delta F$.

Exact sums and the reverse experiment

For this abstract two-level example, take time reversal to be complex conjugation in the energy basis. The reversed drive has $U_R=U^{\mathsf T}$, so its transition probabilities satisfy $P_R(n|m)=P_F(m|n)$. It begins in the Gibbs state of $H_f$ and measures $H_0$ at the end. In a physical model with time-reversal-odd parameters, those parameters must also be reversed; merely replaying an arbitrary pulse is not sufficient.

For a matched branch with nonzero probability,

$$\frac{p_F(m,n)}{p_R(n,m)} =\frac{Z_f}{Z_0}e^{\beta W_{mn}} =e^{\beta(W_{mn}-\Delta F)}.$$

The code stores rows as final indices and columns as initial indices. reverse.T aligns each reverse branch with its forward partner before taking ratios.

import numpy as np

beta = np.log(2.0)
E0 = np.array([0.0, 1.0])
Ef = np.array([0.0, 2.0])
X = np.array([[0, 1], [1, 0]])
theta = np.pi/3
U = np.cos(theta/2)*np.eye(2)-1j*np.sin(theta/2)*X
assert np.allclose(U.conj().T@U, np.eye(2))
P = np.abs(U)**2
assert np.allclose(P.sum(axis=0), 1)
assert np.allclose(P.sum(axis=1), 1)
Z0, Zf = np.exp(-beta*E0).sum(), np.exp(-beta*Ef).sum()
p0, pf = np.exp(-beta*E0)/Z0, np.exp(-beta*Ef)/Zf
forward = P*p0[None, :]
UR = U.T
reverse = np.abs(UR)**2*pf[None, :]
W = Ef[:, None]-E0[None, :]
dF = -np.log(Zf/Z0)/beta
exponential = np.exp(-beta*W)
jarzynski = np.sum(forward*exponential)
mean_work = np.sum(forward*W)
assert np.isclose(forward.sum(), 1)
assert np.isclose(reverse.sum(), 1)
assert np.isclose(mean_work, 0.5)
assert np.isclose(jarzynski, 5/6)
assert np.isclose(jarzynski, Zf/Z0)
mask = (forward > 0) & (reverse.T > 0)
assert np.allclose(np.log(forward[mask]/reverse.T[mask]),
                   beta*(W[mask]-dF))
assert mean_work >= dF
print('mean work:', mean_work, 'free-energy difference:', dF)
print('exponential average:', jarzynski)

# A controlled wrong-sign calculation should fail for this protocol.
wrong_sign = np.sum(forward*np.exp(beta*W))
assert not np.isclose(wrong_sign, Zf/Z0)

rng = np.random.default_rng(20261013)
N = 100000
branches = rng.choice(4, size=N, p=forward.ravel())
samples = exponential.ravel()[branches]
exact_variance = np.sum(forward*(exponential-jarzynski)**2)
se_exact = np.sqrt(exact_variance/N)
se_sample = samples.std(ddof=1)/np.sqrt(N)
assert abs(samples.mean()-jarzynski) < 6*se_exact
print('sample average:', samples.mean(), 'sample SE:', se_sample)

The exact result is $Z_f/Z_0=5/6$, and $\Delta F=\log(6/5)/\log2\approx0.263034$. The wrong-sign exponential average is $41/24\approx1.70833$, which clearly fails this comparison. Keeping such a known failure helps detect an accidentally reversed work convention.

Why the equality holds

Insert the Gibbs weight into the exponential average:

$$\begin{aligned} \langle e^{-\beta W}\rangle &=\frac1{Z_0}\sum_{m,n}|U_{mn}|^2e^{-\beta E_m^f}\\ &=\frac1{Z_0}\sum_m e^{-\beta E_m^f} =\frac{Z_f}{Z_0}. \end{aligned}$$

The second step uses the row sums $\sum_n|U_{mn}|^2=1$, not just the column sums that normalize conditional probabilities. For a two-level transition matrix, nonnegative row and column sums equal to one imply a form $\bigl(\begin{smallmatrix}1-r&r\\r&1-r\end{smallmatrix}\bigr)$ that can come from a unitary. In higher dimensions, being doubly stochastic is necessary but not sufficient to have entries $|U_{mn}|^2$; constructing $U$ explicitly avoids that ambiguity.

If distinct branches share a work value, add their probabilities before plotting the distribution or comparing $P_F(W)/P_R(-W)$. The branch relation then implies the distribution relation because all terms at that work value have the same ratio. Zero-probability branches are excluded from logarithmic ratios.

Sampling and physical limitations

The Monte Carlo estimate is an average of $A=e^{-\beta W}$, whose standard error is $\sqrt{\operatorname{Var}(A)/N}$. In this small model the exact variance is available, so it can check the sample-based estimate. At low temperature, rare negative-work events can dominate the exponential average. If the sample misses them, even its estimated error bar can be misleading; inspect expected branch counts, not just the apparent stability of the sample mean.

The initial energy measurement leaves a Gibbs state unchanged. For a coherent initial state it removes energy-basis coherence, and the standard thermal Jarzynski argument above no longer applies without its preparation assumptions. Likewise, adding dissipative evolution makes a system-only energy difference include heat; it cannot simply be renamed the same isolated-system work. See fluctuation theorems for the additional bath and protocol assumptions.

Changing the common energy zero at both endpoints leaves work unchanged. Changing the two endpoint zeros by different amounts shifts both $W$ and $\Delta F$ by the same difference, preserving the fluctuation relations when done consistently.

Exercises

Check the branch sum. Verify $\langle W\rangle=1/2$ and $\langle e^{-\beta W}\rangle=5/6$ from the table.

Solution

The mean is $0+(1/6)2-(1/12)+(1/4)=1/2$. Since $e^{-\beta W}=2^{-W}$, the exponential average is $1/2+1/24+1/6+1/8=5/6$.

A mean-work bound. Derive the inequality $\langle W\rangle\ge\Delta F$.

Solution

Convexity gives $e^{-\beta\langle W\rangle}\le\langle e^{-\beta W}\rangle=e^{-\beta\Delta F}$. Taking logarithms and dividing by $-\beta<0$ reverses the inequality. This uses $\beta>0$.

Estimate the statistical error. Find the exact variance of $A=e^{-\beta W}$ for the table.

Solution

$\langle A^2\rangle=1/2+1/96+1/3+1/16=29/32$. Subtracting $(5/6)^2$ gives $\operatorname{Var}(A)=61/288$. At $N=100000$, the standard error is about $0.001455$.

A doubly stochastic counterexample. Explain why the $3\times3$ matrix with rows $(1/2,1/2,0)$, $(1/2,0,1/2)$, $(0,1/2,1/2)$ cannot equal $|U_{mn}|^2$ for a unitary $U$.

Solution

Its first and second columns overlap in only the first row. Any putative unitary columns therefore have an inner product consisting of one nonzero term of magnitude $1/2$. Phases cannot cancel it, so the columns cannot be orthogonal, despite all row and column sums being one.

References

• M. Campisi, P. Hänggi, and P. Talkner, “Colloquium: Quantum fluctuation relations: Foundations and applications,” Reviews of Modern Physics 83, 771–791 (2011). doi:10.1103/RevModPhys.83.771.
• G. E. Crooks, “Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences,” Physical Review E 60, 2721–2726 (1999). doi:10.1103/PhysRevE.60.2721.
• C. Jarzynski, “Nonequilibrium equality for free energy differences,” Physical Review Letters 78, 2690–2693 (1997). doi:10.1103/PhysRevLett.78.2690.
• P. Talkner, E. Lutz, and P. Hänggi, “Fluctuation theorems: Work is not an observable,” Physical Review E 75, 050102(R) (2007). doi:10.1103/PhysRevE.75.050102.
• H. Tasaki, “Jarzynski relations for quantum systems and some applications,” arXiv:cond-mat/0009244 (2000). Author manuscript.