Direct numerical simulations and a self-similar analysis of the single-fluid Boussinesq Rayleigh–Taylor instability and transition to turbulence are used to investigate Rayleigh–Taylor turbulence. The Schmidt, Atwood and bulk Reynolds numbers are $Sc\,{=}\,1$, $A\,{=}\,0.01$, $Re \,{\le}\, 3000$. High-Reynolds-number moment self-similarity, consistent with the the energy cascade interpretation of dissipation, is used to analyse the DNS results. The mixing layer width obeys a differential equation with solution $h(t;C_o,h_0)\,{=}\,\fourth C_o Agt^2+ \sqrt{AgC_o}h^{1/2}_0 t+h_0$; the result for $h(t;C_o,h_0)$ is a rigorous consequence of only one ansatz, self-similarity. It indicates an intermediate time regime in which the growth is linear and the importance of a virtual origin. At long time the well-known $h \sim \fourth C_o Agt^2$ scaling dominates. The self-similar analysis indicates that the asymptotic growth rate is not universal. The scalings of the second-order moments, their dissipations, and production–dissipation ratios, are obtained and compared to the DNS. The flow is not self-similar in a conventional sense – there is no single length scale that scales the flow. The moment similarity method produces three different scalings for the turbulence energy-containing length scale, $\ell$, the Taylor microscale, $\la$, and the Kolmogorov dissipation scale, $\eta$. The DNS and the self-similar analysis are in accord showing $\ell \,{\sim}\, Agt^2$, $\la \,{\sim}\, t^{1/2}$ and $\eta \,{\sim}\, (({A^2g^2}/{\nu^3})t)^{-1/4}$ achieving self-similar behaviour within three initial eddy turnovers of the inception of the turbulence growth phase at bulk Reynolds numbers in the range of ${\it Re}\,=\,800$–1000 depending on initial conditions. A picture of a turbulence in which the largest scales grow, asymptotically, as $t^2$ and the smallest scales decrease as $t^{-1/4}$, emerges. As a consequence the bandwidth of the turbulence spectrum grows as $t^{9/4}$ and is consistent with the $R_t^{3/4}$ Kolmogorov scaling law of fully developed stationary turbulent flows. While not all moments are consistent, especially the dissipations and higher-order moments in the edge regions, with the self-similar results it appears possible to conclude that: (i) the turbulence length scales evolve as a power of $h(t;C_o,h_0)$; (ii) $\al$, as demonstrated mathematically for self-similar Rayleigh–Taylor turbulence and numerically by the DNS, is not a universal constant; (iii) there is statistically significant correlation between decreasing $\alpha$ and lower low-wavenumber loading of the initial spectrum.