# How to compute confidence intervals for a scaled reciprocal normal statistic

I have the following statistic of which I want to compute confidence intervals: $$Z = \frac{(X-Y)^{-1}}{\mathbb{PVE}[(X-Y)^{-1}]},$$ where \begin{align} X &\sim\mathcal{N}\left(\sigma_{1}^{2},\frac{2\sigma_{1}^{4}}{n-1}\right),\\ Y &\sim\mathcal{N}\left(\sigma_{2}^{2},\frac{2\sigma_{2}^{4}}{m-1}\right), \end{align} $n$ and $m$ are the sample sizes used to compute $X$ and $Y$ and $\mathbb{PVE}[(X-Y)^{-1}]$ is the Cauchy principal value for the expected value of $(X-Y)^{-1}$ (it's just a constant) defined as $$\mathbb{PVE}[(X-Y)^{-1}] = \frac{\sqrt{2}}{\sigma_{X-Y}}\,\mathcal{D}_{+}\left(\frac{\mu_{X-Y}}{\sqrt{2}\sigma_{X-Y}}\right),$$ where $\mathcal{D}_{+}(z)$ is the dawson function. I have parametric equations for the pdf/cdf of $Z$. Computing exact confidence intervals seems completely out of the question.

For my application, very large sample sizes are the norm (i.e. $n$ and $m$ on the order of several thousands). It is also typical for $\mathrm{Pr}[(X-Y)\leq0]\approx 0$. How can I go about computing confidence intervals of $Z$? What options do I have in this specific application?