The density operator of a quantum system in thermal equilibrium at temperature \(T\) is given by \[\rho_T = \frac{e^{-H / k_{\text{B}} T}}{Z} \,,\] where \(H\) is the system Hamiltonian, \(k_{\text{B}}\) is the Boltzmann constant, and \[Z = \operatorname{Tr} e^{-H / k_{\text{B}} T}\] is the partition function.

This expression for the density matrix can be obtained by extremizing the von Neumann entropy \[S = -k_{\text{B}} \operatorname{Tr} \rho \ln \rho\] over all states \(\rho\) of the same mean energy \[E = \operatorname{Tr} \rho H \,.\]

**Proof:** We want to extremize \(S\) subject to (i) the normalization constraint, \(\operatorname{Tr} \rho = 1\), and (ii) the mean energy constraint, \(\operatorname{Tr} \rho H = E\). To this end, we introduce two Lagrange multipliers, \(\alpha k_{\text{B}}\) and \(\beta k_{\text{B}}\), and perform unconstrained extremization of
\[\mathcal{F} = -k_{\text{B}} \operatorname{Tr} \rho \ln \rho - \alpha k_{\text{B}} \left( \operatorname{Tr} \rho - 1 \right) - \beta k_{\text{B}} \left( \operatorname{Tr} \rho H - E \right) \,.\]
The variation with respect to \(\rho\) yields
\[\delta \mathcal{F} = -k_{\text{B}} \operatorname{Tr} \Big\{ \left( \ln \rho + 1 + \alpha + \beta H \right) \delta \rho \Big\} \,.\]
Functional \(\mathcal{F}\) is stationary at \(\rho = \rho_T\), given by
\[\ln \rho_T + 1 + \alpha + \beta H = 0 \,,\]
or
\[\rho_T = e^{-1 - \alpha - \beta H} \,.\]
The normalization condition, \(\operatorname{Tr} \rho = 1\), requires
\[e^{-1-\alpha} \operatorname{Tr} e^{-\beta H} = 1 \,.\]
This implies that
\[\rho_T = \frac{e^{-\beta H}}{\operatorname{Tr} e^{-\beta H}} \,.\]
The value of \(\beta\) is determined by the energy condition, \(\operatorname{Tr} \rho H = E\), which now reads
\[\frac{\operatorname{Tr} H e^{-\beta H}}{\operatorname{Tr} e^{-\beta H}} = E \,.\]

Finally, the relation between \(\beta\) and the temperature \(T\) can be obtained by calculating the variation of entropy \(S\) and mean energy \(E\) with respect to \(\rho\) around the equilibrium state. Considering \[\rho = \rho_T + \delta \rho \,,\] with \[\operatorname{Tr} \delta \rho = 0\] to ensure normalization, we obtain \[\delta S = -k_{\text{B}} \operatorname{Tr} \Big\{ \left( \ln \rho_T + 1 \right) \delta \rho \Big\} = k_{\text{B}} \beta \operatorname{Tr} H \delta \rho\] and \[\delta E = \operatorname{Tr} H \delta \rho \,.\] Hence, \[\frac{1}{T} \equiv \frac{\delta S}{\delta E} = k_{\text{B}} \beta \,,\] or \[\beta = \frac{1}{k_{\text{B}} T} \,.\]