I've been led to believe (see here and here) that Mahalanobis distance is the same as the Euclidian distance on the PCA-rotated data. In other words, taking multivariate normal data $X$, the Mahalanobis distance of all of the $x$'s from any given point (say $\mathbf{0}$) should be the same as the euclidian distance of the entries of $X^{rot}$ from $\mathbf{0}$, where $X^{rot}$ is the product of the data and the PCA rotation matrix.

1. Is this true?

My code below is suggesting to me that it is not. In particular, it looks like the variance of the mahalanobis distance around the PCA-euclidian distance is increasing in the magnitude of the PCA-euclidian distance. Is this a coding error, or a feature of the universe? Does it have to do with imprecision in an estimate of something? Something that gets squared?

N=1000
cr = runif(1,min=-1,max=1)
A = matrix(c(1,cr,cr,1),2)
e<-mvrnorm(n = N,rep(0,2),A)
mx = apply(e, 2, mean)
sx = apply(e, 2, sd)
e = t(apply(e,1,function(X){(X-mx)/sx}))
plot(e[,1],e[,2])
dum<-rep(0,2)
md = mahalanobis(e,dum,cov(e))

pc = prcomp(e,center=F,scale=F)
d<-as.matrix(dist(rbind(dum,pc$x),method='euclidian',diag=F))
d<-d[1,2:ncol(d)]
plot(d,md^.5)
abline(0,1)

2. Assuming that the answer to the above is true, can one use the PCA-rotated euclidian distance as a stand-in for the mahalanobis distance when P>N?

If not, is there a similar metric that captures multivariate distance, scaled by correlation, and for which distributional results exist to allow the calculation of the probability of an observation?