If you assume that $X_1$ and $X_2$ have the same variance $\sigma^2$, namely that $\Sigma = \sigma^2\begin{pmatrix}
\;1 & \rho\\
\rho & \;1 \end{pmatrix}$, then you can see by diagonalizing $\Sigma$ that $x_1 := X_1+X_2$ and $x_2 :=X_1-X_2$ are independent with variances $2\sigma^2(1 \pm \rho)$ respectively.
Furthermore from the properties of the covariance of the sum we have: $$\text{svar}(X_1) + \text{svar}(X_2) = \frac{1}{2}(\text{svar}(x_1) + \text{svar}(x_2))$$
(I use $\text{svar}$ to distinguish sample variance from variance)
Therefore you can express $V_D$ and $V_p$ in term of the two independent variables $s_1^2 = \text{svar}(x_1)$, $s_2^2 =\text{svar}(x_2)$ such that
$$ V_D = s_2^2 , \quad V_p = \frac{1}{4}(s_1^2 + s_2^2) $$
Where by the well known distribution of the sample variance we know that
$$ \frac{n-1}{2\sigma^2(1+\rho)}s_1^2 \sim \chi^2_{n-1}, \quad \frac{n-1}{2\sigma^2(1-\rho)}s_2^2 \sim \chi^2_{n-1}$$
or equivalently $s_{1,2}^2 \sim \text{Gamma}((n-1)/2,4\sigma^2(1 \pm \rho)/(n-1))$.
Your desired ratio is $$ Z = 4\frac{s_2^2}{s_1^2 + s_2^2} $$
Notice that this implies that $0 \le Z \le 4$, So a $\chi^2$ approximation will probably not work very well (note also that when $\rho \to 1$ $s_2^2 \to 0$ so $Z$ becomes concentrated at 0, and likewise when $\rho \to -1$ $s_1^2 \to 0$ so $Z$ becomes concentrated at 4).
In fact the distribution of such a ratio of gamma variables is known in closed from and is given in this case by : (See e.g. here)
$$Z/4 \sim f(z)$$
$$f(z) = \frac{z^{k-1}(1-z)^{k-1}}{r^k B(k,k)} \left(1 + \frac{1-r}{r}z \right)^{-2k}, \quad 0 < z < 1.$$
Where $k = (n-1)/2$ , $r = \frac{1-\rho}{1+\rho}$ and $B(\cdotp ,\cdotp)$ is the Beta function.
Furthermore you can easily calculate the correlation between $V_D$ and $V_p$ :
$$\text{Cov}(V_D,V_p) = \text{Cov}( s_2^2 ,\frac{1}{4}(s_1^2 + s_2^2)) = \frac{1}{4} \text{var}(s_2^2) = \frac{2\sigma^4(1-\rho)^2}{n-1}$$
$$\text{var}(V_p) = \frac{1}{16} (\text{var}(s_1^2) + \text{var}(s_2^2))= \frac{\sigma^4(1+\rho^2)}{n-1}$$
$$\rho_{V_D,V_p} = \frac{\text{Cov}(V_D,V_p)}{\sqrt{\text{var}(V_D)\text{var}(V_p)}} = \frac{1-\rho}{\sqrt{2(1+\rho^2)}}$$
(Note that when $\rho=1/2$ you get exactly $1/2\sqrt{2*5/4)} = 1/\sqrt{10}$ )
Best Answer
This post elaborates on the answers in the comments to the question.
Let $X = (X_1, X_2, \ldots, X_n)$. Fix any $\mathbf{e}_1\in\mathbb{R}^n$ of unit length. Such a vector may always be completed to an orthonormal basis $(\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_n)$ (by means of the Gram-Schmidt process, for instance). This change of basis (from the usual one) is orthogonal: it does not change lengths. Thus the distribution of
$$\frac{(\mathbf{e}_1\cdot X)^2}{||X||^2}=\frac{(\mathbf{e}_1\cdot X)^2}{X_1^2 + X_2^2 + \cdots + X_n^2} $$
does not depend on $\mathbf{e}_1$. Taking $\mathbf{e}_1 = (1,0,0,\ldots, 0)$ shows this has the same distribution as
$$\frac{X_1^2}{X_1^2 + X_2^2 + \cdots + X_n^2}.\tag{1} $$
Since the $X_i$ are iid Normal, they may be written as $\sigma$ times iid standard Normal variables $Y_1, \ldots, Y_n$ and their squares are $\sigma^2$ times $\Gamma(1/2)$ distributions. Since the sum of $n-1$ independent $\Gamma(1/2)$ distributions is $\Gamma((n-1)/2)$, we have determined that the distribution of $(1)$ is that of
$$\frac{\sigma^2 U}{\sigma^2 U + \sigma^2 V} = \frac{U}{U+V}$$
where $U = X_1^2/\sigma^2 \sim \Gamma(1/2)$ and $V = (X_2^2 + \cdots + X_n^2)/\sigma^2 \sim \Gamma((n-1)/2)$ are independent. It is well known that this ratio has a Beta$(1/2, (n-1)/2)$ distribution. (Also see the closely related thread at Distribution of $XY$ if $X \sim$ Beta$(1,K-1)$ and $Y \sim$ chi-squared with $2K$ degrees.)
Since $$X_1 + \cdots + X_n = (1,1,\ldots,1)\cdot (X_1, X_2, \cdots, X_n) = \sqrt{n}\,\mathbf{e}_1\cdot X$$
for the unit vector $\mathbf{e}_1=(1,1,\ldots,1)/\sqrt{n}$, we conclude that $Z$ is $(\sqrt{n})^2 = n$ times a Beta$(1/2, (n-1)/2)$ variate. For $n\ge 2$ it therefore has density function
$$f_Z(z) = \frac{n^{1-n/2}}{B\left(\frac{1}{2}, \frac{n-1}{2}\right)} \sqrt{\frac{(n-z)^{n-3}}{z}}$$
on the interval $(0,n)$ (and otherwise is zero).
As a check, I simulated $100,000$ independent realizations of $Z$ for $\sigma=1$ and $n=2,3,10$, plotted their histograms, and superimposed the graph of the corresponding Beta density (in red). The agreements are excellent.
Here is the
Rcode. It carries out the simulation by means of the formulasum(x)^2 / sum(x^2)for $Z$, wherexis a vector of lengthngenerated byrnorm. The rest is just looping (for,apply) and plotting (hist,curve).