This topic has piqued my interest on and off for some time now and I'm curious if the methods used here have been discussed somewhere in the literature.
Suppose we have $X\sim\mathcal N(\mu,\sigma)$ and we wish to come up with some meaningful interpretation of the negative integer moments $\mathsf EX^{-n}$, with $n=1,2,\dots$
For $n=1$ we have a means of doing this via the Cauchy principal value. Using this interpretation we find $$ \mathsf EX^{-1}:=\text{p.v.}\int_{\Bbb R}f_X(x)\frac{\mathrm dx}{x}=\frac{\sqrt 2}{\sigma}\mathcal D\left(\frac{\mu}{\sqrt 2\sigma}\right), $$ where $\mathcal{D}(z)=e^{-z^{2}}\int_{0}^{z}e^{t^{2}}\,\mathrm{d}t$ is the Dawson function. What I find particularly interesting is the asymptotic expansion of this expression for $\mathsf EX^{-1}$ gives the same solution as the $\delta$-method when considering $\sigma\nearrow 0$, namely $$ \frac{\sqrt 2}{\sigma}\mathcal D\left(\frac{\mu}{\sqrt 2\sigma}\right)\sim \frac{1}{\mu}+\frac{\sigma^2}{\mu^3}+\mathcal O(\sigma^4). $$ You can verify this is the same solution you get with the $\delta$-method by expanding $g(X)=X^{-1}$ about $\mu$ and then taking the expected value of the first couple terms with respect to the $\mathcal N(\mu,\sigma^2)$ distribution.
What about the case $n=2,3,\dots$? Taking inspiration from Hadamard regularization I defined the auxiliary function $$ M_{X^{-1}}(t):=\text{p.v.}\ \mathsf E(X-t)^{-1}=\frac{\sqrt 2}{\sigma}\mathcal D\left(\frac{\mu-t}{\sqrt 2\sigma}\right). $$ Now make the observation that $$ \left(\frac{1}{(n-1)!}\partial_t^{n-1}(X-t)^{-1}\right)\Bigg|_{t=0}=X^{-n}; $$ hence, $M_{X^{-1}}(t)$ can be viewed as a sort of moment-generating function for the negative moments which leads us to define $$ \mathsf EX^{-n}:=\left(\frac{1}{(n-1)!}\partial_t^{n-1}M_{X^{-1}}(t)\right)\Bigg|_{t=0},\quad n=1,2,\dots $$ As an example, consider the case where $\mu=0$ so that the expression for $\mathsf EX^{-n}$ above gives the "negative central moments". We find $$ \left( \begin{array}{cc} n & (\mathsf EX^{-n})|_{\mu=0}\\ 1 & 0 \\ 2 & -\frac{1}{\sigma ^2} \\ 3 & 0 \\ 4 & \frac{1}{3 \sigma ^4} \\ 5 & 0 \\ 6 & -\frac{1}{15 \sigma ^6} \\ 7 & 0 \\ 8 & \frac{1}{105 \sigma ^8} \\ \end{array} \right) $$ which bear a remarkable resemblance to the positive central moments of the normal distribution. The alternating sign of the even moments concerned me and so I thought this extension did not make sense. However consider the general case $\mu\in\Bbb R$ and define the variance of $X^{-1}$ as $$ \mathsf{Var}X^{-1}=\mathsf EX^{-2}-(\mathsf EX^{-1})^2. $$ I computed $\mathsf{Var}X^{-1}$ in Mathematica using the definition of $\mathsf EX^{-n}$ above and evaluated the asymptotic expansion for $\sigma\nearrow 0$ to find $$ \mathsf{Var}X^{-1}\sim\frac{\sigma^2}{\mu^4}+8\frac{\sigma^4}{\mu^6}+\mathcal O(\sigma^6), $$ which again agrees with the first couple terms you get using the $\delta$-method to "estimate" $\mathsf{Var}X^{-1}$. So it would seem that: $1)$ the moments $\mathsf EX^{-n}$ as defined above are able to assign meaningful values to negative integer moments of the normal distribution and $2)$ the asymptotic expansions of these "moments" agree with what you obtain from the $\delta$-method in the limiting case $\sigma\nearrow 0$.
Have I stumbled into something new (unlikely)? Or has extension of moments for the normal distribution to negative integers already been investigated (more likely but I could not find it)? Also, while I believe I have presented a mildly convincing argument that "my" definition of $\mathsf EX^{-n}$ gives a natural extension of the moments to negative integers it may be that there are situations where this extension does not make sense so I would be interested in any criticism of this approach as well.
Edit:
Here is a simple Mathematica script showing that the asymptotic expansions of these negative moments agree with the $\delta$-method.
M[t_, \[Mu]_, \[Sigma]_] :=
Sqrt[2]/\[Sigma] DawsonF[(\[Mu] - t)/(Sqrt[2] \[Sigma])]
EXn[n_, \[Mu]_, \[Sigma]_] :=
1/(n - 1)! D[M[t, \[Mu], \[Sigma]], {t, n - 1}] /. t -> 0
ord = 9;
EX1[[Mu], [Sigma]] := EXn[1, [Mu], [Sigma]];
EX2[[Mu], [Sigma]] := EXn[2, [Mu], [Sigma]];
VarX1[[Mu], [Sigma]] :=
EX2[[Mu], [Sigma]] - EX1[[Mu], [Sigma]]^2;
Coefficient[Normal[Series[EX1[[Mu], [Sigma]], {[Sigma], 0, ord}]],
Exp[-[Mu]^2/[Sigma]^2], 0];
EX1asym = Expand[Simplify[%, [Sigma] > 0]]
Coefficient[
Normal[Series[VarX1[[Mu], [Sigma]], {[Sigma], 0, ord}]],
Exp[-[Mu]^2/[Sigma]^2], 0];
VarX1asym = Expand[Simplify[%, [Sigma] > 0]]
g = Normal[Series[X^-1, {X, [Mu], ord}]];
EX1delta =
Expand[Expectation[g,
X [Distributed] NormalDistribution[[Mu], [Sigma]]]]
VarX1delta =
Expand[Expectation[g^2,
X [Distributed] NormalDistribution[[Mu], [Sigma]]] -
EX1delta^2]
As you add more terms to these expansions you see that more of the terms agree. So while we get unusual results, i.e. $\mathsf EX^{-2},\,\mathsf{Var}X^{-1}<0$, these generalized negative moments do encode information about the asymptotic properties of $\mathsf EX^{-n}$ when $\mu\neq 0$.