At fixed $V,N$, differentiate $S(U)=\tfrac32 Nk\ln U+\text{const}$:

Inverting gives $U=\tfrac32 NkT$ — exactly equipartition, $\tfrac12 kT$ for each of the three translational degrees of freedom. Temperature is therefore not assumed but emerges as the reciprocal slope of the entropy. For $N=1$ at $U=3\times10^{-21}\,$J this returns $T=2U/(3k)=144.86$ K with $\tfrac32 kT=U$, matching temperature_from_entropy(lambda U: 1.5NK_B*log(U), U).

The first law $dU=T\,dS-P\,dV$ makes $(S,V)$ the natural variables of $U$. Each Legendre transform adds a conjugate product to swap one variable for the other:

Reading off the constructions, $G=(U+PV)-TS=H-TS$ and $G=(U-TS)+PV=F+PV$. Numerically at $(U,T,S,P,V)=(100,300,0.2,10^5,10^{-3})$: $H=200$, $F=40$, $G=140$, with $H-TS=200-60=140=F+PV=40+100$ — so gibbs_free_energy(U,T,S,P,V) equals both enthalpy(U,P,V) − TS and helmholtz_free_energy(U,T,S) + PV.

From $dF=-S\,dT-P\,dV$ the first derivatives are $S=-(\partial F/\partial T)_V$ and $P=-(\partial F/\partial V)_T$. Because mixed second partials of $F$ commute, $\partial^2F/\partial V\,\partial T=\partial^2F/\partial T\,\partial V$, so

For the ideal gas $S=Nk\big(\ln V+\tfrac32\ln T\big)$ gives $(\partial S/\partial V)_T=Nk/V$, and $P=NkT/V$ gives $(\partial P/\partial T)_V=Nk/V$ — identical. The two sides agree to a relative residual $\sim3\times10^{-5}$, so maxwell_relation_residual(S_TV, P_TV, T, V) $\approx 0$.

A reversible cycle returns to its initial state, so its total entropy change is zero and the heats obey $Q_h/T_h=Q_c/T_c$, i.e. $Q_c/Q_h=T_c/T_h$. With energy conservation $W=Q_h-Q_c$, the efficiency is

Run in reverse as a refrigerator, $\text{COP}=Q_c/W=(T_c/T_h)/(1-T_c/T_h)=T_c/(T_h-T_c)$; as a heat pump it delivers $Q_h$, so $Q_h/W=T_h/(T_h-T_c)=1+T_c/(T_h-T_c)=1+\text{COP}_{\text{fridge}}$. Thus carnot_efficiency(300,600) $=1-\tfrac{300}{600}=0.5$ and carnot_cop_heat_pump(270,300) $=\tfrac{300}{30}=10=1+9=$ 1 + carnot_cop_refrigerator(270,300).

Take both reservoirs large enough that their temperatures hold fixed, so each entropy change is $\Delta S=Q_{\text{in}}/T$. The hot reservoir loses $Q$ and the cold one gains it:

For $T_h>T_c$ the numerator is positive, so $\Delta S_{\text{total}}>0$: heat flows hot$\to$cold spontaneously, never the reverse — the second law. As $T_h\to T_c$ the gap $T_h-T_c\to0$ and the transfer becomes reversible, $\Delta S\to0$. Hence total_entropy_change_heat_flow(10, 500, 300) $=10\big(\tfrac1{300}-\tfrac1{500}\big)\approx1.33\times10^{-2}$ J/K $>0$, shrinking toward zero as the reservoir temperatures converge.

Writing $x\equiv\beta\varepsilon=\varepsilon/kT$, the three Boltzmann weights $e^{-\beta E_i}$ are $1,\,e^{-x},\,e^{-2x}$, so

Since $e^{-x}<1$ for every $T>0$, the weights fall with energy and $P_0>P_1>P_2$. As $T\to0$, $x\to\infty$ so $e^{-x},e^{-2x}\to0$, giving $Z\to1$ and $P_0\to1$ — the system freezes into the ground state. As $T\to\infty$, $x\to0$ so every weight $\to1$, $Z\to3$ and each $P_n\to\tfrac13$ — the levels equalize. This is exactly what boltzmann_probability([0,1e-21,2e-21], T) returns: a list summing to $1$ with $P_0>P_1>P_2$.

