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 multiplicativity of the norm of Eisenstein integers, i.e.
$$ N(\alpha)N(\beta) = \alpha\bar\alpha\,\beta\bar\beta = \alpha\beta\,\bar\alpha\bar\beta = \alpha\beta\,\overline{\alpha\beta} = N(\alpha\beta)$$
Remark $\ $ This is an Eisenstein-integer analog of the Brahmagupta–Fibonacci identity for composition of sums of squares, which follows by multiplicativity of Gaussian integer norms
$$\ \begin{eqnarray} N(\alpha)N(\beta) &=& N(\alpha\beta)\\ (a_1^2 + b_1^2)\ (a_2^2 + b_2^2) & =& \!(a_1a_2\pm b_1b_2)^2 + ( 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)\:$.