distribution $N_1$ of Eq. (1) in the main text to second order and the transmission function to first order, and then substituting into Eq. (2), and setting the temperatures $T_1 = T_2 = T$, we obtain $$\begin{aligned}
H &=&  \int_{-\infty}^{+\infty}\!\! \frac{d\omega}{4\pi} \hbar \omega 
\int \frac{d^2 {\bf k}}{(2\pi)^2}  T_{12}({\bf k}, \omega)\Big|_{v_1=0} \frac{1}{2} \frac{\partial^2 N}{\partial \omega^2} 
\bigl( \delta  \omega \bigr)^2 \nonumber \\
&+& \int_{-\infty}^{+\infty}\!\! \frac{d\omega}{4\pi} \hbar \omega 
\int \frac{d^2 {\bf k}}{(2\pi)^2}  \frac{\partial T_{12}({\bf k}, \omega)}{\partial v_1}\Big|_{v_1 =0}   \frac{\partial N}{\partial \omega} 
\bigl( v_1 \delta  \omega \bigr) \nonumber.\end{aligned}$$ The first order term proportional to $\delta \omega$ is 0 because it is an odd function in $k_x$. We will write this long expression as $H = a\, v_1^2 = (a_1 +a_2) v_1^2$ where $a_1$ is from the $(\delta \omega)^2$ term and $a_2$ from the cross term. The coefficient $a_2$ is complicated, hence we focus only on $a_1$. Since the integral over $k_x^2$ and $k_y^2$ factor is the same, we symmetrize the formula about $x$ and $y$ and divide by two. With further simplification of the second derivation of the Bose function, we obtain $$a_1 = \int_{-\infty}^{+\infty} \frac{d\omega}{4\pi} \int \frac{d^2 {\bf k}}{(2\pi)^2}
 \frac{ T_{12}({\bf k}, \omega)\, k^2 (\beta \hbar)^2 \hbar \omega }{ 16 \sinh^2(\frac{\beta \hbar \omega}{2}) 
