Confidence intervals around functions of estimated parameters Suppose I estimate the following model,
$$
 y = \hat{\beta_0} + \hat{\beta_1} X_1 + \hat{\beta_2} X_2\ ,
$$
and am additionally interested in the quantity
$$
 \gamma = \sqrt{\beta_0^2 - 4\beta_1\beta_2}\ .
$$
An estimate of the above, $\hat{\gamma}$, is clearly given by inserting the estimated $\hat\beta_i$'s in the above formula.
But how does one go about computing the confidence interval around such a "compound" quantity that is not directly estimated?
My attempt.
My "empirical" approach is as follows. I generate a multivariate normal distribution of size $n$ using the estimated means $\hat\beta_i$'s and the covariance matrix of the fit. I consequently have a distribution of size $n$ for the $\gamma$. The mean of this distribution $\bar\gamma$ is naturally then close to the $\hat\gamma$ computed above, and I now additionally have a standard deviation of this distribution $\sigma$.
Using $\sigma$ to generate a confidence interval.
If I use the standard definition $\bar\gamma \pm t(\alpha, n-1) \sigma / \sqrt{n}$, where $t$ is the PPF of the $t$-distribution, I run into the following quandary: since I pick $n$ myself, I can in fact make the standard error as small as I like. This is not what I want: I am looking for a "natural" confidence interval which depends solely on the uncertainty in my estimates of the $\beta$'s. Therefore, my idea is the following: $\sigma$ stays positive, and seems to "converge" to its true value for large $n$: hence, 95% of the values in my $\gamma$ distribution lie within the interval $\bar\gamma \pm 2\sigma$.
Is this approach valid?

There's one big problem with your proposed method, although it can be fairly easily remedied. In particular, you compute the estimate of variance of $\gamma$ as $\frac{V(\gamma_i)}{n}$ (where $V(\gamma_i)$ is variance of your samples generated by first sampling from the multi variance normal and then pushing through your function), when the standard error should really just be $V(\gamma_i)$.
But furthermore, I would suggest you use the simulated percentiles to create the CI rather than using the mean + se formula.
And one final note, we note that the sqrt function is not defined for negative values. This means if a random draw gives you $4\beta_1 \beta_2 > \beta_0^2$, your estimate is not real valued. Assuming that this is not what you want to do, I would say you should consider simply dropping the simulated values when this is true. That may seem adhoc, but it perfectly aligns with defining a restricted parameter space where $4\beta_1 \beta_2 \le \beta_0^2$.
