Where does poisson distribution functional form come from? I know that the punctual probability function of a random variable $X$ with a Poisson distribution is:
$$P(X=k)= e^{-\lambda }\frac{\lambda ^{k}}{k!}.$$
Also, I've learned that the formula can be found as a limit of the binomial distribution. However, I was wondering if there is another way to obtain the formula without using the binomial distribution. I guess that the $e^{-\lambda }$ is used to normalize the probability, but I can't work out the idea behind the $\dfrac{\lambda ^{k}}{k!}$ part.
 A: According to the Wikipedia article Poisson distribution 

the sum of two independent random variables is Poisson-distributed, then so are each of those two independent random variables

So assume $\ X_1\ $ is Poisson with parameter $\ \lambda_1\ $ and $\ X_2\ $ is Poisson with parameter $\ \lambda_2.\ $ Let $\ X_0 = X_1 + X_2\ $ be the sum of the two variables.
By definition of sum of random variables we have $\ \lambda_0 = \lambda_1+\lambda_2\ $ and $\ P(X_0 = 0) = P(X_1 = 0)P(X_2 = 0).\ $ Let
$\ f(\lambda_i) = P(X_i=0).\ $ Then $\ f(\lambda_1+\lambda_2) = f(\lambda_1)f(\lambda_2).\ $ This is a well-known Cauchy functional equation with solution 
$\ f(\lambda) = f(1)^\lambda.\ $ Next we look at the value of
$$ P(X_0 = 1) = F(X_1 = 0)P(X_2 = 1) + F(X_1 = 1)P(X_2 = 0). $$
Divide both sides by $\ P(X_0 = 0)\ $ to get
$$ P(X_0 = 1) / P(X_0 = 0) = P(X_2 = 1)/P(X_2 = 0) + P(X_1 = 0)/P(X_1 = 0). $$
If $\ f(\lambda_i) := P(X_i = 1) / P(X_i = 0),\ $ then this a well-known functional equation with solution $\ f(\lambda) = c\lambda.\ $
To continue, we should take advantage of the ordinary generating function of a random variable
$$ G_X(z) := \sum_{n=0}^\infty P(X = n)\ z^n. $$
The special property of this is that $\  G_{X+Y}(z) = G_X(z)G_Y(z).\ $ Thus
$\ G_{X_0}(z) = G_{X_1}(z)G_{X_2}(z).\ $
Let $\ G_{\lambda_i}(z) := G_{X_i}(z),\ $ then
 $$ G_{\lambda_0}(z) = G_{\lambda_1+\lambda_2}(z) = G_{\lambda_1}(z)G_{\lambda_2}(z). $$ This is the Cauchy equation again and the solutions to it  uniquely determine the Poisson distribution. More precisely,
$$ G_\lambda(z) = (\exp(z)/\exp(1))^\lambda = \exp (\lambda z)/\exp (\lambda)
= \sum_{n=0}^\infty = \exp(-\lambda) (\lambda\ z)^n/n!. $$
This now gives us the expected $\ P(X = n) = e^{-\lambda} \lambda^n/n!.$
A: To see that $X$ is a valid probability mass function note that
$$
\sum_{k=0}^\infty P(X=k)=\sum_{k=0}^\infty e^{-\lambda}\frac{\lambda^k}{k!}=e^{-\lambda}\sum_{k=0}^\infty\frac{\lambda^k}{k!}=e^{-\lambda}\times e^{\lambda}=1
$$
and obviously $P(X=k)\geq 0$. 
As you also noted if $X_n\sim\text{Bin}(n, \lambda/n)$, then $X_n\stackrel{d}{\to } Y$ where $Y\sim \text{Poi}(\lambda).$
Another way is to use poisson processes/renewal processes. Let $X_i\stackrel{\text{i.i.d}} {\sim} \text{Exp}(1)$ be exponential random variables with mean one and put $T_n=X_{1}+\dotsb+X_n$ and put $N(t)=\max\{n\mid T_n\leq t\}$. Then we have the relation
$$
N(t)\geq n\iff T_n\leq t\tag{1}
$$
which allows us to compute $P(N(t)=k)$ as  the event $(N(t)=k)=(N(t)\geq k)\setminus (N(t)\geq {k+1})$ and we can use (1). It turns out that after some computation that $N(t)\sim \text{Poi}(t)$.
