What is the mean of this exponential random variable? I have read a paper that says that the following is exponentially distributed
$$ Y=  \bigl| \sum_{i=1}^n \gamma_i^{-\frac{1}{2}} h_i \bigl|^2$$
where $\gamma_i$ are non-negative constants and 
$$h_i \sim \mathcal{CN} (0,1)\,,$$
where $\mathcal{CN}(\mu,\Gamma)$ denotes a circularly symmetric complex normal distribution.
The authors claim that the mean of this exponential distributed $Y$ is $$\sum_{i=1}^n{\gamma_i}^{-1}$$. 
Does anyone know why?
My thoughts are the following
1- The sum of Gaussian is also Gaussian 
2- Magnitude of Gaussian is Rayleigh distributed
3- Taking the square of Rayleigh is exponential
But how do we get that the mean ?
My second part of this question is 
What happens if $$h_{i} \sim \text{Nakagami } m $$
 A: Just addressing the first part for now:
Your line of thinking about why it's exponential is essentially along the right lines, after small modifications to get some details correct (you'll need independence in your step 2 for example, to invoke the Rayleigh).
For constant $c_i$, I think we have that $c_i h_i \sim \mathcal{CN}(0,c_i^2I)$ (where here $\mathcal{CN}$ with two arguments refers to the circularly symmetric complex normal). 
This along with linearity of expectation explains where the $\sum \gamma_i^{-1}$ comes from.

Answer to Q from comments:

how does the scale parameter $γ^{−1}_i$ change the distribution of $|h_i|^2$?

If $h_i^2 \sim \text{Gamma}(\alpha, \beta)$ (where here I mean the shape-scale parameterization, not the shape-rate parameterization), then 
$γ^{−1}_i\,h_i^2\sim\text{Gamma}(\alpha, γ^{−1}_i\beta)$.
(If you want the shape-rate parameterization, you multiply the rate parameter by $\gamma_i$)
A: A standard circularly symmetric complex normal random variable is
of the form $Z = X + iY$ where $X$ and $Y$ are independent zero-mean
normal random variables with variance $\frac 12$. Here $i = \sqrt{-1}$.
We write $Z \sim \mathcal{CN}(0,1)$ and $X, Y \sim \mathcal{N}(0,\frac 12)$.
Note that $$E[|Z|^2] = E[X^2+Y^2]= \frac 12 + \frac 12 =1.$$
Note also that 
$|Z|^2$ is an exponential random variable with parameter (and hence mean)
equal to $1$.
From this, it is easy to deduce that $\sigma Z \sim \mathcal{CN}(0,\sigma^2)$
and that $|\sigma Z|^2$ is an exponential random variable with mean $\sigma^2$
and hence parameter $\frac{1}{\sigma^2}$.
Thus, since the $H_k$ are independent $\mathcal{CN}(0,1)$ random variables, meaning that $H_k = X_k + i Y_k$ where $X_k$ and $Y_k$ are independent 
$\mathcal N(0,\frac 12)$ random
variables, $\gamma_k^{-1/2}H_k$ is a
$\mathcal{CN}(0,\gamma_k^{-1})$ random variable, and so
$$Z = \sum_{k=1}^N \gamma_k^{-1/2}H_k = \sum_{k=1}^N\gamma_k^{-1/2}(X_k+iY_k)
= \bar{X} + i\bar{Y} \sim \mathcal{CN}\left(0,\sum_k \gamma_k^{-1}\right).$$
It follows that $|Z|^2 = \displaystyle \left|\sum_{k=1}^N \gamma_k^{-1/2}H_k\right|$ 
is an exponential random variable with mean $\sum_k \gamma_k^{-1}$.
About the only claim that does not follow from simple applications
of general rules is that if $X$ and $Y$ are independent zero-mean 
normal random variables with common variance $\frac 12\sigma^2$,
then $X^2+Y^2$ has an exponential distribution. Note that
the usual method changing to polar coordinates gives
$$P\{X^2+Y^2 > z\}
=\int_{\sqrt{z}}^\infty \int_0^{2\pi}
 \frac{1}{\pi\sigma^2} r\exp(-r^2/\sigma^2) \,\mathrm d\theta \, \mathrm dr
= \exp(-z/\sigma^2)$$
showing that $X^2+Y^2$ has an exponential distribution with parameter
$\sigma^{-2}$ and hence mean $\sigma^2$.
