I may find a reference, but the derivation is simple once you know that for large ${\cal N}=2^N$ a random vector $X\in\mathbb{C}^{\cal N}$ of unit length has independently distributed complex Gaussian elements of zero mean and variance $1/{\cal N}$. Then you just perform the Gaussian averages to obtain the first and second moment of the energy $E=\langle X|\hat{H}|X\rangle$ in the state $X$,
$$\overline{E}=\mathbb{E}[\langle X|\hat{H}|X\rangle] = {\cal N}^{-1} \,{\rm tr}\, \hat{H}$$
$$\overline{E^2}=\mathbb{E}[\langle X|\hat{H}|X\rangle^2] = {\cal N}^{-2} (\,{\rm tr}\, \hat{H})^2 + {\cal N}^{-2} \,{\rm tr}\, \hat{H}^2$$
$$\Rightarrow {\rm var}\,E={\cal N}^{-2} \,{\rm tr}\, \hat{H}^2=4^{-N}\|\hat{H}\|_F.$$