Hint $\ a^2+ab+b^2 \ $ is the norm of an Eisenstein integer $\,a-b\,\zeta,\,$ so the multiplicativity of integers of this form follows from the $\rm\color{#c00}{multiplicativity}$ of the norm of Eisenstein integers, i.e.
$$ \color{#c00}{N(\alpha)N(\beta)} = \alpha\bar\alpha\,\beta\bar\beta = \alpha\beta\,\bar\alpha\bar\beta = \alpha\beta\,\overline{\alpha\beta} = \color{#c00}{N(\alpha\beta)}\qquad$$
It is an Eisenstein-integer analog of the well-known Brahmagupta–Fibonacci identity below, which arises from norm multiplicativity for quadratic integers $\,a+b\sqrt{- n}\in \Bbb Z[\sqrt{-n}],\,$ i.e.
$$\ \ \begin{eqnarray} N(\alpha)N(\beta) &=& N(\alpha\beta)\\[.2em]
(a_1^2 + nb_1^2)\ (a_2^2 + nb_2^2) & =& \!(a_1a_2\pm nb_1b_2)^2 + n(a_1b_2\mp a_2b_1)^2 \end{eqnarray}$$
The upper $\,\pm\,$ signs are from $\ N(\alpha)\:\!N(\bar\beta) = N(\alpha\bar \beta)\:$.
Because the norm map is multiplicative it preserves many properties related to factorization. For example, in many favorable contexts (e.g. Galois) a number ring enjoys unique factorization iff its monoid of norms does. For further interesting examples see the papers of Bumby, Dade, Lettl, Coykendall cited in this answer.