Approaching the fixed-point problem an SCF solves from the direct minimisation (DM) point of view (for the relation of DM to SCF, see my previous post, the key quantities to check for convergence are actually the total energy and it's derivative, i.e. the occupied-virtual block of the Fock / Kohn-Sham matrix. The latter is after all the gradient and that's what you want to drive to zero in a fixed point. In an SCF context you could theoretically do the same and to the best of my knowledge some codes like PySCF actually use these two criteria for convergence.
Of course there are other choices, which can be done and using the density works the same, since a zero density change implies a zero change in the Fock matrix and vice versa. From that point of view the density is kind of a substitute for checking convergence in the gradient, but of course not the same threshold can be used for both checks and tolerances need to be adjusted between them.
Since clearly the gradient being zero is the whole aim of an SCF, I find it dubious to only check for convergence in the energy as it could accidentally stagnate without being at a fixed point. In turn if you check for the density / gradient change being zero you should be safe. However, as @Tyberius already pointed out, computing the energy is pretty much for free if you compute density and Fock matrix and therefore performing a check on the energy does not cost and shouldn't do harm either.
Talking about convergence speeds: In absolute numbers the energy always converges first, since it depends quadratically on the density (so error 1e-3 in the density usually translates to 1e-6 in the energy). Therefore convergence criteria on energy should be the square of the convergence tolerance in the density / gradient. But notice that many codes are a bit sloppy here and do not print the norm of the density / gradient change in their outputs, but the norm squared, such that both numbers are again on the same scale (e.g. ABINIT does this). So one has to be a little careful when looking at the numbers presented from codes.
Let me point out one last subtlety from a mathematical point of view. When checking the convergence in density / gradient one is faced with the problem that these quantities are vectors or matrices, but of course we want to check convergence in a number. So we need to take the norm of the difference. But which norm should we choose? Most codes to the best of my knowledge use the Frobenius or l2 norms (so just square the elements, add them up and take the square root), but this is just one choice. There are plenty of other norms, see the wikipedia articles on the Lp norms and the matrix norms to get an idea. Which norm is the best one to choose depends a little on which property you are after in your calculation: Total energy, forces / gradients wrt. nuclear positions, partial charges etc, just because these ask different questions about the mathematical properties of the wavefunction.
Now, the take-away from this is not that one should use a different norm depending on what property one is after. This is just impractical and also in calculations with their finite basis set sizes all norms are equivalent up to a constant. But this constant might not be small and usually depends on the size of the basis (it typically grows with larger bases). So depending on what property you are after in the end, the number of digits you can trust in the computed answer differs. For example, even if the density is converged to 6 digits in the Frobenius norm, the forces might only be correct to 5 and the partial charges to 4.
