The homogenization of a heterogeneous mixture of two pure fluids with different densities by molecular diffusion and stirring induced by buoyancy-generated motions, as occurs in the Rayleigh–Taylor (RT) instability, is studied using direct numerical simulations. The Schmidt number, Sc, is varied by a factor of 20, 0.1 ≤ Sc ≤ 2.0, and the Atwood number, A, by a factor of 10, 0.05 ≤ A ≤ 0.5. Initial-density intensities are as high as 50% of the mean density. As a consequence of differential accelerations experienced by the two fluids, substantial and important differences between the mixing in a variable-density flow, as compared to the Boussinesq approximation, are observed. In short, the pure heavy fluid mixes more slowly than the pure light fluid: an initially symmetric double delta density probability density function (PDF) is rapidly skewed and, only at long times and low density fluctuations, does it relax to a Gaussian-like PDF. The heavy–light fluid mixing process asymmetry is relevant to the nature of molecular mixing on different sides of a high-Atwood-number RT layer. Diverse mix metrics are used to examine the homogenization of the two fluids. The conventional mix parameter, θ, is mathematically related to the variance of the excess reactant of a hypothetical fast chemical reaction. Bounds relating θ and the normalized product, Ξ, are derived. It is shown that θ underpredicts the mixing, as compared to Ξ, in the central regions of an RT layer; in the edge regions, θ is larger than Ξ. The shape of the density PDF cannot be inferred from the usual mix metrics popular in applications. For example, when θ, Ξ ≥ 0.6, characteristic of the interior of a fully developed RT layer, the PDFs can have vastly different shapes. Bounds on the fluid composition using two low-order moments of the density PDF are derived. The bounds can be used as realizability conditions for low-dimensional models. For the measures studied, the tightest bounds are obtained using Ξ and mean density. The structure of the flow is also examined. It is found that, at early times, the buoyancy production term in the spectral kinetic energy equation is important at all wavenumbers and leads to anisotropy at all scales of motion. At later times, the anisotropy is confined to the largest and smallest scales: the intermediate scales are more isotropic than the small scales. In the viscous range, there is a cancellation between the viscous and nonlinear effects, and the buoyancy production leads to a persistent small-scale anisotropy.