Differentiate $\ln Z$ for $Z=\sum_i e^{-\beta E_i}$, using $\partial_\beta e^{-\beta E_i}=-E_i e^{-\beta E_i}$:

So the thermodynamic energy is literally the log-derivative of the partition function. The closed form $\sum_i E_iP_i$ computed by internal_energy therefore equals the finite-difference $-d(\ln Z)/d\beta$ for the 4-level system (e.g. $5.69\times10^{-22}$ J at $T=80$ K).

With $U=-\partial_\beta\ln Z$ and $\beta=1/kT$ (so $d\beta/dT=-1/kT^2$), the chain rule gives

The second derivative of $\ln Z$ is the energy variance (its second cumulant), $\partial^2_\beta\ln Z=\langle E^2\rangle-\langle E\rangle^2=\operatorname{Var}(E)$, hence

A response function equals an equilibrium fluctuation — the simplest fluctuation–dissipation relation. The variance form heat_capacity(levels, T) thus reproduces the numerical $dU/dT$ (e.g. $7.79\times10^{-24}$ J/K for the 4-level system at $T=80$ K).

With levels $\{0,\varepsilon\}$ and $x=\beta\varepsilon=\varepsilon/kT$, $Z=1+e^{-x}$ and only the upper level carries energy:

The variance is $\operatorname{Var}(E)=\langle E^2\rangle-\langle E\rangle^2=\varepsilon^2/(e^x+1)-\varepsilon^2/(e^x+1)^2=\varepsilon^2 e^{x}/(e^{x}+1)^2$, so

At low $T$ ($x\to\infty$) this $\sim kx^2e^{-x}\to0$ (the gap is frozen); at high $T$ ($x\to0$) it $\sim kx^2/4\to0$ (both levels equally full). So $C$ vanishes at both ends and bumps between, peaking near $x\approx2.4$, i.e. $kT\approx0.42\,\varepsilon$. Accordingly two_level_heat_capacity(eps, T) is largest near $kT\approx0.42\,\varepsilon$ ($C/k\approx0.44$ at $T\approx30$ K, with $\varepsilon/k=72.4$ K) and small at $T=1$ K and $T=10^8$ K.

Let $y=\beta\hbar\omega=\hbar\omega/kT$, so $\langle E\rangle=\hbar\omega\big(\tfrac12+\tfrac1{e^{y}-1}\big)$. As $T\to0$, $y\to\infty$ and $1/(e^y-1)\to0$, leaving the zero-point energy $\hbar\omega/2$. As $T\to\infty$, $y\to0$ and $e^y-1\approx y$, so $1/(e^y-1)\approx kT/\hbar\omega$ and

classical equipartition. The Einstein heat capacity $C=k\,y^2 e^{y}/(e^{y}-1)^2\to k$ as $y\to0$ (Dulong–Petit) and $\sim k\,y^2 e^{-y}\to0$ for $y\gg1$ (frozen out below $\hbar\omega/k$). Each of three quadratic translational degrees of freedom adds $\tfrac12 kT$, so $U=\tfrac32 NkT$. Numerically ($\hbar\omega/k=76.4$ K): harmonic_oscillator_energy gives $\langle E\rangle/\hbar\omega=0.50$ at $T=10$ K, einstein_heat_capacity$/k\to1.0$ at high $T$ and $\sim6.5\times10^{-3}$ at $T=8$ K, and equipartition_energy(3, N, T) $=\tfrac32 NkT$.

Let $x\equiv(\varepsilon-\mu)/kT$. The three mean occupations are

For $x\gg1$, $e^x\gg1$ so the $\mp1$ is negligible beside $e^x$ and both quantum forms collapse to $1/e^x=e^{-x}$, the classical Maxwell–Boltzmann factor. For finite $x$ the denominators order as $e^x-1<e^x<e^x+1$, so the reciprocals give

