This operation is not mysterious at all! The monoid $(\mathbf N,\circ)$ is isomorphic to a multiplicative submonoid $T$ of the commutative ring $\mathbf Z[\varphi] = \mathbf Z[t]/(t^2-t-1)$, where $\varphi = \tfrac{1+\sqrt{5}}{2}$ is the golden ratio. In particular, it is isomorphic to a submonoid of $\mathbf Z \oplus \mathbf N^{\oplus \mathcal P}$, where $\mathcal P$ is the set of prime ideals in $\mathbf Z[\varphi]$, and the first factor is the powers of the fundamental unit $\varphi$. In fact, the groupification $T^{\text{gp}} \hookrightarrow \mathbf Q(\varphi)^\times \cap \mathbf R_{>0} \cong \mathbf Z \oplus \mathbf Z^{\oplus \mathcal P}$ is an isomorphism; see the corollary below.
I first wrote this as an answer with full details, but a small literature search (starting from Knuth's paper on MathSciNet and looking for forward references) turned up some original sources that are a little slicker than what I wrote. So in the interest of readability, I have removed some of the proofs.
The main reference is the 2-page paper [Arnoux, 1989], explaining the relation to $\mathbf Z[\varphi]$. In addition, [Zhuravlev, 2007] notes that every nonzero element of $\mathbf Z[\varphi]$ can be written (non-uniquely) as $\pm\varphi^{-n}t$ for $t \in T$ and $n \in \mathbf N$, giving the promised computation of $T^{\text{gp}}$ of the corollary below.
Notation. Let $f \colon \mathbf Z[t] \twoheadrightarrow \mathbf Z[\varphi]$ be the quotient map $t \mapsto \varphi$, and let $g \colon \mathbf Z[\varphi] \to \mathbf Z$ be the group homomorphism $a+b\varphi \mapsto b$. Note that $g$ is not a ring homomorphism; in fact $g(\varphi^n) = F_n$ by Binet's formula. Write $h$ for the composition $g \circ f$; this exhibits a polynomial $\sum a_it^i$ as a radix-F expansion of its image $\sum a_i F_i$.
The Fibonacci product $\circ$ has an obvious extension to $\mathbf N \times \mathbf N \to \mathbf N$ by setting $n \circ 0 = n = 0 \circ n$ for all $n \in \mathbf N$ (which does not affect associativity). We view $\mathbf Z[\varphi]$ as a subring of $\mathbf R$ in the obvious way, and it has a conjugation $\overline{(-)} \colon \mathbf Z[\varphi] \to \mathbf Z[\varphi]$ taking $\varphi$ to $1-\varphi = \varphi^{-1} = \tfrac{1-\sqrt{5}}{2}$.
Definition. Define the subset $Z \subseteq t^2\mathbf Z[t]$ of elements of the form $P(t)=\sum_{i=2}^r a_it^i$ with all $a_i \in \{0,1\}$ and $a_ia_{i+1} = 0$ for all $i$. We remove $0$ from this set and add $1$, since this will be our multiplicative unit. Let $T \subseteq \mathbf Z[\varphi]$ be the image of $Z$ under $f$. Then Zeckendorf's theorem shows that $f$ and $g$ give bijections
$$Z \stackrel\sim\to T \stackrel\sim\to \mathbf N.$$
By definition, the map $h$ has the property $h(xy) = h(x)h(y)$ for all $x,y \in Z$. But $Z$ is not closed under multiplication, so this doesn't produce a monoid isomorphism $h \colon Z \stackrel\sim\to \mathbf N$. The key point is:
Proposition [Arnoux, 1989] The set $T$ is closed under multiplication. In particular, $g \colon T \to \mathbf N$ is a monoid isomorphism for the Fibonacci product on $\mathbf N$, i.e.
$$g(xy) = g(x) \circ g(y)\qquad \text{for all } x, y \in T.$$
This also gives a conceptual proof of associativity of $(\mathbf N,\circ)$ (which is not used in the proof). So we see that $\mathbf Z[\varphi]$ interpolates nicely between $\mathbf Z[t]$ and $(\mathbf N,\circ)$.
The proposition can easily be shown by hand using the multiplicative structure of $f$ together with Lemma 3 of Knuth; see the revision history of this post for such a proof. Instead, Arnoux deduces it from the following:
Lemma [Arnoux, 1989]. The set $T \setminus \{1\}$ is given by the elements $t=a+n\varphi$ with $a,n \in \mathbf Z_{>0}$ such that $\overline{\!\ t\ \!} \in (\varphi-2,\varphi-1)$. The map $g^{-1}$ is given by $n \mapsto a_n + n\varphi$, where
$$a_n = \left\lfloor (n+1)\tfrac{-1+\sqrt{5}}{2} \right\rfloor$$
is Hofstadter's $G$-sequence (OEIS A005206).
See Lemmas 2 and 3 in Arnoux. The final statement is not there, but can easily be deduced from the first. This shows that $T$ is closed under multiplication as $(\varphi-2,\varphi-1) \subseteq (-1,1)$. In addition, it gives the clean formula
$$n \circ m = nm + na_m + ma_n.$$
Note that in the first statement, we don't need the assumption $a > 0$. Indeed, if $a \leq 0$ we get $\overline{\!\ t\ \!} = a+n\overline\varphi \leq \overline\varphi < \varphi-2$ since $\overline\varphi < 0$ and $n \geq 1$.
Finally, we need the following observation that appears to be due to [Zhuravlev, 2007]; see Proposition 4.1.
Lemma. Every element $x \in \mathbf Z[\varphi]$ can be written as $\pm\varphi^{-n} \cdot t$ for some $t \in T$ and $n \in \mathbf N$.
Since $T$ and $\varphi$ are positive, the sign agrees with the sign of $x$. The proof (and the whole paper) is notationally heavy (and logically hard to follow), so let me include an argument here.
Proof. If $t = \pm\varphi^n \cdot x$, then $\overline{\!\ t\ \!} \in (\varphi-2,\varphi-1)$ for $n \gg 0$ since $\lvert \overline\varphi \rvert < 1$. If $t = a+n\varphi$, we necessarily have $n \neq 0$, for otherwise $t$ and therefore $x$ is $0$. Wihout loss of generality, we may assume $n > 0$. By the previous lemma (and the discussion after), we get $t \in T$, so $x = \pm\varphi^{-n} \cdot t$ as desired. $\square$
(Notational note: what Zhuravlev calls $\delta(n)$ is related to my $g^{-1}(n)$ via $-\delta(n)/\varphi = \overline{g^{-1}(n)}$. Zhuravlev's Fibonacci sequence is off by $1$ compared to Knuth.)
Corollary. There exists a choice of generators of each prime ideal of $\mathbf Z[\varphi]$ inducing an injection $\psi \colon T \hookrightarrow \mathbf Z \oplus \mathbf N^{\oplus \mathcal P}$. For any such choice, the map $\psi^{\text{gp}} \colon T^{\text{gp}} \to \mathbf Z \oplus \mathbf Z^{\oplus \mathcal P}$ is an isomorphism, and the quotient map $T \to \mathbf N^{\oplus \mathcal P}$ is surjective.
Proof. We saw that $T$ is a submonoid of $(\mathbf Z[\varphi]\setminus\{0\},\times)$. Since $\mathbf Z[\varphi]$ is a real quadratic principal ideal domain, we produce an isomorphism $(\mathbf Z[\varphi]\setminus\{0\},\times) \cong \mathbf Z/2 \oplus \mathbf Z \oplus \mathbf N^{\oplus \mathcal P}$ by picking generators for each prime ideal. By the lemma, we may pick the generator for each prime ideal to be in $T$. Since all elements of $T$ and all chosen representatives are positive, we don't need the sign factor $\mathbf Z/2$, giving an embedding $T \hookrightarrow \mathbf Z \oplus \mathbf N^{\oplus \mathcal P}$. The final two statements follow immediately from the lemma since $\varphi^n \in T$ for all $n \geq 2$. $\square$
While subgroups of free abelian groups are free, the same does not hold for commutative monoids, i.e. there is no unique factorisation into irreducible elements. For instance, the elements $\varphi^n$ for $n \geq 2$ and $n=0$ are in $T$, but $\varphi$ is not (as $g(\varphi) = 1 = g(\varphi^2)$ and $g$ is injective on $T$). So $\varphi^6$ factors both as $(\varphi^2)^3$ or $(\varphi^3)^2$. The result above is probably the most precise you are going to get (also because it depends on infinitely many choices).
Note that negative powers of $\varphi$ are not in the image. In fact, $(\mathbf N,\circ)$ is a sharp monoid: if $a,b \neq 0$, then $a \circ b \neq 0$. I don't know if there exists a choice of generators for which no element picks up a negative power of $\varphi$ (i.e. the image is in $\mathbf N \oplus \mathbf N^{\oplus \mathcal P} \subseteq \mathbf Z \oplus \mathbf N^{\oplus \mathcal P}$), but this seems unlikely to me.
References.
[Arnoux, 1989] P. Arnoux, Some remarks about Fibonacci multiplication. Appl. Math. Lett. 2.4, p. 319-320 (1989). DOI:10.1016/0893-9659(89)90078-5
[Zuravlev, 2007] V. G. Zhuravlev, Sums of squares over the Fibonacci $\circ$-ring. Zap. Nauchn. Semin. POMI 337, p. 165-190 (2006). Translation in J. Math. Sci., New York 143.3, p. 3108-3123 (2007). DOI:10.1007/s10958-007-0195-1
