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$
