Average entropy of quantum system in bipartite pure state for finite temperature

[I got halfway through writing this when I found the paper that answers the question in (essentially) the affirmative. I'll post it anyways in case anyone is interested.]

Background: If a random vector $\psi$ (pure quantum state) is drawn from a $NM$-dimensional vector space for finite $N$ according to the Haar measure, then the entanglement entropy $$S(\rho_N) = \mathrm{Tr}[\rho_N \mathrm{log} \rho_N], \qquad \rho_N = \mathrm{Tr}_M[\psi\psi^\dagger]$$ across a tensor decomposition into $N$ and $M$ dimensional vector spaces (subsystems) is highly likely to be almost the maximum $$S_{\mathrm{max}} = \mathrm{log}_2(\mathrm{min}(N,M)) \,\, \mathrm{bits},$$ for any such choice of decomposition. More precisely, if we fix $\alpha=M/N$ and let $N\to \infty$, then the fraction of the Haar volume of states that have entanglement entropy more than an exponentially small (in $N$) amount away from the maximum is suppressed exponentially (in $N$). In physics this was known as Page's conjecture (later proved), and for mathematicians it is a simple consequence of the concentration of measure phenomenon.

For any given Hermitian operator $H$ on the vector space acting as a Hamiltonian, we can assign an expected value of the energy $\langle H \rangle_\psi = (\psi^\dagger, H \psi)$ to a given vector $\psi$. The Haar measure, then, can be considered a Gibbs probability distribution $p_\psi \propto e^{-\beta\langle H \rangle_\psi}$ in the infinite temperature limit ($\beta \to 0$), i.e. all states are equally likely. The Page conjecture is exactly what you'd expect if a subsystem of a global pure state is maximally mixed, given its size.

My question: is there a way to extend this statement about typical entanglement entropies to cases where the probability distribution is one for finite temperature? More physically, can we prove that the typical entanglement entropy of a system in a pure global state is exponentially close to

$$S_\beta = \mathrm{Tr}[\rho^{(\beta)}_N \mathrm{log} \rho^{(\beta)}_N], \qquad \rho^{(\beta)}_N = \mathrm{Tr}_M[\psi\psi^\dagger]$$

with high probability, according to the Gibbs distribution $p_\psi \propto e^{-\beta\langle H \rangle_\psi}$?

The statement, eq. (2) in the article, is that for arbitrarily small $\epsilon$ there exist $$\eta = \epsilon +\frac{1}{2}\sqrt{\frac{d_s}{d_E^{\mathrm{eff}}}} \quad , \quad \eta' = 4e^{-C d_R \epsilon^2}$$ such that $$\frac{V[\{ \vert \phi\rangle \in \mathcal{H}_R \vert D(\rho_S(\phi),\Omega_S)\ge\eta\}]}{V[\{ \vert \phi\rangle \in \mathcal{H}_R\}]} \le \eta'.$$ Here, $D$ is the trace distance, $\mathcal{H}_R$ is the subspace of the global Hilbert space satisfying the constraint, $\rho_S$ is the reduced state of the system, $d_S$ and $d_R$ are the dimensions of the system and the constrained subspace, and $$d_E^{\mathrm{eff}} = \frac{1}{\mathrm{Tr} \Omega_E^2} \ge \frac{d_R}{d_S}$$ is the effective size of the environment. In these expression, $\Omega_S = \langle \rho_S \rangle_{\mathcal{H}_R}$ and $\Omega_E = \langle \rho_E \rangle_{\mathcal{H}_R}$ are the Haar-average reduced states conditional on the constraint.