Does standard distance follow the 68-95-99.7 rule?

I'm wanting to do a simple standard distance demonstration for my students in R, but I've come across a conundrum. When I simulate the creation of 10,000 points in a spatial normal distribution, nearly all data fall within 2 standard distances of the mean (far more than the expected ~95%). In fact, it seems as though about 63% falls within 1 standard distance and about 98% falls within 2 standard distances (pretty consistent results through several iterations).

From my understanding, points following a spatial normal distribution also follow the 68-95-99.7 rule with respect to standard distance (source). What's going on?

Here's what I'm using to calculate standard distance and create some visualization:

library(plotrix)

## create function for standard distance
st.dist <- function(x,y) {
sqrt(sum(((x - mean(x))^2)/length(x)) + (sum((y - mean(y))^2)/length(y)))
}
## assign values to x
x <- rnorm(10000)

## assign values to y
y <- rnorm(10000)

## plot them
plot(x, y,
xlim = c(mean(x) - 4*st.dist(x,y),
mean(x) + 4*st.dist(x,y)),
ylim = c(mean(y) - 4*st.dist(x,y),
mean(y) + 4*st.dist(x,y)))

## draw standard distances as circles
draw.circle(x = mean(x),
y = mean(y),

draw.circle(x = mean(x),
y = mean(y),

draw.circle(x = mean(x),
y = mean(y),

## create a data frame for the values
df <- data.frame(x=x, y=y, dist=NA)

## calculate distance from mean center to each point
df$$dist <- sqrt((mean(x) - df$$x)^2 + (mean(y) - df$y)^2) Checking the percentage of points within 1, 2, and 3 standard distances of the mean center: > ## percentage of points falling within one standard distance of the mean center > nrow(df[df$$dist > ## percentage of points falling within two standard distances of the mean center > nrow(df[df$$dist<2*st.dist(x,y),])/nrow(df)  0.9804 > > ## percentage of points falling within three standard distances of the mean center > nrow(df[df$dist<3*st.dist(x,y),])/nrow(df)
 0.9998
>

Yet, the x and y values on their own seem to be behaving properly:

## percentage of x coordinates between -1 and 1
> length(x[x<1 & x>-1])/length(x)
 0.6769
>
> ## percentage of x coordinates between -2 and 2
> length(x[x<2 & x>-2])/length(x)
 0.9545
>
> ## percentage of x coordinates between -3 and 3
> length(x[x<3 & x>-3])/length(x)
 0.9964
>
> ## percentage of y coordinates between -1 and 1
> length(y[y<1 & y>-1])/length(y)
 0.6878
>
> ## percentage of y coordinates between -2 and 2
> length(y[y<2 & y>-2])/length(y)
 0.9514
>
> ## percentage of y coordinates between -3 and 3
> length(y[y<3 & y>-3])/length(y)
 0.9977
>

Short answer: definitely not true. Why? For that we need the distribution of standard distance. So let $$(X_1,Y_1), \dotsc, (X_n,Y_n)$$ be iid observations from a bivariate normal distribution $$\mathcal{bivnorm}(\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)$$. Let $$D$$ be the $$n\times 2$$ data matrix (with $$X_i$$, $$Y_i$$ observations as columns) and $$H=I_n - n^{-1} 11^T$$ the $$n\times n$$ centering matrix (a projection matrix), $$1=1_n$$ is the vector with only $$n$$ 1's. Then we can write $$\text{stdist}= \sqrt{\text{trace}(n^{-1} D^T H D))}$$. So what is the distribution? I will come back for a full (I hope) answer, but for now some references: Traces and cumulants of quadratic forms in normal variables.
But for the very specific case $$\rho=0$$ this is simple, then we have $$\text{stdist}= \sqrt{\frac{n-1}{n}(S_x^2 + S_y^2)}=\sqrt{\frac{1}{n}}\sqrt{(n-1)S_x^2 + (n-1)S_y^2}$$ and the sample variances are in this case independent, so the sum under the square root sign in the last factor have a chisquared distribution with $$2n-2$$ degrees of freedom. That will answer your question in that case.