Which polynomials fix the unit circle? Find all polynomials $P(x)$ with real coefficients such that for every $x,y\in \mathbb{R}$ satisfying $x^2+y^2=1$ we have
$$P(x)^2+P(y)^2=1$$
 A: Here is a "simple" solution.

Lemma.
  Consider a polynomial
  $Q \in\mathbb{C}[X]$ 
   such that $|{Q(e^{it})}|=1$ for  $t\in\mathbb{R}$. Then $Q(X)=\lambda X^n$ for some nonnegative integer $n$, and some
   complex $\lambda$ with $|\lambda|=1$.

Proof.
Clearly $Q$ is not zero, so, let $d=\deg Q$, and suppose that
$$
Q(X)=a_0+a_1X+\cdots+a_d X^d.
$$
Also, let  $m=\mathop{\rm val}(Q)=\min\{k\leq d: a_k\ne 0\}$.
Suppose that $m<d$. The coefficient of $e^{i(d-m)t} $ in the expansion of $|Q(e^{it})|^2$ is $a_d\overline{a_m}$, so
$$
a_d\overline{a_m}=\frac{1}{2\pi}\int_0^{2\pi}
|{Q(e^{it})}|^2e^{i(m-d)t}\,dt=
\frac{1}{2\pi}\int_0^{2\pi}
e^{i(m-d)t}\,dt=0
$$
which is absurd since $a_d\ne 0$ and $a_m\ne 0$.
It follows that $d=m$ so we can take $n=d$, $\lambda=a_d$, and the lemma is proved.
$\qquad\square$
$\qquad$ Consider a polynomial $P(X)\in\mathbb{R}[X]$ that satisfies the proposed condition. That is
$$
\forall\, t\in \mathbb{R},\quad \left|P(\cos t)+iP\left(\cos\left(t-\frac{\pi}{2}\right)\right)\right|=1$$
Now, let $d=\deg P$, then 
there are $(b_0,b_1,\ldots,b_d)\in\mathbb{R}^{d+1}$, such that
$$
P(X)=\sum_{k=0}^d b_kT_k(X)
$$
where $T_k$ is the Chebyshev polynomial of the first kind and degree $k$. (because Chebyshev's polynomials of the first kind constitute a basis for $\mathbb{R}[X]$.)
Now, if $q(t)=P(\cos t)+iP\left(\cos\left(t-\frac{\pi}{2}\right)\right)$ then
\begin{align*}
q(t)&=\sum_{k=0}^db_k(\cos(kt)+i\cos(kt-k\pi/2))\\
&=(1+i)b_0+\frac{1}{2}\sum_{k=1}^db_k\left(e^{ikt}+e^{-ikt}\right)
+\frac{i}{2}\sum_{k=1}^db_k\left(i^{-k}e^{ikt}+i^ke^{-ikt}\right)\\
&=(1+i)b_0+\frac{1}{2}\sum_{k=1}^db_k(1+i^{1-k}) e^{ikt}
+\frac{1}{2}\sum_{k=1}^db_k(1+i^{1+k}) e^{-ikt} \\
&=e^{-idt}Q(e^{it})
\end{align*}
where 
$$
Q(X)=\frac{1}{2}\sum_{k=1}^db_k(1+i^{1+k}) X^{d-k}
+(1+i)b_0X^d+\frac{1}{2}\sum_{k=1}^db_k(1+i^{1-k}) X^{d+k}
$$
By assumption we have $|{Q(e^{it})}|^2=1$ and, according to the Lemma, $Q$ must be a monomial. So, we have one of the following cases:


*

*$|{(1+i)b_0}|=1$ and $b_1=\ldots=b_n=0$. This corresponds to $P(X)=\pm\frac{1}{\sqrt{2}}T_0(X)$. (because $b_0$ is real).

*For some $k\in\{1,\ldots,d\}$ we have $|{b_k(1+i^{1+k})}|=2$,  $|{b_k(1+i^{1-k})}|=0$, and the other $b_j$'s are $0$. Since 
$b_k\ne0$ the second condition implies that $k=3\pmod{4}$ and replacing in the first we get $b_k=\pm1$. So in this case $P(X)=\pm T_k(X)$ for some $k=3\pmod{4}$.

*For some $k\in\{1,\ldots,d\}$ we have $|{b_k(1+i^{1-k})}|=2$,  $|{b_k(1+i^{1+k})}|=0$, and the other $b_j$'s are $0$. Again, the second condition implies that $k=1\pmod{4}$ and replacing in the first we get $b_k=\pm1$. So in this case $P(X)=\pm T_k(X)$ for some $k=1\pmod{4}$.
Finally, we have proved that either $P(X)=\pm\frac{1}{\sqrt{2}}$ or $P(X)=\pm T_{k}$ for some odd integer $k$. The converse is trivially true since the converse of the Lemma is trivial. Done.
A: Using some complex analysis it can be done in a shorter way.
Suppose that $P$ is such a nonconstant real polynomial, and consider the complex rational function
$$
Q(z)=P\left(\frac{z+\frac1z}2\right) + iP\left(\frac{z-\frac1z}{2i}\right).
$$
We can see that $Q$ may have a pole at $0$ with order at most $\deg P$ and another pole at $\infty$ and nowhere else.
Since $Q$ maps the unit circle to itself, it must be a Blaschke-polynomial. But all the Blaschke factors must have their poles at $0$ or $\infty$; therefore, $Q$ is a multiple of a power of $z$, with a unit coefficient:
$Q(z)=e^{ia} z^k$
with some integer $k$ and real $a$.
Substituting $z=e^{it}$ we get 
$$
P(\cos t) +i P(\sin t) = \cos(kt+a) +i \sin(kt+a),
$$
so $P(\cos t)= \cos(kt+a)$ and $P(\sin t)=\sin(kt+a)$.
The function $P(\cos t)$ is even, so $a$ is a multiple of $\pi$. Hence,
$P(\cos t)= \pm\cos(kt)$ and $P(\sin t)=\pm\sin(kt)$.
The first equation shows that $P(x)= \pm T_{|k|}$, so $\pm P$ is a Chebyshev-polynomial.
The second equation shows that $k$ must be odd.
A: If I am not wrong, I have a (complicated) solution.
A) Put $x(t)=P(\cos(t))$, $y(t)=P(\sin(t))$. We have that $x,y$ are ${\mathcal C}^{\infty}$, and as $(x(t))^2+(y(t))^2=1$ for all $t$, there exists a fonction $f$, ${\mathcal C}^{\infty}$, such that $x(t)=\cos(f(t))$ and $y(t)=\sin(f(t))$. We get $\sin(t)P^{\prime}(\cos(t))=f^{\prime}(t)\sin(f(t))$ and $\cos(t)P^{\prime}(\sin(t))=f^{\prime}(t)\cos(f(t))$. Using $\sin(f(t))=P(\sin(t))$ and $\cos(f(t))=P(\cos(t))$, and multiplying by $P(\sin(t))$ and $P( \cos(t))$, we get 
$$f^{\prime}(t)=\sin(t)P(\sin(t))P^{\prime}(\cos(t))+\cos(t)P(\cos(t))P^{\prime}(\sin(t))=A(t)$$ 
B) Now we have for $A(t)$ an expression (with finite number of terms):
$$A(t)=a_0+\sum_{k\not =0}(a_k\cos(kt)+b_k\sin(kt))$$
 for some real constants $a_0, a_k, b_k.$ It is clear that $\int_0^{2\pi}A(t)dt=f(2\pi)-f(0)=2\pi a_0$. But we have $\displaystyle P(1)=\cos((f(0))=\cos(f(2\pi))$ and $P(0)=\sin(f(0))=\sin(f(2\pi))$. Thus $f(2\pi)-f(0)\in 2\pi\mathbb{Z}$, and hence $a_0\in \mathbb{Z}$.
We get that $$f(t)=a_0t+c+\sum_{k\not =0}(\alpha_k\cos(kt)+\beta_k\sin(kt))=a_0t+B(t)$$
 for some new constants $c, \alpha_k, \beta_k\in \mathbb{R}$
We have $$P(\cos(t))+iP(\sin(t))=\exp(ia_0t+iB(t))$$ 
C) Now we have $\displaystyle P(\cos(t))+iP(\sin(t)))\exp(-ia_0t)=\exp(iB(t))$.  But as $a_0\in \mathbb{Z}$, $P(\cos(t))+iP(\sin(t)))\exp(-ia_0t)=D(\exp(it),\exp(-it))$ where $D$ is in $\mathbb{C}[x,y]$, and also $iB(t)=E(\exp(it), \exp(-it))$, for $E\in \mathbb{C}[x,y]$. Now the two functions $\displaystyle D(z,1/z)$ and $\displaystyle \exp(E(z,1/z))$ are analytic in $U=\mathbb{C}-\{0\}$, and they are equal on the unit circle. Hence they are equal on $U$.
Write now $\displaystyle D(z,1/z)=\frac{G(z)}{z^M}$ with $M\in \mathbb{Z}$, and $G$ a polynomial in $z$ with $G(0)\not =0$. Suppose that $G$ is not constant. Then there exists $u\in \mathbb{C}$, not $0$, such that $G(u)=0$. But then we get $\exp(E(u, 1/u))=0$, a contradiction. Hence $G$ is a constant $c$.
D) We have proven that $\displaystyle D(z,1/z)=\frac{c}{z^M}$, and replacing $z$ by $\exp(it)$, we have:
$$P(\cos(t))+iP(\sin(t))=c\exp(iNt)$$, with $N\in \mathbb{Z}$. The constant $c$ is clearly of modulus $1$, hence $c=\exp(id)$, $d\in \mathbb{R}$. We have proven that $P(\cos(t))=\cos(Nt+d)$ and $P(\sin(t))=\sin(Nt+d)$, with $N\in \mathbb{Z}$. 
Now the answer by Jean-Claude Arbaut finish the job. 
A: Not (yet) an answer, but perhaps a strategy that could streamline @Jean-Claude's argument.

