Find all polynomials $P$ such that $P(x^2+1)=P(x)^2+1$

One solution is $P(x) = x^2 + 1$:

$$ P(x^2 + 1) = (x^2 + 1)^2 + 1 = P(x)^2 + 1$$

Another solution is $P(x) = x^4 + 2x^2 + 2$:

$$ \begin{align} P(x^2 + 1) &= (x^2 + 1)^4 + 2(x^2 + 1)^2 + 2 \\&= x^8 + 4 x^6 + 8 x^4 + 8 x^2 + 5 \\&= (x^4 + 2x^2 + 2)^2 + 1 \\&= P(x)^2 + 1 \end{align}$$

How did I find these? Turns out it's obvious after rewriting the equation! Let $Q(x) = x^2+1$. Then the equation is

$$ P(Q(x)) = Q(P(x)) $$

or more succinctly, $P \circ Q = Q \circ P$. It's now clear that $P=Q$ is one solution, $P = Q \circ Q$ is another solution, $P = Q \circ Q \circ Q$ is yet another, and so forth.

Note that $P(x) = x$ is in this family as well, as the identity function is the empty repeated composition (just like the empty sum is $0$ and the empty product is $1$).

This argument is easily adapted to show the set of solutions is a monoid under composition.

I haven't worked out the complete solution though. Alas I don't know much about the monoid of polynomials under composition.

This answer to a similar question cites a 1922 theorem of Ritt that contains a characterization of polynomials that commute under composition; in particular, we can conclude the repeated compositions of $Q$ are indeed the entire solution space, excluding nonconstant polynomials. The remaining solutions are thus the two functions $P(x) = \beta$ where $\beta$ is a solution to $\beta = Q(\beta)$.

A full characterization of solutions which is quite short and sweet.

Plugging in $x = \pm a$, we see that $P(x) = \pm P(-x)$ for each $x$. Clearly one sign holds infinitely often. By looking at the roots of the polynomial $P(x) - P(-x)$ or $P(x) + P(-x)$, we easily see that $P(x)$ is either odd or even. If $P(x)$ is odd, then note that $P(0) = 0$. Then $P(1) = 1$, $P(2) = 2$, $P(5) = 5$, etc. and by induction we can get infinitely many values $a$ such that $P(a) = a$. By looking at the roots of the polynomial $P(x) - x$ we see that $P(x) = x$ for all $x$ then.

Now suppose $P$ is even. This means $P(x) = Q(x^2)$ for some polynomial $Q$. Let $R(x) = x^2 + 1$. Remark there clearly exists a polynomial $S(x)$ such that $P(x) = S(x^2+1)$, just shift $Q$. Now remark that $P \circ R = R \circ P$. Hence $S \circ R \circ R = R \circ S \circ R \implies S \circ R = R \circ S$, so $S$ is a solution to the functional equation. By an easy induction on degree, it follows the only even solutions to the equations are iterations of $R$, so we are done. (EDIT : Also there is the constant solution $P(x) = c$ where $c = c^2 + 1$, because the induction starts on degree 1 which I forgot)

There is another method (besides that posted here) to prove that solutions have the structure $\sum_{k=0}^n\,a_{2k}x^{2k}$ without constructing formulas to

calculate the coefficients as in my other answer.

We have $$P(x)^2-P(-x)^2=(P(x^2+1)-1)-(P((-x)^2+1)-1)=0$$

At least one of $P(x)-P(-x)$ and $P(x)+P(-x)$ must have infinitely many zeros and therefore must be identical to $0$.

