We use two very large cosmological simulations to study how the density profiles of relaxed Λ cold dark matter dark haloes depend on redshift and on halo mass. We confirm that these profiles deviate slightly but systematically from the NFW form and are better approximated by the empirical formula, d log ρ/d log r∝rα, first used by Einasto to fit star counts in the Milky Way. The best-fitting value of the additional shape parameter, α, increases gradually with mass, from α∼ 0.16 for present-day galaxy haloes to α∼ 0.3 for the rarest and most massive clusters. Halo concentrations depend only weakly on mass at z= 0, and this dependence weakens further at earlier times. At z∼ 3 the average concentration of relaxed haloes does not vary appreciably over the mass range accessible to our simulations (M≳ 3 × 1011 h−1 M⊙). Furthermore, in our biggest simulation, the average concentration of the most massive, relaxed haloes is constant at 〈c200〉∼ 3.5–4 for 0 ≤z≤ 3. These results agree well with those of Zhao et al. and support the idea that halo densities reflect the density of the universe at the time they formed, as proposed by Navarro, Frenk & White. With their original parameters, the NFW prescription overpredicts halo concentrations at high redshift. This shortcoming can be reduced by modifying the definition of halo formation time, although the evolution of the concentrations of Milky Way mass haloes is still not reproduced well. In contrast, the much-used revisions of the NFW prescription by Bullock et al. and Eke, Navarro & Steinmetz predict a steeper drop in concentration at the highest masses and stronger evolution with redshift than are compatible with our numerical data. Modifying the parameters of these models can reduce the discrepancy at high masses, but the overly rapid redshift evolution remains. These results have important implications for currently planned surveys of distant clusters.