We can explicitly give the requested implicit definition for multiplication: It is the unique function on $(\mathbb{N},S,0)$ satisfying: \begin{align} 0a&=0\\ (S0)a&=a\\ ab&=ba\\ a(bc)&=(ab)c\\ (Sa)b=(Sa)c &\to b=c\\ (Sa)Sa &= S(aSSa)\\ (SS0)Sa &= SS(aSS0)\\ (S(aSS0))(S(ab))&= S(aSS(bS(aSS0))) \end{align} The last three identities would be more familiarly written as \begin{align} (a+1)^2 &= 1+a(2+a)\\ 2(1+a)&=2+2a\\ (1+2a)(1+ab)&=1+a(2+b(1+2a)) \end{align} We refer to them as the square, even and odd identities: they give the key facts about squares and multiplication by even and odd numbers. The square identity defines the square of $a+1$ in terms of multiplication by $a$, the odd identity defines (a special case of) multiplication by $2a+1$ in terms of multiplication by $a$, and we can use them in inductions.
(The odd identity is a special case of $1+a(1+b(1+c(1+a(1+b))))=(1+a+ab)(1+abc)$, which may be easier to understand.)
As is usual in these matters, we look at numerals of the form $S^n0$. For each positive integer $n$, that numeral is a term in the language of the model. We quantify over $n$'s outside the model.
We prove by strong induction that for all $m$ and $n$, the axioms imply $S^m0\ S^n0=S^{mn}0$. This determines the values of the multiplication function for the whole domain of the model.
The cases $m=0$ and $m=1$ are immediate.
The case $m=2$ follows as an induction from the even identity.
The composite case $m=jk$ with $j<m$, $k<m$ is $$S^{jk}0\ S^n0 = S^j0\ S^k0\ S^n0 = S^j0\ S^{kn}0 = S^{jkn}0$$ which follows from the cases $m=j$, $m=k$ and $m=j$.
This leaves the case where $m$ is an odd prime, and we take $m=1+2j$.
-
- If $n=0,1,2$ or $n=jk$ then we use the same argument as for $m$, so we may assume wlog that $n$ is also an odd prime.
-
- If $m=n$, then $$S^m0\ S^m0 = S(S^{m-1}0\ S^{m+1}0)=S(S^{m^2-1}0)=S^{m^2}0$$ by the square identity, the inductive hypothesis for $m-1$, and simplifying notation. So we may also assume wlog that $n<m$.
-
- If $n=1+jk$ (the mod-1 case), then abbreviating $a=S^j0$, $b=S^k0$ gives \begin{align} S^{1+2j}0\ S^{1+jk}0 &=S(aSS0)\ S(ab) \ \ \ \text{by inductive hypotheses}\\ &=S(aSS(bS(aSS0))) \ \ \ \ \text{by the odd identity}\\ &=S(aSS(bS(S^{2j}0))) \ \ \ \ \text{by inductive hypothesis for }m=j\\ &=S(aSS(b(S^{1+2j}0))) \ \ \ \ \text{by simplifying notation}\\ &=S(aSS(S^{(1+2j)k}0)) \ \ \ \ \text{by inductive hypothesis for }m=k\\ &=S(aS^{2+(1+2j)k}0) \ \ \ \ \text{by simplifying notation}\\ &=S(S^{2j+(1+2j)jk}0) \ \ \ \ \text{by inductive hypothesis for }m=j\\ &=S^{(1+2j)(1+jk)}0 \ \ \ \ \text{by simplifying notation}\\ \end{align}
-
- If $n=i+jk$, where $1<i<j$, then because $n$ is prime, $i$ must be relatively prime to $j$. So let $h$ be a multiplicative inverse of $i$ mod $j$, with $hi=gj+1$. Then \begin{align} S^{1+2j}0\ S^{h(i+jk)}0 &= S^{(1+2j)h(i+jk)}0\ \ \ \ \text{by the mod-1 case}\\ S^{1+2j}0\ S^{i+jk}0 &= S^{(1+2j)(i+jk)}0 \ \ \ \ \text{by cancellation}\\ S^m0\ S^n0 &= S^{mn}0 \end{align}
-
- In any case $S^m0\ S^n0 = S^{mn}0$, QED.