Thermal state of light with two modes The thermal state corresponding to a photon mode $a$ can be written as 
$\rho_{a}(\beta) = \frac{e^{-\beta a^\dagger a}}{Tr(e^{-\beta a^\dagger a})}$. Since it is easy (I suppose) to talk in terms of the photon number state, I am interested in following:
1) How can one write a similar state for two modes $a_1$ and $a_2$?, That is, what would be $\rho_{a_1a_2}(\beta)?$
2) How can these states (one mode and multimode) be represented in terms of photon number states $|n>$?
 A: A general thermal state has the form
$$
\rho = \frac{e^{-\beta H} }{Z} 
$$
where $H$ is the Hamiltonian and $Z = \mathrm{Tr}\, e^{-\beta H}$ is the partition function.
For a single photon mode 
$$
H = \omega a^\dagger a
$$
which, for $\omega = 1$ gives the expression in your question
For a system with 2 modes the Hamiltonian is simply the sum of the Hamiltonians for each mode so
$$
H =\omega_1 a_1^\dagger a_1 + \omega_2 a_2^\dagger a_2
$$
and so 
$$ \rho = \frac{e^{-\beta(\omega_1 a_1^\dagger a_1 + \omega_2 a_2^\dagger a_2)}}{Z} 
$$
To express this in terms of the photon number basis we place the resolution of the identity, $\sum_{n_1, n_2}|n_1,n_2\rangle\langle n_1,n_2|$ before and after the density matrix
\begin{align}
\rho &= \frac{1}{Z}\sum_{n_1, n_2, n_1',n_2'}|n_1,n_2\rangle\langle n_1,n_2| e^{-\beta(\omega_1 a_1^\dagger a_1 + \omega_2 a_2^\dagger a_2)}|n_1',n_2'\rangle\langle n_1',n_2'| \\
&= \frac{1}{Z}\sum_{n_1, n_2}|n_1,n_2\rangle e^{-\beta(\omega_1n_1+\omega_2n_2)}\langle n_1,n_2|
\end{align}
