  We generalize the concept of optical scattering matrix ($S$-matrix) to
characterize harmonic generation and frequency mixing in planar metasurfaces in
the limit of undepleted pump approximation. We show that the symmetry
properties of such nonlinear $S$-matrix are determined by the microscopic and
macroscopic symmetries of the metasurface. We demonstrate that for description
of degenerate frequency mixing processes such as optical harmonic generation,
the multidimensional $S$-matrix can be replaced with a reduced two-dimensional
$S$-matrix. We show that for metasurfaces possessing specific point group
symmetries, the selection rules determining the transformation of the reduced
nonlinear $S$-matrix are simplified substantially and can be expressed in a
compact form. We apply the developed approach to analyse chiral harmonic
generation in nonlinear metasurfaces with various symmetries including
rotational, in-plane mirror, and out-of-plane mirror symmetries. For each of
those symmetries, we confirm the results of the developed analysis by full-wave
numerical calculations. We believe our results provide a new paradigm for
engineering nonlinear optical properties of metasurfaces which may find
applications in active and nonlinear optics, biosensing, and quantum
information processing.
