This question is inspired by this thread. However, in this question, I take an arbitrary field instead of $\mathbb{R}$ and drop the assumption that $P(0)$ must be $0$.

Let $K$ be a field. Determine all $P(X)\in K[X]$ such that $$P\big(X^2+1\big)=\big(P(X)\big)^2+1\,.$$ In other words, if $Q(X):=X^2+1$, we are to determine all $P(X)\in K[X]$ that commute with $Q(X)$ with respect to polynomial composition: $$P\big(Q(X)\big)=Q\big(P(X)\big)\,.$$

First suppose that the characteristic of $K$ is not equal to $2$. I have observed that some solutions are given by

(i) $P(X)=\frac{1+\sqrt{-3}}{2}$ and $P(X)=\frac{1-\sqrt{-3}}{2}$ if $\sqrt{-3}\in K$, and

(ii) $P(X)=Q_n(X)$ for $n=0,1,2,3,\ldots$,

where $Q_0(X):=X$ and $Q_n(X):=Q\big(Q_{n-1}(X)\big)$ for all $n=1,2,3,\ldots$. Are these all of the solutions?

If the characteristic of $K$ is equal to $2$, then the given functional equation is equivalent to $$P\left((X+1)^2\right)=\left(P(X)+1\right)^2\,.$$ If $K=\mathbb{F}_2$, then it follows that $$P(X+1)=P(X)+1\,,$$ whence $P(X)$ must take the form $$P(X)=X^{2^{r_1}}+X^{2^{r_2}}+\ldots+X^{2^{r_k}}\text{ or }P(X)=1+X^{2^{r_1}}+X^{2^{r_2}}+\ldots+X^{2^{r_k}}$$ where $k$ is an odd positive integer and $r_1,r_2,\ldots,r_k$ are pairwise distinct nonnegative integers.

Can anybody prove or disprove whether my list is complete when the characteristic of $K$ is not equal to $2$? Bill Dubuque's answer in this thread shows that, if $K$ is a subfield of $\mathbb{C}$, then the list is indeed complete. This should imply that, if $K$ is of characteristic $0$ and the cardinality of $K$ is that of the continuum, then the list is complete (as $K$ can be embedded into $\mathbb{C}$).

Can someone try to make a characterization of the solutions when $K$ is an arbitrary field of characteristic $2$? This case seems to have weird polynomial solutions. For example, when $K=\mathbb{F}_4=\mathbb{F}_2[\alpha]$ with $\alpha^2+\alpha+1=0$, we have the following solutions $P(X)=\alpha$, $P(X)=X$, $P(X)=X+1$, $P(X)=X^2$, $P(X)=X^2+1$, $P(X)=X^2+X+\alpha$, and $P(X)=X^3+\alpha\,X^2+\alpha\,X$.

*UPDATES:*

(1) By balancing the coefficients, it can be shown that, in characteristic not equal to $2$, there is at most one solution of degree $d>0$, and this solution lies within $\mathbb{F}_p[X]$, where $\mathbb{F}_p$ is the prime field of $K$ (with $\mathbb{F}_0:=\mathbb{Q}$). Hence, it suffices to show that the list is complete when $K$ is a prime field. Ergo, for characteristic $0$, we are done.

(2) In characteristic $2$, via balancing the coefficients, we see that all equations involved are quadratic in the coefficients, whence the solutions must lie within $\mathbb{F}_4[X]$. Consequently, we may take $K$ to be $\mathbb{F}_2$ or $\mathbb{F}_4$.

(3) In odd characteristic $p$, the list is indeed incomplete. Note that $P(X)=X^{p^k}$ is a solution for all $k=0,1,2,\ldots$. Let $R_p(X):=X^p$. Is it true that all solutions must take the form $$P(X)=\left(Q^{\circ k}\circ R_p^{\circ l}\right)(X)$$
for some $k,l\in\mathbb{Z}_{\geq 0}$? Here, $\left(F^{\circ 0}\right)(X):=X$ and $$\left(F^{\circ l}\right)(X):=\left(\underbrace{F\circ F\circ \ldots \circ F}_{F\text{ occurs }l\text{ times}}\right)(X)$$ for each $F(X)\in K[X]$ and $l\in\mathbb{Z}_{> 0}$.

(4) It turns out that the characterization in (3) is not complete either, at least not in characteristic $3$. The polynomials
$$X^{5}+X^3-X\,,\;\;X^7-X^5-X^3-X\,,\text{ and }X^{11}+X^9-X^7+X^5+X^3+X$$
are solutions in characteristic $3$ which do not arise from the formula in (3). The case of characteristic $3$ may be the only special case where (3) does not give a full characterization.

(5) Combined with Hurkyl's answer, we know all solutions in characteristics $0$, $2$, and $3$. So far, I have not been able to find any solutions outside (3) in characteristics greater than $3$.