Is there an orthonormal set of polynomials whose derivatives are orthogonal? Let $T:\mathbb{R}[x]\rightarrow \mathbb{R}[x]$ be the linear map $p(x)\mapsto \frac{d}{dx} p(x)$. Is there an orthonormal basis $b_1,b_2,...$ for the polynomial ring $\mathbb{R}[x]$ with respect to the $L^2([0,1])$ norm such that the elements $T(b_1), T(b_2),...$ are orthogonal?
Edit: Perhaps I should add some motivation. It often happens in linear algebra that methods in the finite dimensional case carry over to the infinite dimensional case. For instance, Fourier series are a consequence of finding an orthonormal basis of $L^2([-\pi,\pi])$. The SVD is an example that as far as I know doesn't have that many applications in pure math, despite it being extremely useful in applied math. If you know how the SVD works, you'll recognize that this is exactly what I'm asking for here, in the simplest example of an interesting linear operator on $\mathbb{R}[x]$ that I know, $\frac{d}{dx}$.
 A: The system $\{p_n\}_{n=1}^\infty$ contains a non-constant polynomial. We may assume that $p_1$ is non-constant.
Then
$$p_1(x)=a_nx^n+\ldots +a_1x+a_0,\ n\ge 1, a_n\neq  0$$
Consider the linear functional $\varphi$ acting on polynomials according to the formula
$$\varphi(q)=\int\limits_0^1 p_1'(x)q'(x)\,dx$$
Assume $\varphi(q)=0.$ We have
$$q(x)=a_1p_1(x)+a_2p_{n_2}(x)+\ldots+a_kp_{n_k}(x), \ 2\le n_1<\ldots<n_k$$ and
$$q'(x)=a_1p_1'(x)+a_2p_{n_2}'(x)+\ldots+a_kp_{n_k}'(x)$$
By orthogonality we get
$$0=\varphi(q)=a_1\int\limits_0^1[p_1'(x)]^2\,dx$$
Therefore $a_1=0.$
This implies that $\ker\varphi\subset {\rm span}\{p_n\,:\, n\ge 2\}.$
By orthogonality we get $p_n\in \ker\varphi$ for $n\ge 2.$ Thus
$$\ker\varphi={\rm span}\{p_n\,:\, n \ge 2\}$$
The equality means that
$$\varphi(q)=\int\limits_0^1 p_1'(x)q'(x)\,dx=0 \iff \psi(q):=\int\limits_0^1 p_1(x)q(x)\,dx=0$$
The nonzero functionals $\varphi$ and $\psi$ are linearly dependent as their kernels are equal.
Assume $\varphi=\lambda \psi,$ $\lambda\neq 0.$
Applying integration by parts gives
$$\varphi(q)=q(1)p_1'(1)-q(0)p_1'(0)-\int\limits_0^1 p_1''(x)q(x)\,dx$$ Hence
$$q(1)p_1'(1)-q(0)p_1'(0)-\int\limits_0^1 p_1''(x)q(x)\,dx=\lambda \int\limits_0^1 p_1(x)q(x)\,dx$$
The equality extends from the space of polynomials to $ C[0,1].$
Therefore we obtain the equality of signed measures
$$p_1'(1)\delta_1-p_1'(0)\delta_0-p_1''(x)\,dx=\lambda p_1(x)\,dx$$
This implies  $p_1''=-\lambda p_1.$ As $\lambda\neq 0$ we get  $p_1=0$ which gives a contradiction.
A: The simplest example is provided by the Chebyshev polynomials of the first and the second kind. They are orthogonal on $[-1,1]$ with respect to the weights $(1-x^2)^{-1/2}$ and
$(1-x^2)^{1/2}$ respectively. So they do not meet the requirement of being orthogonal with respect to the Lebesgue measure.
The Chebyshev polynomials of the first kind are defined by $$T_n(x)=\cos(n\arccos x)$$ while those of the second kind by
$$U_n(x)={\sin[(n+1)\arccos x]\over \sin(\arccos x)}$$
Orthogonality follows from orthogonality of trigonometric functions by applying the substitution $t=\arccos x.$
We have
$${d\over dx}\,T_{n+1}(x)=(n+1)U_n(x)$$
By affine change of variables $x=2u-1$ we can move the example to the interval $[0,1].$
A: Maybe this could be useful:
https://digitalcommons.unl.edu/cgi/viewcontent.cgi?referer=&httpsredir=1&article=1032&context=mathfacpub
in particular the Theorem (page 880) that states the follows:
If $\{ \phi_n \} $ and $ \{ \phi_n' \} $ are orthogonal systems of polynomials, then $ \{ \phi_n \} $ may be reduced to the classical polynomials of Jacobi, Laguerre, or Hermite by means of linear transformation on $x$.
A: Looking at the references in the paper mentioned by 3mOo, we find  a paper by Krall with the following result:
theorem . If $\{\phi_n(x)\}$ is a set of orthogonal polynomials with
the weight function $p(x)$ in the finite interval $(a, b)$, and if we assume that the derivatives $\{\phi_n^\prime(x)\}$ also form a set of orthogonal
polynomials in a certain interval $(c, d)$ (infinite or not), with a
non-negative weight function $q(x)$, then $\{\phi_n(x)\}$ is a set of Jacobi
polynomials.
The Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$ are orthogonal on $[-1,1]$ with weight function $(1-x)^\alpha(1+x)^\beta$. So if we want a constant weight, we set $\alpha=\beta=0$.
$P_1^{(0,0)}(x) = x$
$P_3^{(0,0)}(x) = \frac{5}{2}x^3-\frac{3}{2}x$
$\frac{d}{dx}P_1^{(0,0)}(x) = 1$
$\frac{d}{dx}P_3^{(0,0)}(x) = \frac{15}{2}x^2-\frac{3}{2}$
In this case, we see that the derivatives are not orthogonal with the same weight on the same interval.
On the other hand, we have
$\frac{d}{dx} P_n^{(0,0)}(x) = \frac{n+1}{2}P_{n-1}^{(1,1)}$. Hence the derivatives are orthogonal on $[-1,1]$ with weight $(1-x)(1+x)$.
A: As a few people pointed out, if we assume that the orthonormal basis $p_0, p_1, \dots$ is graded, i.e. $\mathrm{deg}(p_n)=n$, then the result follows by an old paper of Hahn (1935). I argue that's overkill, since the following elementary argument works:

First observe that $p_1',p_2',\dots$ have degrees
$0,1,\dots$, hence form an orthogonal basis $B$. The trick: consider
the polynomial $q(x) = x-x^2$. Then $\int_0^1q(x)p_n'(x)dx = [q(x)p_n(x)]_0^1-\int_0^1q'(x)p_n(x)dx$ by the chain rule. The choice
of $q$ makes $[q(x)p_n(x)]_0^1=0$ and $q'(x)$ is a linear polynomial,
hence orthogonal to each $p_n$ for $n\ge 2$, so the second term vanishes too. It follows that
$\int_0^1q(x)p_n'(x)dx=0$ for $n\ge 2$. Since $q$ is nonzero, it must
be a multiple of the remaining vector $p_1'$ in the basis $B$. But
$p_1'$ is constant while $q$ isn't, a contradiction.

I'm happy with this argument but unfortunately it doesn't fully answer my original question: there could be a surprising basis of polynomials of degrees all over the place satisfying the conditions.
Intuitively, I think it should be possible to construct an infinite family (not necessarily spanning $\mathbb{R}[x]$) of orthogonal polynomials with orthogonal derivatives, due to the following heuristic "process". Start with a list $p_0,...,p_k$ satisfying the conditions. Then add a polynomial $p_{k+1}$ of degree $d$ and unknown coefficients. The orthogonality conditions give $2(k+1)$ homogeneous linear equations in the $d+1$ unknown coefficients, so if $d$ is large enough, there should typically be nonzero solutions. Pick one and repeat the process.
In the above argument, one could conceivably bake in the property that $x^n$ lies in the span of the first $N=N(n)$ polynomials, perhaps by adding a condition of certain minors vanishing. To me, this seems unlikely work.
I don't think anyone came close to answering the question, so I will not award the bounty for now.
A: Ok so let's assume this orthonormal basis exists call these polynomials $p_n$. By hypothesis we have:
$$\langle p_n,p_m\rangle= \delta_{nm}$$
and also:
$$\langle p'_n,p'_m\rangle=\delta_{nm}$$
Since both of these are basis of $\mathbb R[x]$ we are in escence asseting that the derivative map is a change of basis, but since the derivative has a kernel this is certainly not possible. For example we know that there must exist some coefficients $a_n$ such that:
$$1=\sum_n a_n p_n$$
since the $p_n$ form a basis, but by taking a derivative we also have:
$$0=\sum_n a_n p'_n$$
so the second set could not also be a basis.
