Motivated by the ever-growing demand for \emph{green} wireless communications and the advantages of \emph{cell-free} (CF) massive multiple-input multiple-output (MIMO) systems, we focus on the design of their downlink for optimal \emph{energy efficiency} (EE). To address this fundamental topic, we assume that each access point (AP) is deployed with multiple antennas and serves multiple users on the same time-frequency resource while the APs are Poisson point process (PPP) distributed, which approaches realistically their opportunistic spatial randomness. Relied on tools from stochastic geometry, we derive a lower bound on the downlink average achievable spectral efficiency (SE). Next, we consider a realistic power consumption model for CF massive MIMO systems. These steps enable the formulation of a tractable optimization problem concerning the downlink EE per unit area, which results in the analytical determination of the optimal pilot reuse factor, the AP density, and the number of AP antennas and users that maximize the EE. Notably, the EE per unit area and not just the EE is the necessary metric to describe CF systems, where we meet multi-point transmission. Hence, we provide useful design guidelines for CF massive MIMO systems relating to fundamental system variables towards optimal EE. Among the results, we observe that a lower pilot reuse factor enables a decrease of the interference, and subsequently, higher EE up to a specific value. Overall, it is shown that the CF massive MIMO technology is a promising candidate for next-generation networks achieving simultaneously high SE and EE per unit area.