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} \,.$