Suppose that I have $N$ i.i.d. random variables with a distribution $Q$, which has mean around 1. That is $R_1, R_2,\ldots,R_N \sim Q$. I would like $\prod_{i=1}^N R_i$ to approximate a normal random variable.
What are some candidate distributions for $Q$ so that I get a reasonable approximation as $N$ gets larger?
So long a partial answer, using the Mellin transform. We use the definition of Mellin transform of a probability density on the non-negative real axis $$ \DeclareMathOperator{\M}{\mathcal{M}} \M(X)(s) = \int_0^\infty x^s f(x) \; dx $$ where $f$ is the density function of the random variable $X$.
Your random variables are defined on the real axis, but are symmetric there, so we reduce to the non-negative axis by multiplication with Rademacher variables, that is, random variables taking the values $\pm 1$ with equal probability. So write $R_i = S_i T_i$, where the factors are independent, $S_i$ (sign) is Rademacher and $T_i$ is standard half-normal (the absolute value of a standard normal).
Then $\prod_1^n R_i = \prod_1^n S_i \prod_1^n T_i$ and since the product of $n$ independent rademacher variables again is Rademacher, we get $$ \prod_1^n R_i = S \prod_1^n T_i $$ where $S$ is Rademacher. Since a standard normal random variable is the independent product of a Rademacher variable and a half-normal random variable, your question is reduced to ask for a distribution such that the product of $n$ independent copies is half-normal.
The Mellin transform of a half-normal variable $N$ is $$ \M(N)(s) = \int_0^\infty x^s \sqrt{2/\pi} e^{-x^2/2}\; dx \\ = \frac{2^{s/2} \Gamma(\frac{s+1}{2})}{\sqrt{\pi}} $$ and the Mellin transform of $T_i$ is then the $n$th root of this $$ \M(T_i)(s) = \left( \frac{2}{\pi} \right)^{s/2n} \Gamma\left(\frac{s+1}{2}\right)^{1/n}. $$ But the question remains how to invert this, exactly or approximately ...
