If you can assume that the two distributions are roughly normal or lognormal, this would be my approach.
Assume that $\ln(X)$ and $\ln(Y)$ are binormally distributed with means $m,n$, correlation $r$, and standard deviations $s,t$. Then there are formulas $ES(2.5\%, m, s)$ and $ES(2.5\%, n, t)$ for the expected shortfalls.
Meanwhile let $M,N$ be variables for the sample means of $\ln(X)$ and $\ln(Y)$. Let $Means(m,n,r,s,t)$ be the distribution of $M$ and $N$ under the given parameters, which is also binormal. Let $M_0$, $N_0$ be the values of these variables in this dataset.
Similarly let $R,S,T$ be variables for the sample correlation and standard deviations of $\ln(X)$ and $\ln(Y)$. Let $Devs(m,n,r,s,t)$ be the distribution of $R$, $S$ and $T$ under the given parameters, which also has an explicit approximation (recently in papers by Joarder, apparently originally going back to a 1915 paper by Fisher). Let $R_0$, $S_0$ , $T_0$ be the values of these variables in this dataset.
Now let the test statistic be $f(M,N,S,T) = ES(2.5\%, M,S) - ES(2.5\%, N,T)$. We assume that estimates of $M$ and $N$ are roughly independent of estimates of $R$, $S$ and $T$. So a reasonable first null hypothesis is that $f$ is roughly normal with variance the sum of
\begin{align}
Var_{M,N} &=
Var[f(M, N, S_0, T_0)\, |\ M,N \sim Means(M_0, N_0, R_0, S_0, T_0)]\\
Var_{R,S,T}&=Var[f(M_0, N_0, S, T)\, |\, R,S,T \sim Devs(M_0, N_0, R_0, S_0, T_0)]
\end{align}
These variances can be calculated numerically using the formulas above, and the final test is a z-test for $f(M_0, N_0, S_0, T_0)$ under the distribution
$\mathcal{N}(0,Var_{M,N}+Var_{R,S,T})$.
