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.
On
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!.$
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)$.