How do we find the canonical ensemble density matrix for two spins? A compound system is constructed by two coupling spins, and the Hamiltonian is 

$$
\hat H = -J\hat\sigma_1·\hat\sigma_2 - \mu_\mathbf{B}\big( \hat\sigma_{1z}+\hat\sigma_{2z} \big)B.
$$

So, how can we obtain the canonical ensemble density matrix, i.e. $\hat\rho = Z^{-1}\exp(-\beta\hat H)$? Can anyone give me some hints or references?
$$EDIT: Question two$$
Suppose that the first spin represents system while the other represents environment. Their Hamiltonian are $\hat H_1 =-\mu_\mathbf{B} \hat\sigma_{1z}B$ and $\hat H_2 =-\mu_\mathbf{B} \hat\sigma_{2z}B$. Then $\hat H_{12} = -J\hat\sigma_1·\hat\sigma_2$ represents the interaction between system and env. My question is that what's the reduced density matrix of the first spin?
 A: There are four state kets $|\uparrow_1\uparrow_2\rangle,|\uparrow_1\downarrow_2\rangle,|\downarrow_1\uparrow_2\rangle,|\downarrow_1\downarrow_2\rangle$ representing four possible configuration of the two spins, the eigenvectors of the Hamiltonian is the linear combination of them. There will be four eigenvectors and four eigenenergy.
You can also obtain the eigen-energy in the matrix way. Just write out the matrix form of the Hamiltonian, note that $\sigma_1\cdot\sigma_2=\sigma_{1x}\sigma_{2x}+\sigma_{1y}\sigma_{2y}+\sigma_{1z}\sigma_{2z}$ should be understood as the direct product of matrix. Similarly, $\mu B \sigma_{1z}$ should be understand as $\sigma_{1z}\otimes I_2$.
Once you get the eigenvalue, your matrix form of $\hat \rho$ is a diagonal matrix $[e^{-\beta E_1}/Z,e^{-\beta E_2}/Z,e^{-\beta E_3}/Z,e^{-\beta E_4}/Z]$
Regarding your second question
Once you obtained the eigenvalue $E_1,E_2,\cdots$ and corresponding eigenvector $|1\rangle,|2\rangle\cdots$, your density matrix operator is written as $\hat\rho=e^{-\beta E_1}/Z |1\rangle\langle1|+e^{-\beta E_2}/Z |2\rangle\langle2|+\cdots$ 
Expand $|1\rangle,|2\rangle\cdots$ using the four vector $|\uparrow_1\uparrow_2\rangle,|\uparrow_1\downarrow_2\rangle,|\downarrow_1\uparrow_2\rangle,|\downarrow_1\downarrow_2\rangle$, you can rewrite your density matrix $\hat\rho$. The reduced density matrix $\hat\rho_r=\mathrm{Tr}_2(\hat\rho)=\langle\uparrow_2|\hat\rho|\uparrow_2\rangle+\langle\downarrow_2|\hat\rho|\downarrow_2\rangle$
