Finding probability distribution from given moments

Let us denote $q := ghk$ and your variable by $X.$ You can consider the moment-generating function of the variable $X,$ $M_X(t) = \mathbb{E}[e^{tX}] = \sum_{0 \leq k} \mathbb{E}[X^k] \frac{t^k}{k!},$ which in your case would be $$M_X(t) = \sum_{0 \leq k} 2^{-2k} \binom{2k}{k} \frac{(tq)^{2k}}{(2 k)!} = \sum_{k \geq 0} 2^{-2 k} \frac{(2 k)!}{(k!)^2} \frac{(t q)^{2k}}{(2k)!} = \sum_{k \geq 0} (\frac{t^k q^k}{2^k k!})^2.$$ Note that we already get that $M_X(t) \leq e^{tq},$ so this series is actually well defined in this case (it is not always the case with moment generating functions),. In fact, this is just a modified Bessel function of the first kind. More precisely, $M_X(t) = I_0(t q),$ so your intuition was right that we would get some kind of Bessel function somehow.

P.S.: We can determine the PDF of $X$ by using the inversion formula. Note that by the same reasoning as above, we deduce that the characteristic function of $X$ is given by $\varphi_X(t) = J_0(q t),$ a Bessel function of the first kind. It follows from the Inversion Theorem that $X$ admits a probability density function, which is given simply as the Fourier Transform of the characteristic function of $X.$ This has been computed before on this very site (see for example here or here). We have that $$\mathcal{F}(J_0)(\xi) = \chi_{(-1, 1)}(\xi)\frac1\pi \frac{1}{\sqrt{1 - \xi^2}}.$$ However, we are interested in a scaling of this function, so performing a simple substitution in the Fourier Transform integral, we get the expression of the PDF of $X:$ $$f_X(\xi) = \chi_{(-q, q)}(\xi) \frac{1}{\pi} \frac{1}{\sqrt{q^2 - \xi^2}}.$$

The arcsine distribution stretched to have support on $[-ghk, ghk]$ has the moments you are looking for. So, with $-ghk < x < ghk$, the probability density function is $$f(x)=\frac{1}{\pi \sqrt{{{(ghk)}^{2}}-{{x}^{2}}}}$$ and the cumulative distribution function is $$F(x)=\mathbb P(X\le x) = \frac1\pi \cos^{-1}\left(\frac{-x}{ghk}\right)$$

As a quick check in R of the first twelve moments, setting $ghk=2$ gives the central binomial coefficients:

qstretchedarcsine <- function(prob, ghk){ - ghk * cos(prob * pi) }
thisdata <- qstretchedarcsine(ppoints(10^6), ghk=2)
round(colMeans(outer(thisdata, 1:12, "^")), 12)
#    0   2   0   6   0  20   0  70   0 252   0 924

See https://math.stackexchange.com/a/2629048/6460 for the relationship between the arcsine distribution and Bessel functions