I'm not sure about a reference for this result, but it is possible to derive it relatively easily, so I hope that suffices. One way to approach this problem is to look at it as a problem involving a quadratic form taken on a normal random vector. The pooled sample variance can be expressed as a quadratic form of this kind, and these quadratic forms are generally approximated using the chi-squared distribution (with exact correspondence in some cases).
Derivation of the result: In order to show where your assumptions come into the derivation, I will do the first part of the derivation without assuming equal variances for the two groups. If we denote your vectors by $\mathbf{X} = (X_1,...,X_n)$ and $\mathbf{Y} = (Y_1,...,Y_n)$ then your stipulated problem gives the joint normal distribution:
$$\begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \end{bmatrix}
\sim \text{N} (\boldsymbol{\mu}, \mathbf{\Sigma} )
\quad \quad \quad
\boldsymbol{\mu} =
\begin{bmatrix} \mu_X \mathbf{1} \\ \mu_Y \mathbf{1} \end{bmatrix}
\quad \quad \quad
\mathbf{\Sigma} =
\begin{bmatrix} \sigma_X^2 \mathbf{I} & \rho \sigma_X \sigma_Y \mathbf{I} \\
\rho \sigma_X \sigma_Y \mathbf{I} & \sigma_Y^2 \mathbf{I} \end{bmatrix}.$$
Letting $\mathbf{C}$ denote the $n \times n$ centering matrix, you can write the pooled sample variance in this problem as the quadratic form:
$$\begin{align}
S_\text{pooled}^2
&= \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \end{bmatrix}^\text{T}
\mathbf{A} \begin{bmatrix} \mathbf{X} \\ \mathbf{Y} \end{bmatrix}
\quad \quad \quad
\mathbf{A} \equiv \frac{1}{2(n-1)} \begin{bmatrix} \mathbf{C} & \mathbf{0} \\ \mathbf{0} & \mathbf{C} \end{bmatrix}. \\[6pt]
\end{align}$$
Now, using standard formulae for the mean and variance of quadradic forms of normal random vectors, and noting that $\mathbf{C}$ is an idempotent matrix (i.e., $\mathbf{C} = \mathbf{C}^2$), you have:
$$\begin{align}
\mathbb{E}(S_\text{pooled}^2)
&= \text{tr}(\mathbf{A} \mathbf{\Sigma}) + \boldsymbol{\mu}^\text{T} \mathbf{A} \boldsymbol{\mu} \\[6pt]
&= \text{tr} \Bigg( \frac{1}{2(n-1)} \begin{bmatrix} \sigma_X^2 \mathbf{C} & \rho \sigma_X \sigma_Y \mathbf{C} \\ \rho \sigma_X \sigma_Y \mathbf{C} & \sigma_Y^2 \mathbf{C} \end{bmatrix} \Bigg) + \mathbf{0} \\[6pt]
&= \frac{1}{2(n-1)} \text{tr} \Bigg( \begin{bmatrix} \sigma_X^2 \mathbf{C} & \rho \sigma_X \sigma_Y \mathbf{C} \\ \rho \sigma_X \sigma_Y \mathbf{C} & \sigma_Y^2 \mathbf{C} \end{bmatrix} \Bigg) \\[6pt]
&= \frac{1}{2(n-1)} \Bigg[ n \times \frac{n-1}{n} \cdot \sigma_X^2 + n \times \frac{n-1}{n} \cdot \sigma_Y^2 \Bigg] \\[6pt]
&= \frac{\sigma_X^2 + \sigma_Y^2}{2}, \\[12pt]
\mathbb{V}(S_\text{pooled}^2)
&= 2 \text{tr}(\mathbf{A} \mathbf{\Sigma} \mathbf{A} \mathbf{\Sigma}) + 4 \boldsymbol{\mu}^\text{T} \mathbf{A} \mathbf{\Sigma} \mathbf{A} \boldsymbol{\mu} \\[6pt]
&= 2 \text{tr} \Bigg( \frac{1}{4(n-1)^2} \begin{bmatrix} \sigma_X^2 \mathbf{C} & \rho \sigma_X \sigma_Y \mathbf{C} \\ \rho \sigma_X \sigma_Y \mathbf{C} & \sigma_Y^2 \mathbf{C} \end{bmatrix}^2 \Bigg) + \mathbf{0} \\[6pt]
&= \frac{1}{2(n-1)^2} \text{tr} \Bigg( \begin{bmatrix} (\sigma_X^4 + \rho^2 \sigma_X^2 \sigma_Y^2) \mathbf{C} & (\sigma_X^2 + \sigma_Y^2) \rho \sigma_X \sigma_Y \mathbf{C} \\ (\sigma_X^2 + \sigma_Y^2) \rho \sigma_X \sigma_Y \mathbf{C} & (\sigma_Y^4 + \rho^2 \sigma_X^2 \sigma_Y^2) \mathbf{C} \end{bmatrix} \Bigg) \\[6pt]
&= \frac{1}{2(n-1)^2} \Bigg[ n \times \frac{n-1}{n} \cdot (\sigma_X^4 + \rho^2 \sigma_X^2 \sigma_Y^2) + n \times \frac{n-1}{n} \cdot (\sigma_Y^4 + \rho^2 \sigma_X^2 \sigma_Y^2) \Bigg] \\[6pt]
&= \frac{1}{2(n-1)} \Bigg[ (\sigma_X^4 + \rho^2 \sigma_X^2 \sigma_Y^2) + (\sigma_Y^4 + \rho^2 \sigma_X^2 \sigma_Y^2) \Bigg] \\[6pt]
&= \frac{\sigma_X^4 + \sigma_Y^4 + 2 \rho^2 \sigma_X^2 \sigma_Y^2}{2(n-1)}. \\[12pt]
\end{align}$$
Using the equal variance assumption we have $\sigma_X = \sigma_Y = \sigma$ so the moments reduce to:
$$\mathbb{E} \bigg( \frac{S_\text{pooled}^2}{\sigma^2} \bigg) = 1
\quad \quad \quad
\mathbb{V} \bigg( \frac{S_\text{pooled}^2}{\sigma^2} \bigg) = \frac{1+\rho^2}{n-1}.$$
It is usual to approximate the distribution of the quadratic form by a scaled chi-squared distribution using the method of moments. Equating the first two moments to that distribution gives the variance requirement $\mathbb{V}(S_\text{pooled}^2/\sigma^2) = 2/\nu$, which then gives the degrees-of-freedom parameter:
$$\nu = \frac{2(n-1)}{1+\rho^2}.$$
Bear in mind that the degrees-of-freedom parameter here depends on the true correlation coefficient $\rho$, and you may need to estimate this using the sample correlation in your actual problem.
