arXiv:2610.37387v1 Announce Type: new Abstract: We develop an $\mathcal{H}_2$ model reduction framework for polynomial bilinear systems with quadratic outputs. We define the associated $\mathcal{H}_2$ norm and show its consistency with the established norm in the original (non-polynomial) BQO case. We derive first-order necessary conditions for $\mathcal{H}_2$ optimal model reduction of these polynomial systems and use them to construct an efficient iterative algorithm for computing reduced BQO models that approximately satisfy these conditions. Numerical results on two examples demonstrate the accuracy and effectiveness of the proposed reduction approach.