Asume that $P(-x)=-P(x)$. This means that $P$ is an odd function and $P(0)=0$. Let us define $$Q(x)=x^2+1$$ then $P(Q^n(0)=Q^n(0)$ and so $$ P(x)-x=0, \quad x=0,\,Q(0),\,Q^2(0),\ldots$$ for infintely many $x$, so $P(x)=x$.

$P$ must be an even polynomial if not $P(x)=x$. But the even polynomials are exactly the polynomial that contain only even powers of x.

The polynomials $Q^n$ are even polynomials that satisfy the functional equation. They have a degree of $2^n$

In contrast to my other answer I was not able to prove with this method that there is at most one polynomial of a certain degree that satisfy the functional equation.

Proof
a polynomial is even $\Longleftrightarrow$ the polynomial contains only even powers of $x$

The even polynomial $P(x)$ can be expressed as sum $f(x)+g(x)$ were $f(x)$ is a polynomial that only contains even powers of $x$ and $g(x)$ is a polynomial that contains only odd powers of $x$. We have $$P(x)-f(x)=g(x)$$ The left side is an even polynomial, the right side an odd one, so $P(x)=f(x)$. So polynomials that contain only even powers of $x$ are exactly the even polynomials.

EDIT:

The following completes the proofs:

Lemma $Q(x)=x^2+1$,$P(Q(x)=Q(P(x))$ and $P(x)$ is even. Then there is a polynomial $T(x)$ with $P=T \circ Q = Q \circ T$

So a polynomial $P$ of degree $n$ can be reduced to a polynomial with degree $\frac{n}{2}$ that also satisfies the functional equation. This process can repeated until we

arrive at a polynomial with odd degree. But the only polynomial with odd degree that satisfies the functional equation is $T(x)=x$.

Therefore the functional equation $$Q \circ P = P \circ Q $$ has exactly the following polynomial solutions:

$$ -\frac{\sqrt{3}i-1}{2} \\ \frac{\sqrt{3}i+1}{2} \\ x \\ Q(x)\\ Q^2(x) \\ Q^3(x) \\ \ldots $$

Proof of the Lemma
We define $R(x)=\sqrt[+]{x-1}$, then $$R(Q(x))=Q(R(x))=x,\quad \forall x>1$$ and $$R(x)^2=x-1, \quad \forall x>1$$ We have $$ Q \circ P = P \circ Q$$ and therefore $$ Q \circ P \circ R = P \circ Q \circ R = P = P \circ R \circ Q , \quad \forall x>1$$ For $x> 1$ the polynomial $$T(x)=P(R(x))=\sum_{k=0}^{n}a_{2n}(x-1)^{n}$$ satisfies $$P=T \circ Q = Q \circ T$$ But if a polynomial equation is satisfied for infinite many $x$ it is satisfied for all $x$.

The coefficients of the polynomials can be efficiently calculated by the formulae $(3)$ and $(4)$ of my other post

Let $P(y)=\sum_{0\le r\le n}a_ry^r$

So, $$P(1+x^2)=\sum_{0\le r\le n}a_r(1+x^2)^r=a_0+a_1(1+x^2)+a_2(1+\binom 21x^2+x^4)+\cdots +a_{n-1}(1+\binom {n-1}1x^2+\binom {n-1}2x^4+\cdots+\binom {n-1}{n-2}x^{2(n-2)}+\binom {n-1}{n-1}x^{2(n-1)}) +a_n(1+\binom n1x^2+\binom n2x^4+\cdots+\binom n{n-1}x^{2(n-1)}+\binom n nx^{2n})$$

$$=x^{2n}a_n+x^{2n-2}(a_n\binom n{n-1}+a_{n-1})+x^{2n-4}(a_n\binom n{n-2}+a_{n-1}\binom {n-1}{n-2}+a_{n-2})+x^2(a_n\binom n1+a_{n-1}\binom{n-1}1+\cdots+a_2\binom21+a_1)+\sum_{0\le r\le n}a_r$$

and $$\{P(x)\}^2+1=\{\sum_{0\le r\le n}a_rx^r\}^2+1$$ $$=a_n^2x^{2n}+x^{2n-1}2a_na_{n-1} +x^{2n-2}(a_{n-1}^2+2a_na_{n-2})+x^{2n-3}2(a_na_{n-3}+a_{n-1}a_{n-2}) +x^{2n-4}(a_{n-2}^2+2a_na_{n-4}+2a_{n-1}a_{n-3})+\cdots+x^2(a_1^2+2a_0a_2)+\sum_{0\le r\le n}a_r^2+1$$

Comparing the coefficients of the different powers of $x$

$r=n\implies a_n=a_n^2\implies a_n=1$ as $a_n\ne0$

$r=n-1\implies 2a_na_{n-1}=0\implies a_{n-1}=0$

$r=n-2\implies a_n\binom n{n-1}+a_{n-1}=a_{n-1}^2+2a_na_{n-2}\implies a_{n-2}=\frac n2$

$r=n-3\implies 2(a_na_{n-3}+a_{n-1}a_{n-2})=0\implies a_{n-3}=0$

$r=n-4\implies a_n\binom n{n-2}+a_{n-1}\binom {n-1}{n-2}+a_{n-2}=a_{n-2}^2+2a_na_{n-4}+2a_{n-1}a_{n-3}\implies a_{n-4}=\frac {n^2}8=\frac1{2!}\left(\frac n2\right)^2$

$r=n-5\implies 2(a_na_{n-5}+a_{n-1}a_{n-4}+a_{n-2}a_{n-3})=0\implies a_{n-5}=0$

$r=n-6\implies 2(a_na_{n-6}+a_{n-1}a_{n-5}+a_{n-2}a_{n-4})+a_{n-3}^2=a_n\binom n{n-3}+a_{n-1}\binom{n-1}{n-3}+a_{n-2}\binom{n-2}{n-3}+a_{n-3}\implies a_{n-6}=\frac1{3!}\left(\frac n2\right)^3-\frac n3$