[proofplan]
We pull the function back from the Bernstein ellipse to the annulus $\rho^{-1}<|w|<\rho$ by the Joukowski map $w \mapsto (w+w^{-1})/2$. The pullback has a Laurent expansion, and the symmetry $w \mapsto w^{-1}$ identifies the Chebyshev coefficients with twice the positive Laurent coefficients. Cauchy's coefficient estimate on the circle $|w|=\rho$ gives the exponential coefficient decay, and the uniform error estimate follows by summing the geometric tail and using $|T_n(x)|\le 1$ on $[-1,1]$.
[/proofplan]
[step:Pull back the problem to an annulus]
Define the open annulus
\begin{align*}
A_\rho := \{w \in \mathbb{C} : \rho^{-1}<|w|<\rho\}.
\end{align*}
Define the Joukowski map $J: A_\rho \to \Omega_\rho$ by
\begin{align*}
J(w):=\frac{1}{2}(w+w^{-1}).
\end{align*}
By the definition of $\Omega_\rho$, this map has image $\Omega_\rho$. Define the pulled-back function
\begin{align*}
F: A_\rho \to \mathbb{C}, \qquad F(w):=f(J(w)).
\end{align*}
Since $J$ is holomorphic on $A_\rho$ and $f$ is holomorphic on $\Omega_\rho$, the composition $F$ is holomorphic on $A_\rho$. Since $f$ extends continuously to $\overline{\Omega}_\rho$, the function $F$ extends continuously to the closed annulus
\begin{align*}
\overline{A}_\rho := \{w \in \mathbb{C} : \rho^{-1}\le |w|\le \rho\}.
\end{align*}
Furthermore, the restriction $J|_{\overline{A}_\rho}: \overline{A}_\rho \to \mathbb{C}$ is continuous, and $J(A_\rho)=\Omega_\rho$ by the definition of the Bernstein ellipse. Since every point of $\overline{A}_\rho$ is a limit of points of $A_\rho$, continuity gives $J(\overline{A}_\rho)\subset \overline{J(A_\rho)}=\overline{\Omega}_\rho$. Therefore the hypothesis on $f$ gives
\begin{align*}
|F(w)| \le M
\end{align*}
for every $w \in \overline{A}_\rho$.
Finally, $J(w)=J(w^{-1})$ for every $w \in A_\rho$, so
\begin{align*}
F(w)=F(w^{-1})
\end{align*}
for every $w \in A_\rho$.
[/step]
[step:Identify Chebyshev coefficients with symmetric Laurent coefficients]
By the Laurent expansion theorem for holomorphic functions on an annulus, $F$ has a Laurent expansion
\begin{align*}
F(w)=\sum_{k\in\mathbb{Z}} c_k w^k
\end{align*}
converging locally uniformly on $A_\rho$, where each $c_k \in \mathbb{C}$ is the $k$th Laurent coefficient. By the Laurent coefficient formula applied on the unit circle, which is contained in $A_\rho$ because $\rho>1$, for every $k \in \mathbb{Z}$,
\begin{align*}
c_k=\frac{1}{2\pi i}\oint_{|w|=1} F(w)w^{-k-1}\,dw.
\end{align*}
Parametrize the positively oriented unit circle by
\begin{align*}
\gamma: [0,2\pi]\to \mathbb{C}, \qquad \gamma(\theta):=e^{i\theta}.
\end{align*}
Then $d\gamma_\theta=i e^{i\theta}\, d\mathcal{L}^1(\theta)$, and $J(e^{i\theta})=\cos\theta$. Hence, for $n\ge 1$,
\begin{align*}
c_n=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)e^{-in\theta}\,d\mathcal{L}^1(\theta).
\end{align*}
Similarly,
\begin{align*}
c_{-n}=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)e^{in\theta}\,d\mathcal{L}^1(\theta).
\end{align*}
For $k=0$, the same parametrized Laurent coefficient formula gives
\begin{align*}
c_0=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)\,d\mathcal{L}^1(\theta)=a_0,
\end{align*}
where the final equality is the zeroth Chebyshev coefficient convention used in the expansion. Adding the two identities for $c_n$ and $c_{-n}$ and using $e^{in\theta}+e^{-in\theta}=2\cos(n\theta)$ gives
\begin{align*}
c_n+c_{-n}=\frac{1}{\pi}\int_0^{2\pi} f(\cos\theta)\cos(n\theta)\,d\mathcal{L}^1(\theta)=a_n.
\end{align*}
The symmetry $F(w)=F(w^{-1})$ forces $c_n=c_{-n}$ for every $n\in\mathbb{Z}$. Indeed, the Laurent expansion of $F(w^{-1})$ is
\begin{align*}
F(w^{-1})=\sum_{k\in\mathbb{Z}} c_k w^{-k}=\sum_{k\in\mathbb{Z}} c_{-k}w^k,
\end{align*}
and the uniqueness theorem for Laurent coefficients on an annulus gives $c_k=c_{-k}$. Therefore, for every $n\ge 1$,
\begin{align*}
a_n=2c_n.
\end{align*}
[guided]
The purpose of the annulus is that Chebyshev polynomials become ordinary powers after the substitution $x=\frac{1}{2}(w+w^{-1})$. Since $F$ is holomorphic on $A_\rho$, it has a Laurent expansion
\begin{align*}
F(w)=\sum_{k\in\mathbb{Z}} c_k w^k,
\end{align*}
with locally [uniform convergence](/page/Uniform%20Convergence) on $A_\rho$. The Laurent coefficient formula gives
\begin{align*}
c_k=\frac{1}{2\pi i}\oint_{|w|=1} F(w)w^{-k-1}\,dw
\end{align*}
for every $k\in\mathbb{Z}$. The circle $|w|=1$ lies inside $A_\rho$ because $\rho>1$.
Now parametrize the unit circle by the map
\begin{align*}
\gamma: [0,2\pi]\to\mathbb{C}, \qquad \gamma(\theta):=e^{i\theta}.
\end{align*}
This parametrization is positively oriented, and its complex differential is $d\gamma_\theta=i e^{i\theta}\,d\mathcal{L}^1(\theta)$. Since
\begin{align*}
J(e^{i\theta})=\frac{1}{2}(e^{i\theta}+e^{-i\theta})=\cos\theta,
\end{align*}
we obtain, for $n\ge 1$,
\begin{align*}
c_n=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)e^{-in\theta}\,d\mathcal{L}^1(\theta).
\end{align*}
The same calculation with $-n$ in place of $n$ gives
\begin{align*}
c_{-n}=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)e^{in\theta}\,d\mathcal{L}^1(\theta).
\end{align*}
The zeroth coefficient is identified by taking $k=0$ in the same formula:
\begin{align*}
c_0=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)\,d\mathcal{L}^1(\theta)=a_0.
\end{align*}
This is the normalization of the zeroth Chebyshev coefficient used by the truncated expansion in the theorem statement. Adding the two formulas for $c_n$ and $c_{-n}$ converts complex exponentials into cosines:
\begin{align*}
c_n+c_{-n}=\frac{1}{2\pi}\int_0^{2\pi} f(\cos\theta)(e^{-in\theta}+e^{in\theta})\,d\mathcal{L}^1(\theta).
\end{align*}
Since $e^{-in\theta}+e^{in\theta}=2\cos(n\theta)$, this becomes
\begin{align*}
c_n+c_{-n}=\frac{1}{\pi}\int_0^{2\pi} f(\cos\theta)\cos(n\theta)\,d\mathcal{L}^1(\theta)=a_n.
\end{align*}
It remains to use the symmetry of the Joukowski map. Because $J(w)=J(w^{-1})$, we have $F(w)=F(w^{-1})$ on $A_\rho$. The Laurent expansion of $F(w^{-1})$ is
\begin{align*}
F(w^{-1})=\sum_{k\in\mathbb{Z}} c_k w^{-k}.
\end{align*}
Reindexing by replacing $k$ with $-k$ gives
\begin{align*}
F(w^{-1})=\sum_{k\in\mathbb{Z}} c_{-k}w^k.
\end{align*}
Since Laurent expansions on an annulus are unique, the coefficient of $w^k$ in $F(w)$ equals the coefficient of $w^k$ in $F(w^{-1})$, so $c_k=c_{-k}$ for every $k\in\mathbb{Z}$. Therefore
\begin{align*}
a_n=c_n+c_{-n}=2c_n
\end{align*}
for every $n\ge 1$.
[/guided]
[/step]
[step:Estimate the Laurent coefficients by circles inside the annulus]
Fix $n\ge 1$ and let $r$ be a real number with $1<r<\rho$. Since the circle $|w|=r$ is contained in the holomorphic domain $A_\rho$, the Laurent coefficient formula applied on this circle gives
\begin{align*}
c_n=\frac{1}{2\pi i}\oint_{|w|=r} F(w)w^{-n-1}\,dw.
\end{align*}
Parametrize this circle by the map
\begin{align*}
\gamma_r: [0,2\pi]\to\mathbb{C}, \qquad \gamma_r(\theta):=r e^{i\theta}.
\end{align*}
Then $d\gamma_{r,\theta}=ir e^{i\theta}\,d\mathcal{L}^1(\theta)$. Because $r<\rho$, every point with $|w|=r$ lies in $\overline{A}_\rho$, so the previously proved bound $|F(w)|\le M$ applies on this circle. Therefore
\begin{align*}
|c_n|\le \frac{1}{2\pi}\int_0^{2\pi} Mr^{-n-1}r\,d\mathcal{L}^1(\theta).
\end{align*}
Evaluating the integral gives
\begin{align*}
|c_n|\le Mr^{-n}.
\end{align*}
This estimate holds for every $r\in(1,\rho)$. Letting $r\uparrow\rho$ and using the continuity of the function $r\mapsto Mr^{-n}$ on $(0,\infty)$ yields
\begin{align*}
|c_n|\le M\rho^{-n}.
\end{align*}
Since $a_n=2c_n$, we conclude that
\begin{align*}
|a_n|\le 2M\rho^{-n}
\end{align*}
for every $n\ge 1$.
[/step]
[step:Recover the Chebyshev expansion on the interval]
For every $\theta\in\mathbb{R}$, the locally uniform Laurent expansion on the unit circle gives
\begin{align*}
f(\cos\theta)=F(e^{i\theta})=\sum_{k\in\mathbb{Z}} c_k e^{ik\theta}.
\end{align*}
More precisely, [Laurent series](/page/Laurent%20Series) converge normally on compact subannuli of their annulus of holomorphy. Since the unit circle is a compact subset of $A_\rho$, the series $\sum_{k\in\mathbb{Z}} c_k e^{ik\theta}$ converges absolutely and uniformly for $\theta\in\mathbb{R}$, so regrouping the positive and negative modes is justified. Using $c_{-k}=c_k$, $a_k=2c_k$ for $k\ge 1$, and $c_0=a_0$, this becomes
\begin{align*}
f(\cos\theta)=a_0+\sum_{k=1}^{\infty} c_k(e^{ik\theta}+e^{-ik\theta}).
\end{align*}
Because $e^{ik\theta}+e^{-ik\theta}=2\cos(k\theta)$ and $T_k(\cos\theta)=\cos(k\theta)$, we obtain
\begin{align*}
f(\cos\theta)=a_0+\sum_{k=1}^{\infty} a_kT_k(\cos\theta).
\end{align*}
Thus, for every $x\in[-1,1]$, choosing $\theta\in[0,\pi]$ with $x=\cos\theta$ gives
\begin{align*}
f(x)=a_0+\sum_{k=1}^{\infty} a_kT_k(x).
\end{align*}
[/step]
[step:Sum the Chebyshev tail uniformly on $[-1,1]$]
Let $N\in\mathbb{N}_0$. For $x\in[-1,1]$, choose $\theta\in[0,\pi]$ with $x=\cos\theta$. Then
\begin{align*}
|T_n(x)|=|\cos(n\theta)|\le 1
\end{align*}
for every $n\in\mathbb{N}_0$. Therefore,
\begin{align*}
|f(x)-S_Nf(x)|\le \sum_{n=N+1}^{\infty}|a_n||T_n(x)|.
\end{align*}
Using the coefficient estimate from the previous step,
\begin{align*}
|f(x)-S_Nf(x)|\le \sum_{n=N+1}^{\infty}2M\rho^{-n}.
\end{align*}
The geometric series satisfies
\begin{align*}
\sum_{n=N+1}^{\infty}\rho^{-n}=\frac{\rho^{-(N+1)}}{1-\rho^{-1}}=\frac{\rho^{-N}}{\rho-1}.
\end{align*}
Hence
\begin{align*}
|f(x)-S_Nf(x)|\le \frac{2M\rho^{-N}}{\rho-1}.
\end{align*}
Taking the supremum over $x\in[-1,1]$ gives
\begin{align*}
\|f-S_Nf\|_\infty\le \frac{2M\rho^{-N}}{\rho-1}.
\end{align*}
This is the asserted uniform convergence estimate.
[/step]