  Recent gravitational wave observations show evidence for the presence of
higher harmonics, thus possibly indicating that these waves were generated in
the inspiral of compact objects with asymmetric mass ratios. Signals with
higher harmonics contain a trove of information that can lead to a better
estimation of system parameters and possibly to more stringent tests of general
relativity. Gravitational wave model that include higher harmonics, however,
have only been developed within general relativity, while models to test
theory-agnostic deviations from general relativity have been purely based on
the signal's dominant mode. We here extend the parameterized post-Einsteinian
framework to include the $\ell=2, 3$ and $4$ higher harmonics to first
post-Newtonian order, therefore providing a ready-to-use Fourier-domain
waveform model for tests of general relativity with higher harmonics. We find
that the deformations to the higher harmonics of the Fourier phase can be
easily mapped to the deformation of the dominant harmonic, while the
deformations to the higher-harmonics of the Fourier amplitude in general cannot
in a theory-agnostic way. Nonetheless, we develop a simple ansatz for the
deformations of the waveform amplitude (through a re-scaling deformation of the
time-domain amplitude) that both minimizes the number of independent amplitude
deformations parameters and captures the predictions of all known modified
theories to date.
