  The Perdew-Zunger (PZ) self-interaction correction (SIC) is an established
tool to correct unphysical behavior in density functional approximations. Yet,
PZ-SIC is well-known to sometimes break molecular symmetries. An example of
this is the benzene molecule, for which PZ-SIC predicts a symmetry-broken
electron density and molecular geometry, since the method does not describe the
two possible Kekul\'e structures on an even footing, leading to local minima
[Lehtola et al, J. Chem. Theory Comput. 2016, 12, 3195]. PZ-SIC is often
implemented with Fermi-L\"owdin orbitals (FLOs), yielding the FLO-SIC method,
which likewise has issues with symmetry breaking and local minima [Trepte et
al, J. Chem. Phys. 2021, 155, 224109].
  In this work, we propose a generalization of PZ-SIC - the ensemble PZ-SIC
(E-PZ-SIC) method - which shares the asymptotic computational scaling of PZ-SIC
(albeit with an additional prefactor). E-PZ-SIC is straightforwardly applicable
to various molecules, merely requiring one to average the self-interaction
correction over all possible Kekul\'e structures, in line with chemical
intuition. We showcase the implementation of E-PZ-SIC with FLOs, as the
resulting E-FLO-SIC method is easy to realize on top of an existing
implementation of FLO-SIC. We show that E-FLO-SIC indeed eliminates symmetry
breaking, reproducing a symmetric electron density and molecular geometry for
benzene. The ensemble approach suggested herein could also be employed within
approximate or locally scaled variants of PZ-SIC and their FLO-SIC versions.
