Functional equation $4f(x^2+y^2)=(f(x)+f(y))^2$
Solution 1:
A beginning.
As you said, $$4f_0=4f_0^2\quad\wedge\quad 4f_1=(f_0+f_1)^2$$ implies $$(f_0,f_1)\in\bigl\{(0,0),(1,1),(0,4)\bigr\}$$ and suggests the three solutions $$\phi_0(x):\equiv0,\quad \phi_1(x):\equiv1,\quad\phi_2(x):=4x^2\ .$$ The conjecture is that the $\phi_i$ are the only solutions. One easily verifies the following
Lemma. If $f:\>{\mathbb N}_{\geq0}\to{\mathbb N}_{\geq0}$ is a solution, and $f(z_1)=\phi_i(z_1)$, $f(z_2)=\phi_i(z_2)$ for any two different elements of a triple $T:=\{x,y,x^2+y^2\}$ then one also has $f(z_3)=\phi_i(z_3)$ for the third element of $T$.
It follows that each of the solutions $\phi_i$ propagates through the set of triples as far as it is "connected" to the basic triples $\{0,0,0\}$, $\{0,1,1\}$. How does this connection come about? It comes from pairs $(x,y)$, $(u,v)$ tied to each other by $x^2+y^2=u^2+v^2$. If we have such an incidence and already know that $$f(x)=\phi_i(x),\quad f(y)=\phi_i(y),\quad f(u)=\phi_i(u)$$ then the Lemma allows to conclude that $f(v)=\phi_i(v)$ as well. Therefore it remains to prove that the set of such incidences is rich enough to make the "complex" of good triples $T$ universal.
I did some computer search and found that there is no number $x\leq50$ whose $f$-value is not enforced by $(f_0,f_1)$. To this end we need a list of all incidences with square sum $\leq5000$. E.g., $$2125=3^2+46^2=10^2+45^2=19^2+42^2=30^2+35^2\ .$$