Solution (source code)

= Solution

Take the two magnetic moments along $z$ to be $\pm\mu$, so $H=-\mu B_z\sigma_z$ and the two <energies> are $\mp\mu B_z$. The <quantum Gibbs state> at inverse <temperature> $\beta=1/(k_BT)$ is
$$
\boxed{\rho_E=\frac{e^{\beta\mu B_z\sigma_z}}{2\cosh(\beta\mu B_z)}
=\frac1{2\cosh(\beta\mu B_z)}\begin{pmatrix}e^{\beta\mu B_z}&0\\0&e^{-\beta\mu B_z}\end{pmatrix}
=\frac{I+\tanh(\beta\mu B_z)\sigma_z}{2}.}
$$
Its <Bloch vector> is $(0,0,\tanh\beta\mu B_z)$. At $B_z=0$, or in the high-temperature limit, it approaches $I/2$; at $T\to0^+$ with nonzero $\mu B_z$, it approaches the lower-energy <pure state>. Here $\mu$ is the magnitude or signed coefficient multiplying $\sigma_z$ in the magnetic <Hamiltonian>, rather than a moment with an additional unstated factor of two.