If you wanted to get a closed-form solution for the posterior distribution then you would use a conjugate prior, in which case $\mu$ and $\sigma^2$ would not be independent. Since you are using an MCMC algorithm to produce an approximation to the posterior distribution, there is no need to have a conjugate prior. I believe the prior in your code is intended to approximate a flat prior, $p(\mu,1/\sigma^2) \propto 1$.
But suppose you wanted to use an MCMC algorithm to simulate draws from the posterior using the conjugate prior just for the heck of it. Then you would write the code for the joint distribution for $(\mu,\sigma^2)$, which is not the product of the two marginal distributions. (You are correct that the marginal prior distribution for $\mu$ would be $t$.)
Yes it does, the Generalized Beta prime distribution with shape parameter equal to 1.
We can get there fairly easily by integrating $\beta$ out of the joint distribution of $x$ and $\beta$:
$$f(x,\beta|\alpha, \alpha_0, \beta_0) = {\beta^{\alpha}x^{\alpha-1} \over \Gamma(\alpha)}e^{-\beta x}{\beta_0^{\alpha_0}\beta^{\alpha_0-1} \over \Gamma{\alpha_0}}e^{-\beta_0\beta}$$
Rearranging terms and ignoring everything that isn't related to either $\beta$ or $x$ (as they will all be handled by the constant of integration at the end) results in a great deal of simplification:
$$f(x,\beta|\cdot) \propto x^{\alpha-1}\beta^{\alpha+\alpha_0-1}e^{-(\beta_0+x)\beta}$$
We can integrate out $\beta$ easily enough by noting that the two terms involving $\beta$ are those of a Gamma-distributed variate with shape parameter $\alpha + \alpha_0$ and rate parameter $\beta_0 + x$, so the integral must equal the inverse of the constant of integration of such a distribution:
$$x^{\alpha-1}\int \beta^{\alpha+\alpha_0-1}e^{-(\beta_0+x)\beta}\text{d}\beta = {x^{\alpha-1} \Gamma(\alpha+\alpha_0) \over(\beta_0 + x)^{\alpha + \alpha_0}} \propto x^{\alpha-1} (\beta_0 + x)^{-\alpha - \alpha_0}$$
A slight rearrangement of terms and some minor algebra gets us to:
$$f(x|\cdot) \propto \left({x \over \beta_0}\right)^{\alpha-1} \left(1 + {x\over \beta_0}\right)^{-\alpha - \alpha_0}$$
which clearly matches the formula for the GBPD
(with shape parameter $p=1$) as given by Wikipedia and reproduced here:
$$f(x;\alpha,\beta,p,q) = \frac{p \left(\frac x q \right)^{\alpha p-1} \left(1+ \left(\frac x q \right)^p\right)^{-\alpha -\beta}}{qB(\alpha,\beta)}$$
Best Answer
I believe that the problem is with your guess that the inverse gamma is so easily extended to the multivariate case.
From the distributions appendix of Gleman et al, Bayesian Data Analysis (3rd Edition), 582
I'm uncertain whether you'd like to proceed in your own investigation with this hint, or if you'd like me to spill the beans and post a full solution. (Though, turning to page 73 of the same text, we find the particular underlying algebra that you're interested in.)