Here's my understanding of your question: you have a circular covariance matrix, $$C = \left( \begin{array}{cc} \sigma^{2} & 0 \\ 0 & \sigma^{2} \end{array} \right)$$ together with a position vector $ \overrightarrow{X} = [x,y]$ and a 2D Gaussian pdf defined by those two items, i.e., $$p(x,y) = \frac{1}{2\pi\sigma^{2}} e^{-\frac{X^{T} C^{-1} X}{2}} dx dy = \frac{1}{2\pi\sigma^{2}} e^{-\frac{(x^{2} + y^{2})}{2 \sigma^{2}}} dx dy$$ and you want to know what encircling radius $r$, in polar coordinates, will enclose 90% of the volume inside the pdf.
If that is indeed your question, you may answer it by first transforming to polar coordinates, and then performing a pair of integrals. Using the fact that this transformation requires a Jacobian (see example 3 in the link) and the fact that $r^{2} = x^{2} + y^{2}$, then in polar coordinates, the equivalent expression becomes: $$p(r,\theta) = \frac{1}{2 \pi \sigma^{2}} e^{-r^{2}/2 \sigma^{2}} r dr d \theta $$ You can marginalize this over $ \theta $ to get $$ p(r) = \int_{\theta=0}^{\theta=2\pi} \frac{1}{2 \pi \sigma^{2}} e^{-r^{2}/2 \sigma^{2}} r dr d \theta = \frac{r}{\sigma^{2}} e^{-r^{2} / 2 \sigma^{2}} dr $$ which you may recognize as a Rayleigh distribution. The Rayleigh distribution has a cdf: $$ p_{c}(r) = \int_{0}^{r} \frac{z}{\sigma^{2}} e^{-z^{2} / 2 \sigma^{2}} dz = 1 - e^{-r^{2} / 2 \sigma^{2}} $$
You can solve for the 90% encircling radius by setting: $$ 0.9 = 1 - e^{-r^{2} / 2 \sigma^{2}}$$ or $$ r = \sqrt{-2 \ln (0.1)} \sigma \approx 2.15 \sigma $$ so your first instinct appears to have been the correct one.
Phrased a little differently, the marginal approach is wrong because a square measuring $2 \times 1.64 \sigma$ on each side will cover a smaller area than a circle of radius $2.15 \sigma $ and will enclose a volume consisting of less than 90% of the probability. Thus, the marginal approach is wrong.
Alternatively, you can think of it this way: if you assign a coverage area of $\pm 1.64 \sigma $ in each dimension, then the probability for each coordinate of the $(x, y)$ arrival vector to lie separately within the coverage is indeed 0.9, but the joint probability for both components to fall simultaneously within the coverage area is smaller than that: it's $0.9^{2} = 0.81$.
As you mentioned in your post we know the distribution of the estimate of $\widehat{r_{true}}$ if we are given $\mu$ so we know the distribution of the estimate $\widehat{r^2_{true}}$ of the true $r^2$.
We want to find the distribution of $$\widehat{r^2} = \frac{1}{N}\sum_{i=1}^N (x_i-\overline{x})^T(x_i-\overline{x})$$ where $x_i$ are expressed as column vectors.
We now do the standard trick
$$\begin{eqnarray*}
\widehat{r^2_{true}} &=& \frac{1}{N}\sum_{i=1}^N(x_i - \mu)^T(x_i-\mu)\\
&=& \frac{1}{N}\sum_{i=1}^N(x_i-\overline{x} + \overline{x} -\mu)^T(x_i-\overline{x} + \overline{x}-\mu)\\
&=&\left[\frac{1}{N}\sum_{i=1}^N(x_i - \overline{x})^T(x_i-\overline{x})\right] + (\overline{x} - \mu)^T(\overline{x}-\mu) \hspace{20pt}(1)\\
&=& \widehat{r^2} + (\overline{x}-\mu)^T(\overline{x}-\mu)
\end{eqnarray*}
$$
where $(1)$ arises from the equation
$$\frac{1}{N}\sum_{i=1}^N(x_i-\overline{x})^T(\overline{x}-\mu) = (\overline{x} - \overline{x})^T(\overline{x} - \mu) = 0$$
and its transpose.
Notice that $\widehat{r^2}$ is the trace of the sample covariance matrix $S$ and $(\overline{x}-\mu)^T(\overline{x}-\mu)$ only depends only on the sample mean $\overline{x}$. Thus we have written
$$\widehat{r_{true}^2} = \widehat{r^2} + (\overline{x}-\mu)^T(\overline{x}-\mu)$$
as the sum of two independent random variables. We know the distributions of the $\widehat{r^2_{true}}$ and $(\overline{x} - \mu)^T(\overline{x}-\mu)$ and so we are done via the standard trick using that characteristic functions are multiplicative.
Edited to add:
$||x_i-\mu||$ is Hoyt so it has pdf
$$f(\rho) = \frac{1+q^2}{q\omega}\rho e^{-\frac{(1+q^2)^2}{4q^2\omega} \rho^2}I_O\left(\frac{1-q^4}{4q^2\omega} \rho^2\right)$$
where $I_0$ is the $0^{th}$ modified Bessel function of the first kind.
This means that the pdf of $||x_i-\mu||^2$ is
$$f(\rho) = \frac{1}{2}\frac{1+q^2}{q\omega}e^{-\frac{(1+q^2)^2}{4q^2\omega}\rho}I_0\left(\frac{1-q^4}{4q^2\omega}\rho\right).$$
To ease notation set $a = \frac{1-q^4}{4q^2\omega}$, $b=-\frac{(1+q^2)^2}{4q^2\omega}$ and $c=\frac{1}{2}\frac{1+q^2}{q\omega}$.
The moment generating function of $||x_i-\mu||^2$ is
$$\begin{cases}
\frac{c}{\sqrt{(s-b)^2-a^2}} & (s-b) > a\\
0 & \text{ else}\\
\end{cases}$$
Thus the moment generating function of $\widehat{r^2_{true}}$ is
$$\begin{cases}
\frac{c^N}{((s/N-b)^2-a^2)^{N/2}} & (s/N-b) > a\\
0 & \text{else}
\end{cases}$$
and the moment generating function of $||\overline{x} - \mu||^2$ is
$$\begin{cases}
\frac{Nc}{\sqrt{(s-Nb)^2-(Na)^2}} = \frac{c}{\sqrt{(s/N-b)^2-a^2}} & (s/N-b) > a\\
0 & \text{ else}
\end{cases}$$
This implies that the moment generating function of $\widehat{r^2}$ is
$$\begin{cases}
\frac{c^{N-1}}{((s/N-b)^2-a^2)^{(N-1)/2}} & (s/N-b) > a\\
0 & \text{ else}.
\end{cases}$$
Applying the inverse Laplace transform gives that $\widehat{r^2}$ has pdf
$$g(\rho) = \frac{\sqrt{\pi}Nc^{N-1}}{\Gamma(\frac{N-1}{2})}\left(\frac{2\mathrm{i} a}{N\rho}\right)^{(2 - N)/2} e^{b N \rho} J_{N/2-1}( \mathrm{i} a N \rho).$$
Best Answer
You could use a truncated normal distribution. It's just a normal distribution that you only consider an interval for. You need to rescale it to make sure that the pdf integrates to 1. But this sounds to me to be exactly what you're looking for.