Context
I am interested in the Cayley transform on the Stiefel manifold $\mathcal V_{np}$ at the point $I_{np} = (I_p, 0_{p,n-p})^\top$. For a matrix $X$ of the form $$ X = \begin{pmatrix} A & -B^\top\\ B & 0 \end{pmatrix} $$ with $A\in \mathbb R^{p\times p}$ such that $A^\top = -A$ and $B\in \mathbb R^{p\times(n-p)}$, the Cayley transform $C: \mathbb R^{n\times n}\to \mathcal V_{np}$ is defined by $$ C(X) = (I - X/2)^{-1}(1 + X/2) I_{np} $$ In practice, $X$ is parameterized linearly by a vector $\varphi$ with dimension $d_\mathcal V=p(p-1)/2 + p(n-p)$, and we are actually interested in the map $\bar C(\varphi) = C(X_\varphi)$.
My question
In order to make changes of variable relying on $\bar C$, a recent paper [1] uses the generalized Jacobian determinant: $$ J(\varphi) = | d_\varphi \bar C^* d_\varphi \bar C | \,. $$ I would like to prove that $J(-\varphi) = J(\varphi)$. I can easily verify this property numerically by approximating the Jacobian with finite differences. It is also relatively intuitive if we think of the Stiefel manifold as a generalization of the sphere.
What I have so far
The paper mentioned above [1] also showed in particular that $$ J(\varphi) = |\Gamma_\mathcal V^\top \left[\left[ (I_n - X_\varphi)^{-1} I_{np}I_{np}^\top (I_n - X_\varphi)^{-\top} \right] \otimes \left[ (I_n - X_\varphi)^{-\top} (I_n - X_\varphi)^{-1} \right]\right]\Gamma_\mathcal V| $$ (with $\Gamma_\mathcal V$ denoting the linear map $\varphi \mapsto X_\varphi$ ; it is thus a rectangular $n^2 \times d_\mathcal V$ matrix).
Using the property of the Kronecker product, this gives in particular $$ J(\varphi) = |\Gamma_\mathcal V^\top \left[ (I_n - X_\varphi)^{-1} \otimes (I_n - X_\varphi)^{-\top}\right] \left[ I_{np}I_{np}^\top \otimes I_n\right] \left[ (I_n - X_\varphi)^{-\top} \otimes (I_n - X_\varphi)^{-1} \right]\Gamma_\mathcal V|\,. $$ Using the fact that $X_\varphi$ is skew-symmetric, and using commutation matrices, I can show that $$ J(-\varphi) = |\Gamma_\mathcal V^\top \left[ (I_n - X_\varphi)^{-1} \otimes (I_n - X_\varphi)^{-\top}\right] \left[ I_n \otimes I_{np}I_{np}^\top \right] \left[ (I_n - X_\varphi)^{-\top} \otimes (I_n - X_\varphi)^{-1} \right]\Gamma_\mathcal V|\,, $$ i.e. only the Kronecker product in the middle differs between the matrices in $J(\varphi)$ and $J(-\varphi)$. Numerically, the two formulas above match even when replacing $X_\varphi$ by a random skew-symmetric matrix, which leads me to hope that relatively simple linear algebra may allow to conclude.
Any help would be appreciated!
[1] Jauch, M., Hoff, P. D., & Dunson, D. B. (2020). Random orthogonal matrices and the Cayley transform. Bernoulli, 26(2), 1560–1586. https://doi.org/10.3150/19-BEJ1176
Differential writes \begin{eqnarray} d\mathbf{C} &=& +\frac12 \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} (d\mathbf{X}) \left(\mathbf{C}+ \mathbf{I}_{np} \right) \\ &=& \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} (d\mathbf{X}) \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} \mathbf{I}_{np} \end{eqnarray} since \begin{eqnarray*} \mathbf{C}(\mathbf{X})+\mathbf{I}_{np} &=& \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} \left \lbrace \left( \mathbf{I} + \mathbf{X}/2 \right) + \left( \mathbf{I} - \mathbf{X}/2 \right) \right \rbrace \mathbf{I}_{np} \\ &=& 2 \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} \mathbf{I}_{np} \end{eqnarray*}
From here use the vec trick \begin{eqnarray} d\mathbf{c} &=& \left \lbrace \mathbf{I}_{np}^T \left( \mathbf{I} - \mathbf{X}/2 \right)^{-T} \otimes \left( \mathbf{I} - \mathbf{X}/2 \right)^{-1} \right \rbrace d\mathbf{x} \end{eqnarray} The Jacobian writes $\mathbf{J}=\mathbf{U}\otimes \mathbf{V}$ where $\mathbf{U}=\mathbf{I}_{np}^T \left( \mathbf{I} - \mathbf{X}/2 \right)^{-T}$ and $\mathbf{V}=\left( \mathbf{I} - \mathbf{X}/2 \right)^{-1}$.
Looking at the Jacobian determinant, we must compute the determinant of the matrix $$ \mathbf{J}^T \otimes \mathbf{J} = (\mathbf{U}\otimes \mathbf{V})^T (\mathbf{U}\otimes \mathbf{V}) = \mathbf{U}^T \mathbf{U} \otimes \mathbf{V}^T \mathbf{V} $$
It is easy to see that $\mathbf{V}^T \mathbf{V}$ evaluated for $-\varphi$ is equal to $\mathbf{V}^T \mathbf{V}$ evaluated for $\varphi$. and similar for $\mathbf{U}^T \mathbf{U}$. This shows that the Jacobian determinant is unchanged if $\varphi \rightarrow \varphi$.
To do that, we can use the fact that $\mathbf{X}^T=-\mathbf{X}$ so that $\mathbf{V}(-\varphi) = \left[ \mathbf{I} - \mathbf{X}(-\varphi)/2 \right]^{-1} = \left( \mathbf{I} + \mathbf{X}(\varphi)/2 \right]^{-1} = \left[ \mathbf{I} - \mathbf{X}(\varphi)/2 \right]^{-T} = \mathbf{V}^T $.