Expected value of the projective metric between random orthogonal Stiefel matrices in $\mathbb{R}^{N \times k}$ equals $1 - \frac{k}{N}$

Question

This is a call for help to the random-matrix-theory savvy people. I've observed the below equality in experiments, and have been looking for a proof in the RMT literature but couldn't find one. I'd appreciate some pointers!

---

Definitions:

Consider the set of orthogonal matrices in the compact Stiefel manifold:

$\mathcal{M}_{N, k} := \{o \in \mathbb{R}^{N \times k} : \quad o^\top o = I_k \}$

We also define the projector $P_i$ of a matrix $o_i \in \mathcal{M}_{N, k}$ as:

$P_i = o_i o_i^\top$

And the principal angles $\theta_{i\rightarrow j}$ between 2 matrices $o_i, o_j$ (see e.g. Absil, Mahony, Sepulchre, 2003 or Qiu, Zhang, Li, 2005, s.2) as:

$SVD[o_i^\top o_j] = U \quad cos(\theta_{i\rightarrow j}) \quad V^T$

Finally, consider the following (squared and normalized) Frobenius projective metric between any 2 matrices $(o_i, o_j)$ in $\mathcal{M}_{N, k}$ (see e.g. Edelman, Arias, Smith, 1998, section 4.3):

$d_{pF}^2(o_i, o_j) = \frac{1}{k} \lVert P_i - P_j \rVert_F^2 = \frac{1}{k} \lVert sin^2(\theta_{i\rightarrow j}) \rVert_F^2 = 1 - \frac{1}{k} \lVert cos ^2(\theta_{i\rightarrow j}) \rVert_F^2 = 1 - \frac{1}{k} \lVert o_i^\top o_j \rVert_F^2 \in [0, 1]$

---

Setup:

1. Draw a random pair $(o_1, o_2) \sim (\mathcal{M}_{N, k} \times \mathcal{M}_{N, k})$ uniformly from the Haar measure (e.g. using the QR method)
2. Compute $d_{pF}^2(o_1, o_2)$
3. Repeat steps 1 and 2 for different random pairs, and average the results. We observe that $\frac{1}{k} \lVert o_i^\top o_j \rVert_F^2$ converges to $\frac{k}{N}$. Or equivalently, $\mathbb{E}[d_{pF}^2] = 1 - \frac{k}{N}$

---

Question:

I imagine this is a well known result in the literature. Where could I find it? Alternatively, how could I prove it?

The approach I thought of so far seems rather cumbersome, and I am not sure if it would be correct:

1. Take the joint distribution $\mathcal{X}$ of principal angles for a Haar-uniform random orthogonal matrix in $\mathbb{R}^{N \times N}$ (presented e.g. in Rummler, 2002).
2. Truncate $\mathcal{X}$ to $k$ entries, and calculate the normalized squared cosine, i.e. $\mathcal{Y} = \frac{1}{k} cos^2(\mathcal{X}_k)$
3. Then, $\mathbb{E}[\mathcal{Y}]$ should equal $\frac{k}{N}$ (?)

Again, this is probably not necessary since I assume this result is known and I'm just missing the reference. Thanks in advance!