Let $U_i \sim Po(1)$ be independent, for $i=1,2,\ldots$ Set
$X_n = U_1 + U_2 + \cdots + U_n$
Then $X_n \sim Po(n)$
Why is that so? Is this valid for all distributions? Is there any formula we can use to get a general rule for all probabilities? I.e can we say that if the sum of independent Exp($1$) variables have the distribution Exp($n$) ?
This phenomenon is not valid for all distributions, but it does hold true for some. For example, the sum of independent Gamma distributions with the same rate parameter is also Gamma with the shape parameter being the sum of the individual shapes. The sum of independent Normal distributions is also Normal. The sum of independent Binomial variables with the same trial probability of success is also Binomial.
Mathematically, we would write these examples as follows: $$X_i \sim \operatorname{Gamma}(a_i, b); \quad f_{X_i}(x) = \frac{b^{a_i} x^{a_i - 1} e^{-b x}}{\Gamma(a_i)}$$ implies $$T_n = \sum_{i=1}^n X_i \sim \operatorname{Gamma}\left(\sum_{i=1}^n a_i, b\right).$$ For the Normal distribution, $$X_i \sim \operatorname{Normal}(\mu_i, \sigma_i^2) \quad \text{implies} \quad T_n \sim \operatorname{Normal}\left(\sum_{i=1}^n \mu_i, \sum_{i=1}^n \sigma_i^2\right).$$ And for the Binomial case, $$X_i \sim \operatorname{Binomial}(n_i, p) \quad \text{implies} \quad T_n = \operatorname{Binomial}\left(\sum_{i=1}^n n_i, p\right).$$ These are all provable using moment generating functions, or by an induction argument, as is also the case with the Poisson distribution. But it is not true for the exponential distribution example. Rather, because $$\operatorname{Exponential}(\lambda) = \operatorname{Gamma}(1,\lambda),$$ we can see from the above that the sum of $n$ IID exponential variables with common rate $\lambda$ will be Gamma distributed with shape $n$ and rate $\lambda$.
The proof for the Poisson case using MGFs is straightforward. We first calculate the MGF of a single Poisson variable $X$ with rate $\lambda$: $$M_X(t) = \operatorname{E}[e^{tX}] = \sum_{x=0}^\infty e^{tx} e^{-\lambda} \frac{\lambda^x}{x!} = \sum_{x=0}^\infty e^{\lambda(e^t-1)} e^{-\lambda e^t} \frac{(\lambda e^t)^x}{x!} = e^{\lambda(e^t-1)}.$$ Therefore, if $X_i \sim \operatorname{Poisson}(\lambda_i)$ are independent, the MGF of their sum is $$M_{T_n}(t) = \prod_{i=1}^n M_{X_i}(t) = \prod_{i=1}^n e^{\lambda_i (e^t - 1)} = \exp \left( (e^t-1) \sum_{i=1}^n \lambda_i \right),$$ which is clearly Poisson with rate $$\lambda_T = \sum_{i=1}^n \lambda_i.$$