Forty-five years after the point de départ [Hohenberg and Kohn, Phys Rev, 1964, 136, B864] of density functional theory, its applications in chemistry and the study of electronic structures keep steadily growing. However, the precise form of the energy functional in terms of the electron density still eludes us-and possibly will do so forever [Schuch and Verstraete, Nat Phys, 2009, 5, 732]. In what follows we examine a formulation in the same spirit with phase space variables. The validity of Hohenberg-Kohn-Levy-type theorems on phase space is recalled. We study the representability problem for reduced Wigner functions, and proceed to analyze properties of the new functional. Along the way, new results on states in the phase space formalism of quantum mechanics are established. Natural Wigner orbital theory is developed in depth, with the final aim of constructing accurate correlation-exchange functionals on phase space. A new proof of the overbinding property of the Müller functional is given. This exact theory supplies its home at long last to that illustrious ancestor, the Thomas-Fermi model. 2 Even before the results in Ref.[4] the problem has been regarded as well nigh intractable. Nevertheless, there are the formal solution by Garrod and Percus [5], still an appropriate reference for N-representability; Ayers' reformulation of the latter [6]; and other recent remarkable progress both on it and on pending questions of the representability problem for 1-matrices [7][8][9][10]. We return to the Ayers method in the next section. 4 We still find most illuminating the treatment of the relations between phase space representations given by Cahill and Glauber [29] long ago.Also the form factors are easily expressed in terms of d 1 . As dequantization and reduction commute, we may also use formulas analogous to (3) for reduced density matrices. At each order, the reduced Wigner distributions contain the same information as the corresponding reduced density matrices. It is convenient to have "inversion formulas" to recover the latter matrices from the Wigner functions. It