Parameterize $(x,y)$ by $(\cos\theta,\sin\theta)$, and then invoke complex exponential relations:
$$x \to \cos\theta \to \frac12\left( e^{i\theta}+ e^{-i\theta}\right) \qquad
y \to \sin\theta \to \frac{1}{2i} \left( e^{i\theta}-e^{-i\theta}\right)$$
These allow us to write the target relation as a "polynomial" (with negative exponents) in $e^{i\theta}$:
$$Q(e^{i\theta}) := P(x)^2 + P(y)^2 - 1$$
Note that $Q$ is identically zero for all $\theta$. In particular, taking our polynomial $P$ to have degree $n$, and writing
$$P(z) := \sum_{k=0}^{n} a_n z^n$$
we have that $Q$ vanishes at the "$n$-th roots of unity", where $\theta_k := 2\pi k/n$ and $k=0$, $1$, $\dots$, $n-1$.
The relations $Q(\;\exp(i\theta_k)\;) = 0$ form a system of $n$ non-linear equations in $n+1$ unknowns, $a_k$. That's just underdetermined. However, if we take as given (based on other answers) that $P$ must be an odd polynomial ---with $n = 2m-1$ and $a_{\text{even}} = 0$--- then the reduced system
$$Q(\;\exp(i\theta_{\text{odd}})\;) = 0$$
has $m$ equations in $m$ unknowns, $a_{\text{odd}}$. In theory, this system is solvable; by @Jean-ClaudeArbaut's argument, the solution is unique (up to sign, since both $P$ and $-P$ satisfy the target relation), giving coefficients of Chebyshev polynomials of the first kind. Analysis of the system may give a more-direct proof of this fact, although I haven't had much luck so far.

