For polynomials in the monomial basis like $p_n(x) = \sum_{k=0}^N a_k x^k $, the product of 2 polynomials is can be either found though the convolution of the 2 corresponding polynomial vectors or with FFT/IFFT.
I wonder, if there exists a "numerical recipe" to compute the product of 2 polynomials like $ p_n(x) = \sum_{k=0}^N c_k T_k(x)$ (i.e. represented in the Chebyshev basis).
Yes, they boil down to sums of things that look just like convolution or correlation of the coefficients. See equations 2.6, 2.7, and 2.8 here for the exact formulas:
http://www2.mathematik.hu-berlin.de/~gaggle/S09/AUTODIFF/projects/papers/baszenski_fast_polynomial_multiplication_and_convolutions_related_to_the_discrete_cosine_transform.pdf
Furthermore, the sums can be efficiently calculated using FFT-related techniques, which also seems to be outlined in that paper, though I didn't look at the details.
Edit: To make this a complete answer that does not require visiting the reference, here are the relevant formulas:
For $n\in\mathbb{N}$, $p_n$ and $q_n$ are $n^\text{th}$ order polynomials written as Chebyshev series in the $T_k$ polynomials: $$ p_n = \frac{a_0}{2}+\sum_{k=1}^n a_k T_k,\;\;\;q_n = \frac{b_0}{2}+\sum_{k=1}^n b_k T_k $$ The coefficients $a_k$ and $b_k$ are real. Then the product of the two polynomials has order $2n$ and is denoted $r_{2n}$. The Chebyshev series form of that product is $$ r_{2n} = \frac{c_0}{2}+\sum_{k=1}^{2n} c_k T_k $$ And the $c_k$ coefficients can be calculated from the original $a_k$ and $b_k$ coefficients as: $$ 2c_k = \left\{ \begin{array}{lcl} \displaystyle a_0b_0 + 2\sum_{l=1}^n a_l b_l, & {} & k=0 \\ \displaystyle \sum_{l=0}^k a_{k-l}b_l + \sum_{l=1}^{n-k}\left(a_l b_{l+k}+a_{l+k}b_l\right), & {} & k=1,\ldots,n-1 \\ \displaystyle \sum_{l=n-k}^n a_{k-l} b_l, & {} & k=n,\ldots,2n \end{array} \right. $$