  In this work, we consider the $\mathcal{H}_2$ optimal model reduction of
dynamical systems that are linear in the state equation and up to quadratic
nonlinearity in the output equation. As our primary theoretical contributions,
we derive gradients of the squared $\mathcal{H}_2$ system error with respect to
the reduced model quantities and, from the stationary points of these
gradients, introduce Gramian-based first-order necessary conditions for the
$\mathcal{H}_2$ optimal approximation of a linear quadratic output (LQO)
system. The resulting $\mathcal{H}_2$ optimality framework neatly generalizes
the analogous Gramian-based optimality framework for purely linear systems.
Computationally, we show how to enforce the necessary optimality conditions
using Petrov-Galerkin projection; the corresponding projection matrices are
obtained from a pair of Sylvester equations. Based on this result, we propose
an iteratively corrected algorithm for the $\mathcal{H}_2$ model reduction of
LQO systems, which we refer to as LQO-TSIA (linear quadratic output two-sided
iteration algorithm). Numerical examples are included to illustrate the
effectiveness of the proposed computational method against other existing
approaches.