bosons bunch (no exclusion), fermions are suppressed (Pauli). Tie to check: at $(\varepsilon-\mu)/kT=10$ the three functions agree to $\sim10^{-4}$, all $\approx e^{-10}=4.540\times10^{-5}$.

Write $\langle n\rangle_{\rm FD}=1/(e^x+1)$ with $x=(\varepsilon-\mu)/kT$. For every real $x$, $e^x>0$, so the denominator obeys $e^x+1>1$ and is finite, giving

The ceiling $\langle n\rangle\le1$ is the Pauli principle — at most one fermion per state. The extremes are $\langle n\rangle\to1$ for $\varepsilon\ll\mu$ ($x\to-\infty$, filled) and $\langle n\rangle\to0$ for $\varepsilon\gg\mu$ ($x\to+\infty$, empty). At the Fermi level $\varepsilon=\mu$, $x=0$ and $e^0=1$, so

Tie to check: fermi_dirac(mu,mu,T) $=0.5$, and the occupation stays in $[0,1]$ for all $\varepsilon$.

At $T=0$ the chemical potential is the Fermi energy, $\mu(0)=\varepsilon_F$, so $\langle n\rangle_{\rm FD}=1/\!\left(e^{(\varepsilon-\varepsilon_F)/kT}+1\right)$. Fix $\varepsilon\neq\varepsilon_F$ and let $T\to0^+$, so $1/kT\to+\infty$. If $\varepsilon<\varepsilon_F$ the exponent $\to-\infty$, $e^{(\varepsilon-\varepsilon_F)/kT}\to0$ and $\langle n\rangle\to1$; if $\varepsilon>\varepsilon_F$ the exponent $\to+\infty$ and $\langle n\rangle\to0$:

The smooth Fermi function sharpens into a step whose width $\sim kT$ shrinks to zero — every state below $\varepsilon_F$ filled, every one above empty (the Fermi sea). Tie to check: fermi_dirac_T0(0.9eF, eF) $=1$ and fermi_dirac_T0(1.1eF, eF) $=0$.

The Planck density is $u(\omega)=\dfrac{\hbar\omega^3}{\pi^2c^3}\dfrac{1}{e^{\hbar\omega/kT}-1}$. Put $x=\hbar\omega/kT$; at low frequency $x\ll1$, so $e^x-1=x+\tfrac{x^2}{2}+\cdots\approx x$ and $1/(e^x-1)\approx1/x=kT/\hbar\omega$. Then

the Rayleigh–Jeans law — note that $\hbar$ has cancelled. It rises as $\omega^2$ without bound, so $\int_0^\infty u\,d\omega\propto\int\omega^2\,d\omega$ diverges (the ultraviolet catastrophe). Planck's denominator $e^{\hbar\omega/kT}$ rescues it: modes with $\hbar\omega\gg kT$ are exponentially frozen out. Tie to check: planck_energy_density $\approx$ rayleigh_jeans for $\omega=10^{10}$–$10^{12}$ rad/s at $T=5778$ K (ratio $0.99993$ at $10^{11}$).

Write $u\propto x^3/(e^x-1)$ with $x=\hbar\omega/kT$. The peak is where $du/d\omega=0$:

whose nonzero root is $x=2.821439$ — Wien's displacement law, $\omega_{\max}\propto T$. For the total density substitute $x=\hbar\omega/kT$ ($d\omega=(kT/\hbar)\,dx$) and use $\int_0^\infty x^3/(e^x-1)\,dx=\pi^4/15$:

The radiated flux is $\sigma T^4$ with $\sigma=ac/4=\pi^2k^4/60\hbar^3c^2$. Tie to check: wien_peak_x() $\approx2.8214$, the numerical $\int u\,d\omega$ equals radiation_energy_density(T), and stefan_boltzmann_constant() $\approx5.670\times10^{-8}$.

Mean-field theory replaces each of the $q$ neighbors of a spin $s_i$ by the average $\langle s\rangle=m$, so $-J\sum_{\langle ij\rangle}s_is_j$ acting on $s_i$ becomes an effective field $h_{\rm eff}=qJm+\mu B$. A single Ising spin in that field has