Example. $n = 3$, so that $P(z) = a_1 z + a_3 z^3$ (with $a_3 \neq 0$).
Define $\omega := \exp(2i\pi/3)$, and our system is
$$\begin{align}
Q(\omega^1) = 0 &\quad\to\quad a_3 \left( 4 a_1 + 3 a_3 \right)\left(1+\omega^2\right) - \omega \left( 16 - 16 a_1^2 - 24 a_1 a_3 - 10 a_3^2 \right) = 0 \\
Q(\omega^3) = Q(1) = 0 &\quad\to\quad (a_1+a_3)^2 = 1
\end{align}$$
The latter equation gives $a_1 = \pm 1 - a_3$; substitution into the the first equation gives
$$( 1 - \omega )^2 \; a_3\;(a_3\mp 4) = 0 \qquad\to\qquad a_3 = \pm 4\quad\text{(since $a_3 \neq 0$)}$$
Consequently,
$$(a_1, a_3)\;\in\;\left\{\;(-3,4),\;(3,-4)\;\right\} \qquad\to\qquad P = \pm T_3$$
where $T_3$ is the Chebyshev polynomial.

With $n=5$, the equations are already too complicated for me to want to write down. Nevertheless, we can use, say, the method of resultants to eliminate the "lesser" $a_k$s, to get a unique(ish) $a_5 = \pm 16$, and so forth.

