A black hole evaporates by Hawking radiation. Each mode of that radiation is thermal. If the total state is nevertheless to be pure, modes must be entangled. Estimating the minimum size of this entanglement has been an important outstanding issue. We develop a new theory of constrained random symplectic transformations, based on the assumptions that the total state is pure and Gaussian with given marginals. In the random constrained symplectic model we then compute the distribution of mode-mode correlations, from which we bound mode-mode entanglement. Modes of frequency much larger than [k(B)T(H)(t)/h] are not populated at time t and drop out of the analysis. Among other relatively thinly populated modes (earlytime high-frequency modes and/or late modes of any frequency), we find correlations and hence entanglement to be strongly suppressed. Relatively highly populated modes (early-time low-frequency modes) can, on the other hand, be strongly correlated, but a detailed analysis reveals that they are nevertheless very unlikely to be entangled. Our analysis hence establishes that restoring unitarity after a complete evaporation of a black hole does not require any significant quantum entanglement between any pair of Hawking modes. Our analysis further gives exact general expressions for the distribution of modemode correlations in random, pure, Gaussian states with given marginals, which may have applications beyond black hole physics.