Self-consistency gives $m=\tanh[(T_c/T)m+b]$ with $T_c\equiv qJ/k$ and $b\equiv\mu B/kT$. At $b=0$ and small $m$, $\tanh u\approx u$ gives $m\approx(T_c/T)m$; a nonzero branch can appear only when the slope reaches $1$, i.e. $T=T_c=qJ/k$. Tie to check: critical_temperature(4, J) $=4J/k$, and mean_field_residual at the solved $m$ is $\approx0$.

Solve $m=\tanh(am)$, $a=T_c/T$, by intersecting the line $y=m$ with the curve $y=\tanh(am)$. Both pass through the origin, where the curve has slope $a$. For $a\le1$ ($T\ge T_c$) the curve starts no steeper than the line and, being concave for $m>0$, stays below it — so $m=0$ is the only root. For $a>1$ ($T<T_c$) the curve leaves the origin steeper than the line yet is bounded by $1$, so it must bend back and cross $y=m$ again at some

the spontaneous magnetization. As $T\to0$, $a=T_c/T\to\infty$ and $\tanh(am)\to1$ for any fixed $m>0$, forcing saturation:

Tie to check: spontaneous_magnetization(330,300) $=0$ (paramagnet), (150,300) $\approx0.957$, and (5,300) $>0.99\to1$.

Near $T_c$ the magnetization is small, so expand $m=\tanh(am)$, $a=T_c/T$, using $\tanh u=u-\tfrac13u^3+\cdots$:

Just below $T_c$, $a=T_c/T\to1^+$, so $a^3\to1$ and $a-1=\dfrac{T_c-T}{T}\approx\dfrac{T_c-T}{T_c}=1-\dfrac{T}{T_c}$. Therefore

so the order-parameter exponent in $m\sim(T_c-T)^{\beta}$ is $\beta=\tfrac12$ (mean field). Tie to check: critical_exponent_beta(Tc, "meanfield") $\approx0.5$, read off as the log–log slope of $m$ vs $(T_c-T)$.

The exact transfer-matrix magnetization per spin is

Switch off the field: $\sinh(\beta\cdot0)=0$, so the numerator vanishes, while the denominator becomes $\sqrt{0+e^{-4\beta J}}=e^{-2\beta J}$, strictly positive for every finite $T$ (finite $\beta$). Hence

no spontaneous magnetization — thermal fluctuations destroy long-range order in one dimension. Order sets in only as $T\to0$ ($\beta\to\infty$, $e^{-2\beta J}\to0$) or when a field $h\neq0$ makes $\sinh(\beta h)>0$. Tie to check: ising_1d_magnetization(T, J, 0) $=0$ at $T=10,100,500$ K, but nonzero once a field is applied.

Minimize $F(m)=a(T-T_c)m^2+bm^4$ ($a,b>0$):

The roots are $m=0$ and $m^2=\dfrac{a(T_c-T)}{2b}$. For $T>T_c$ the bracket is positive for all $m$, so $m=0$ is the only extremum — a single well. For $T<T_c$ the factor $T_c-T>0$ makes a real pair appear,

the double-well minima; meanwhile $m=0$ becomes a maximum since $\partial^2F/\partial m^2\big|_0=2a(T-T_c)<0$. The exponent is again $\beta=\tfrac12$. Tie to check: below $T_c$, landau_free_energy(m0,…) $<$ landau_free_energy(0,…), and since $m_0\propto\sqrt{T_c-T}$, doubling $(T_c-T)$ scales $m_0$ by $\sqrt2$.

Set $\alpha=\tfrac{m}{2kT}$, so $f(v)=4\pi(\alpha/\pi)^{3/2}v^2e^{-\alpha v^2}$. The peak satisfies $\frac{d}{dv}\big(v^2e^{-\alpha v^2}\big)=0$, i.e. $2v=2\alpha v^3$, giving

The mean and mean-square use $\int_0^\infty v^3e^{-\alpha v^2}\,dv=\tfrac{1}{2\alpha^2}$ and $\int_0^\infty v^4e^{-\alpha v^2}\,dv=\tfrac{3}{8}\sqrt\pi\,\alpha^{-5/2}$:

