  We construct nuclear energy density functionals in terms of derivatives of
densities up to sixth, next-to-next-to-next-to-leading order (N3LO). A
phenomenological functional built in this way conforms to the ideas of the
density matrix expansion and is rooted in the expansions characteristic to
effective theories. It builds on the standard functionals related to the
contact and Skyrme forces, which constitute the zero-order (LO) and
second-order (NLO) expansions, respectively. At N3LO, the full functional with
density-independent coupling constants, and with the isospin degree of freedom
taken into account, contains 376 terms, while the functionals restricted by the
Galilean and gauge symmetries contain 100 and 42 terms, respectively. For
functionals additionally restricted by the spherical, space-inversion, and
time-reversal symmetries, the corresponding numbers of terms are equal to 100,
60, and 22, respectively.
