Uniform rate of convergence of $\sum_{k=0}^{n} {2n+1\choose{k}}\left(a^{k+1}(1-a)^{(2n+1)-k}+a^{(2n+1)-k}(1-a)^{k+1}\right)$ to $\min(a,1-a)$ This question is related to this one. There, it is proved that
$$\forall a \in [0,1], f_{2n+1}(a)\to\min(a,1-a), n\to\infty$$
where
$$\forall n\in\mathbb{N}, \forall a\in[0,1], f_{2n+1}(a) :=\sum_{k=0}^{n} {2n+1\choose{k}}\left(a^{k+1}(1-a)^{(2n+1)-k}+a^{(2n+1)-k}(1-a)^{k+1}\right).$$
However, the aforementioned result is qualitative in spirit, and now I'm looking for a quantitative version of this results. Specifically, can we show that
$$\sup_{a\in[0,1]} |f_{2n+1}(a)-\min(a,1-a)| \to 0, n\to +\infty?$$
If this is the case, at which rate this quantity approaches zero?
 A: By symmetry (w.r.t. $a\mapsto 1-a$), we may consider $a\in[0,1/2]$. From my answer to the LQ, $$f_{2n+1}(a)-a=A_n(1-2a)\int_0^a x^n(1-x)^n\,dx,\qquad A_n:=(2n+1)\binom{2n}{n},$$ so that the maximum is attained at the solution (see below) $a=a_n$ of $$(1-2a)a^n(1-a)^n=2\int_0^a x^n(1-x)^n\,dx.$$
Now we substitute $x=a(1-t)$ and put $a=\lambda/(1+\lambda)$, which gives the equation $$\frac{1-\lambda}{2\lambda}=\int_0^1(1-t)^n(1+\lambda t)^n\,dt;$$ now we see that the solution $\lambda=\lambda_n\in(0,1)$ exists (and is unique). Further, $$\int_0^1(1-t)^n(1+\lambda t)^n\,dt<\int_0^1(1-t^2)^n\,dt=2^{2n}/A_n=\mathcal{O}(n^{-1/2})$$ implies $\lambda_n=1-\mathcal{O}(n^{-1/2})$, and our maximum value is $$M_n:=\max_{0\leqslant a\leqslant 1/2}\big(f_{2n+1}(a)-a\big)=\frac{A_n\lambda_n^n(1-\lambda_n)^2}{2(1+\lambda_n)^{2n+2}}=\mathcal{O}(n^{-1/2}).$$

The exact coefficients are harder to compute. The result is as follows. Let $\alpha$ be the solution of $$\alpha=\int_0^\infty e^{-t^2-2\alpha t}\,dt\color{LightGray}{\left[=\frac{\sqrt\pi}{2}e^{\alpha^2}\operatorname{erfc}\alpha\right]}.$$ Then $\lambda_n\asymp 1-2\alpha/\sqrt{n}$, and $M_n\asymp\beta/\sqrt{n}$ with $\beta=\alpha^2 e^{-\alpha^2}/\sqrt\pi$.
Numerically, $2\alpha\approx 1.06319377\ldots$ and $\beta\approx 0.12018779\ldots$
A: In @metamorphy's nice and elegant answer appears the interesting equation (I cannot resist an equation !)
$$\alpha=\int_0^\infty e^{-t^2-2\alpha t}\,dt$$ which is
$$\frac{\sqrt{\pi }}{2}  e^{\alpha ^2} \text{erfc}(\alpha )-\alpha=0$$ which, as mentioned, requires some  numerical method (with Newton, is is just a piece of cake starting with $\alpha_0=0$.
What can also be done is a series axpansion around $\alpha=0$ and use series reversion to get things like
$$\alpha=\frac{\sqrt \pi} 4 \left(1+\frac{\pi }{24}+\frac{31 \pi ^2}{7680}+\frac{323 \pi ^3}{645120}+\frac{104161 \pi
   ^4}{1486356480}+\frac{6893851 \pi ^5}{653996851200}+\frac{2715451111 \pi
   ^6}{1632376140595200}+\cdots\right)$$