For all (odd) $n$, our system of equations always contains the simple relation
$$Q(\omega^n) = Q(1) = 0 \quad\to\quad a_1+a_3+\cdots+a_n = \pm 1$$
where $\omega = \exp(2\pi i/n)$ is the corresponding principal root of unity.
Moreover, the remaining equations use identical "coefficients" on the powers of $\omega$, just permuted. For instance, with $n = 5$, if we have
$$Q(\omega^1) = 0 \qquad\to\qquad 
b_0 (1+\omega) + b_1 ( \omega^2 + \omega^4 ) + b_2 \omega^3 = 0$$
for appropriate $b_k$s, then
$$Q(\omega^3) = 0 \qquad\to\qquad
b_0 ( 1 + \omega^3 ) + b_1 ( \omega + \omega^2 ) + b_2 \omega^4 = 0$$ 
And, of course, there's a good deal of structure in the $b_k$s, which arise from linear combinations of binomial expansions of powers of $(\omega+\omega^{-1})$ and $(\omega - \omega^{-1})$.
Marshaling these facts in just the right way could possibly make the Chebyshev connection more immediate.
A: The following solution is for $1=P_n(x)+P_n(y)$, not for $1=(P_n(x))^2+(P_n(y))^2$.
Let $x=\cos z$ and $y=\sin z$ so that $x^2+y^2=\cos^2 z+\sin^2 z=1$ can be trivially satisfied.
Let $P_n(z)=\sum_{k=0}^{n} a_{2k} z^{2k}$.
$$1=P_2(\sin z)+P_2(\cos z)=2a_0+a_2 \implies a_2=1-2a_0\tag{1}$$
$$1=P_6(\sin z)+P_6(\cos z)=(1/8)(-8 + 16a_0 + 8a_2 + 6a_4 + 5a_6 + (2a_4 + 3a_6)\cos(4x))$$
$$\implies a_4 = -3 (-1 + 2 a_0 + a_2), a_6 = 2 (-1 + 2 a_0 + a_2)\tag{2}$$
$$1=P_{10}(\sin z)+P_{10}(\cos z)=(1/128)(-128 + 256a_0 + 
      63a_{10} + 128a_2 + 96a_4 + 80a_6 + 70a_8 + 4(15a_10 + 8a_4 + 12a_6 + 
           14a_8)\cos(4x) + (5a_10 + 2a_8)\cos(8x))$$
$$\implies a_6 = -2 (-5 + 10 a_0 + 5 a_2 + 2 a_4), a_8 = 
 5 (-3 + 6 a_0 + 3 a_2 + a_4), a_{10} = -2 (-3 + 6 a_0 + 3 a_2 + a_4)\tag{3}$$
Thus there seem to exist infinite many such polynomials for $P_2(z)$, $P_6(z)$, $P_{10}(z)$, or $P_{2+4k}(z),k=0,1,2,...$. 
A: The following solution is for $1=(P_n(x))^2+(P_n(y))^2$.
Let $x=\cos z$ and $y=\sin z$ so that $x^2+y^2=\cos^2 z+\sin^2 z=1$ can be trivially satisfied.
Let $t=\arctan(z/2)$ so that $x=\cos z=2t/(1+t^2)$ and $y=\sin z=(1-t^2)/(1+t^2)$.
Let $P_n(z)=\mu \sum_{k=0}^{n} a_{k} z^{k},\mu=\pm 1$.
We first consider $n=3$:
$$\left(P_3(\frac{2t}{1+t^2})\right)^2+\left(P_3(\frac{1-t^2}{1+t^2})\right)^2-1=0\tag{1}$$
We now multiply (1) by $(1+t^2)^6$ and obtain:
$$\sum_{n=0}^{12}b_{n}t^{n}=0 \tag{2}$$
where $b_n$ are functions of bilinear combinations of $a_k,k=0,1,2,3$.
Setting $b_n=0,n=0,1,2,...,12$, we finally obtain nonzero coefficients of the 3 sets of solutions:
(A) $a_0=\frac{1}{\sqrt{2}}$.
(B) $a_1=1$.
(C) $a_1=3,a_3=-4$.
The last sets in (C) is new (not mentioned in the comments section).
EDIT:
We next consider $n=5$.
$$\left(P_5(\frac{2t}{1+t^2})\right)^2+\left(P_5(\frac{1-t^2}{1+t^2})\right)^2-1=0\tag{3}$$
We now multiply (3) by $(1+t^2)^{10}$ and obtain:
$$\sum_{n=0}^{20}b_{n}t^{n}=0 \tag{4}$$
where $b_n$ are functions of bilinear combinations of $a_k,k=0,1,2,...,5$.
Setting $b_n=0,n=0,1,2,...,20$, we finally obtain the nonzero coefficients of the extra  set of solutions:
(D) $a_1=5,a_3=-20,a_5=16$.
EDIT:
Because only solutions in (A) contain $a_{2n},n=0,1,2,...$, all other solutions in (B),(C), (D) contains $a_{2n+1},n=0,1,2,...$, we are lead to define odd polynomials $Q_{2n+1}(x)$ as:
$Q_{2n+1}(z)=\sum_{k=0}^{n} a_{2k+1} z^{2k+1}$.
Using the same method as shown above for $P_3(x),P_5(x)$, we obtain nonzero coefficients of the new solutions in terms of $Q_{2n+1}(z)$:
(E) $a_1=7,a_3=-56,a_5=112,a_7=-64$.
(F) $a_1=9,a_3=-120,a_5=432,a_7=-576,a_9=256$.
