When $X_k$ is a random variable with Binomial distribution $B(k;p)$ its mean is $kp$ and variance is $kp(1-p)$.

So, in order to normalise the distribution, we define $Y_k=(X_k-kp)/\sqrt{kp(1-p)}$.

Putting $r=kp+x\sqrt{kp(1-p)}$, see that the distribution of $Y_k$ is given by

$$P(Y_k=x)=\left(\genfrac{}{}{0ex}{}{k}{r}\right)p^r(1-p{)}^{k-r}$$

sage: def cenbinom(k,p,y):
...       r = round(k*p+y) # we will put y=x*sqrt(k*p*(1-p))
...       return N(binomial(k,r)*p^r*(1-p)^(k-r))

For $\mid x\mid \le A$ and $k$ sufficiently large the de Moivre Laplace theorem says that

$$P(Y_k=x)\simeq \frac{1}{\sqrt{2\pi kp(1-p)}}{e}^{-x^2/2}$$

We try to see this in the following images.

sage: p=0.8;a=sqrt(p*(1-p))
sage: normpl=plot(exp(-x^2/2),(x,-3,3),color='orange')
sage: plseq=[normpl+point(map(lambda x:(x,k*cenbinom(k^2,p,x*k*a)), srange(-3,3,1/(k*a))),color='gray') for k in range (2,40)]

sage: planim=animate(plseq,xmin=-3,xmax=3,ymin=0,ymax=1,figsize=[6,4])

sage: planim.show(delay=50)

This property is not limited to the case of a coin flip (with uneven probability).

Let $T_i$ be a sequence of independent identically distributed variables with expectation $m$ and variance ${\sigma }^2$.

Let $Y_n={\sum }_{i=1}^n{T}_i$. Then $E(Y_n)=nm$ and ${\sigma }^2(Y_n)=n{\sigma }^2$.

Then, the distribution of $X_n=\frac{Y_n-nm}{\sqrt{n{\sigma }^2}}$ approaches the normal distribution as $n$ approaches infinity. This is called the Central Limit Theorem.

sage: n=50;m=200
sage: ylist=[sum(random() for _ in range(n)) for _ in range(m)] # we take m=200 samples of Y_n

The function random() is a uniformly distributed random variable in the interval $[0,1]$. Hence its mean is $1/2$ and variance is $1/12$.

sage: xlist=map(lambda y: (y-n/2)/N(sqrt(n/12)), ylist)

We count the number of values of $x$ less than or equal to a given value to get the distribution function.

sage: distfvals=[(t,len([x for x in xlist if x<=t])/m) for t in srange(-3,3,0.01)]

sage: grf1=point(distfvals,color='gray')

We calculate the normal distribution values for the same range of values for $t$.

sage: c=N(1/sqrt(2*pi));var('t,x')
sage: def distnorm(x):
...       s=N(sqrt(2))
...       return 1/2 + erf(x/s)/2
sage: grf2=plot(distnorm(x),(x,-3,3),color='orange')

sage: g=grf1+grf2
sage: g.show(xmin=-3,xmax=3,ymin=0,ymax=1,figsize=[6,3])