In this article, we investigate a balancing-based model order reduction method for dynamical systems with a linear dynamical equation and a quadratic output function. With this aim, we propose a new algebraic observability Gramian for the system based on the Hilbert space adjoint theory. We then show that the proposed Gramian satisfies a particular type of generalized Lyapunov equation and we investigate its connection to energy functionals. It allows one to find states that are hard to control and hard to observe via an appropriate balancing transformation. Truncation of such states yields reduced-order models. Finally, based on $\mathscr {H}_2$ energy considerations, we derive error bounds, depending on the neglected singular values. We demonstrate the efficiency of the proposed methodology using a numerical example.