Dividing all three by $v_p$ gives $1:\sqrt{4/\pi}:\sqrt{3/2}$, with $kT/m$ cancelling so the ratio is $T$-independent. For N$_2$ ($m=4.652\times10^{-26}$ kg) at $300$ K they evaluate to $422,\,476,\,517$ m/s, matching most_probable_speed, mean_speed, rms_speed.

With $\alpha=\tfrac{m}{2kT}$ the prefactor is $4\pi(\alpha/\pi)^{3/2}$, and the Gaussian moments are $\int_0^\infty v^2e^{-\alpha v^2}\,dv=\tfrac{\sqrt\pi}{4}\alpha^{-3/2}$ and $\int_0^\infty v^4e^{-\alpha v^2}\,dv=\tfrac{3\sqrt\pi}{8}\alpha^{-5/2}$. Normalization is then

and the second moment is

Hence $\langle\tfrac12 mv^2\rangle=\tfrac12 m\cdot\tfrac{3kT}{m}=\tfrac32 kT$ — equipartition, $\tfrac12 kT$ per translational degree of freedom. Numerically integrating maxwell_speed_pdf reproduces $\int f\,dv=1$, $\langle v^2\rangle=3kT/m$, and $\tfrac12 m\langle v^2\rangle=$ mean_kinetic_energy(T) $=6.21\times10^{-21}$ J at $300$ K.

Model the molecule as a disk of cross-section $\sigma=\pi d^2$ (two centres collide when their separation reaches $d$). In time $t$ it sweeps a volume $\sigma\langle v_\text{rel}\rangle t$ and strikes $n\sigma\langle v_\text{rel}\rangle t$ partners, so

Averaging two independent Maxwellians gives $\langle v_\text{rel}\rangle=\sqrt2\,\langle v\rangle$, so $\lambda=1/(\sqrt2\,n\sigma)$ and the collision frequency $z=\langle v\rangle/\lambda =\sqrt2\,n\sigma\langle v\rangle$. At $1$ atm and $300$ K, $n=P/kT\approx2.45\times10^{25}\,\text{m}^{-3}$ and $\sigma=\pi(3.7\times10^{-10})^2$:

These reproduce mean_free_path(n, 3.7e-10) $\approx67$ nm and collision_rate $=\langle v\rangle/\lambda$.

Place the hole at $x=0$; only molecules with $v_x>0$ escape, so the flux is

with $g$ the one-dimensional Maxwellian. The integral is $\int_0^\infty v_x g\,dv_x=\sqrt{kT/2\pi m}$, and since $\langle v\rangle=\sqrt{8kT/\pi m}$ this equals $\tfrac14\langle v\rangle$, giving

Because $\langle v\rangle\propto m^{-1/2}$, at fixed $n,T$ the effusion rate scales as $1/\sqrt m$: a lighter gas effuses faster, with rate ratio $\sqrt{m_2/m_1}$ (Graham's law). This is exactly effusion_flux(n, m, T) $=\tfrac14 n\,$mean_speed(m, T).

Hold the Brownian particles in equilibrium under a weak force $F=-dU/dx$. Two currents balance: a drift current $j_\text{drift}=n\mu F$ (mobility $\mu$ set by the drift speed $v=\mu F$) and Fick diffusion $j_\text{diff}=-D\,dn/dx$. Zero net flux requires

Equilibrium also fixes the barometric law $n\propto e^{-U/kT}$, so $\frac{dn}{dx}=-\frac{1}{kT}\frac{dU}{dx}\,n=\frac{F}{kT}n$. Substituting,

Diffusion (fluctuation) is tied to mobility (the response to drag) through $T$ — a fluctuation-dissipation theorem. Thus einstein_relation_diffusion(mu, T) $=\mu kT$ is linear in both $\mu$ and $T$ (doubling $T$ doubles $D$); at $\mu=10^{11}$ s/kg, $T=300$ K it returns $4.14\times10^{-10}\,\text{m}^2/\text{s}$.