Random variable $X$ has gamma distribution with parameters $\alpha = \frac{1}{3}, \ \beta = n,$ i.e. its pdf $p(y)=\frac{\frac{1}{3^n}}{\Gamma(n)}y^{n-1}e^{-\frac{y}{3}}$ when $y>0$, and $0$ otherwise. Find $\lim_{n \longrightarrow \infty} P(|X_n-3n|<x\sqrt{n}),$ and plot its graph.
$P(|X_n-3n|<x\sqrt{n}) = P(X_n \in (-x\sqrt{n}+3n, \ x\sqrt{n}+3n)) = $
$$ =\int_{-x\sqrt{n}+3n}^{x\sqrt{n}+3n} p(y) dy = \frac{1}{(n-1)!}\int_{-x\sqrt{n}+3n}^{x\sqrt{n}+3n} \left( \frac{y}{3} \right)^{n-1} e^{-\frac{y}{3}}d \left( \frac{y}{3} \right) = \left/ t:=\frac{y}{3} \right/ = $$
$$ =\frac{1}{(n-1)!} \left( \int_{0}^{-x\sqrt{n}+3n}t^{n-1} e^{-t}dt + \int_{-x\sqrt{n}+3n}^{x\sqrt{n}+3n} t^{n-1} e^{-t}dt -\int_{0}^{-x\sqrt{n}+3n}t^{n-1} e^{-t}dt \right) = \frac{1}{(n-1)!} \left( \int_{0}^{x\sqrt{n}+3n}t^{n-1} e^{-t}dt - \int_{0}^{-x\sqrt{n}+3n}t^{n-1} e^{-t}dt \right) = \frac{1}{(n-1)!} \left( \gamma(n, x\sqrt{n}+3n) - \gamma(n, -x\sqrt{n}+3n) \right) = \sum_{k=0}^{\infty} \frac{e^{-x\sqrt{n}-3n}(x\sqrt{n}+3n)^{k+n} - e^{x\sqrt{n}-3n}(-x\sqrt{n}+3n)^{k+n}}{\Gamma (n+k+1)} = $$ ($\frac{1}{(n-1)!}$ and $\Gamma (n)$ from $\gamma$s cancelled each other out.)
By Stirling's approximation,
$$ =\sum_{k=0}^{\infty} \frac{e^{-x\sqrt{n}-3n}(x\sqrt{n}+3n)^{k+n} - e^{x\sqrt{n}-3n}(-x\sqrt{n}+3n)^{k+n}}{\sqrt{2\pi (n+k)}(n+k)^{n+k}e^{-(n+k)}}. $$
I am stuck. Is this the wrong route? Have I made any initial mistakes?
Would be very grateful for any help on how to approach it.
There's nothing wrong per se, but you seem to have missed the main idea here, which is to invoke the CLT. Your notation is somewhat different from standard notation for the Gamma (in which the role of $\alpha$ and $\beta$ are exchanged), but I'll stick to it for the following.
Its a standard (and important) result that if $X \sim \mathrm{\Gamma}(\alpha, \beta_X), Y \sim \mathrm{\Gamma}(\alpha, \beta_Y),$ and $X$ and $Y$ are independent, then $X+Y \sim \Gamma(\alpha, \beta_X + \beta_Y)$ (see below).
So, let $Z_1, Z_2, \cdots$ be iid $\Gamma(1/3,1)$ random variables. Then the law of $X_n$ is the same as the law of $S_n := Z_1 + Z_2 + \cdots + Z_n.$ But, as a sum of iid random variables with finite variance, $S_n$ admits a central limit behaviour, i.e., $$ \frac{S_n - n \mathbb{E}[Z_1]}{\sqrt{n \mathrm{Var}(Z_1)}} \rightsquigarrow \mathcal{N}(0,1).$$ Here, notice that $\mathbb{E}[Z_1] = 3, \mathrm{Var}(Z_1) = 9.$ So, $$ P(|X_n - 3n|< x \sqrt{n}) = P\left(\frac{|X_n - 3n|}{\sqrt{9n}} < \frac{x}{3} \right) \\= P\left( \frac{|S_n - 3n|}{\sqrt{9n}} < \frac{x}{3} \right) = \Phi(x/3) - \Phi(-x/3) + o(1),$$ where $\Phi$ is the standard Gaussian cdf, and taking limits yields the answer $\Phi(x/3) - \Phi(-x/3) = 2\Phi(x/3) - 1$.
Easiest way to show the additivity is via characteristic functions. Assuming you know (or can show) that the characterisitc function of $\Gamma(\alpha, \beta)$ is $$ \varphi(t;\alpha, \beta) = \left(1 - \frac{it}{\alpha}\right)^{-\beta},$$ observe that since $X$ and $Y$ are independent, we have $$\varphi_{X+Y}(t) = \varphi_X(t)\varphi_Y(t) = \left(1 - \frac{it}{\alpha}\right)^{\beta_X} \cdot \left(1 - \frac{it}{\alpha}\right)^{\beta_Y} \\ = \left(1 - \frac{it}{\alpha}\right)^{\beta_X + \beta_Y} = \varphi(t;\alpha, \beta_X+\beta_Y),$$ and we're done.
Its worth noting that for $\beta = 1,$ the corresponding Gamma distribution is just an exponential distribution. So, for integer $\beta,$ $\Gamma(\alpha, \beta)$ is the distribution of a sum of $\beta$ iid exponential random variables. This is sometimes how the Gamma distribution is motivated, and is a reasonable way to remember the additivity.