Controllability and observability energy functions play a fundamental role in model order reduction and are inherently connected to optimal control problems. For linear dynamical systems the energy functions are known to be quadratic polynomials and various low-rank approximation techniques allow for computing them in a large-scale setting. For nonlinear problems computing the energy functions is significantly more challenging. In this paper, we investigate a special class of nonlinear systems that have a linear state and a polynomial output equation. We show that the energy functions of these systems are again polynomials and investigate under which conditions they can effectively be approximated using low-rank tensors. Further, we introduce a new perspective on the well-established balanced truncation method for linear systems which then readily generalizes to the nonlinear systems under consideration. This new perspective yields a novel energy-based model order reduction procedure that accurately captures the input-output behavior of linear systems with polynomial outputs via a low-dimensional reduced order model. We demonstrate the effectiveness of our approach via two numerical experiments.