Do the principal values of higher moments exist for the Cauchy distribution? It is known that the Cauchy distribution has undefined moments, and that the expectation has a principal Cauchy value $\operatorname{PV}\left( \mathbb{E} [X] \right)$ of zero.
I wonder if $\operatorname{PV}\left( \mathbb{E} \left[ (X - \operatorname{PV}\left( \mathbb{E} [X] \right))^n \right] \right)$ exists for all $n$, or if there is a highest principal Cauchy moment for the Cauchy distribution.
Do all the higher moment expectations have defined principal Cauchy values for the Cauchy distribution?

The first simplification of the above is to note that $\operatorname{PV} \left(\mathbb{E}[X]\right) = 0 \implies \operatorname{PV}\left( \mathbb{E} \left[ (X - \operatorname{PV}\left( \mathbb{E} [X] \right))^n \right] \right) = \operatorname{PV}\left( \mathbb{E} \left[ X^n \right] \right)$.
From the definition of the Cauchy distribution and the Cauchy principal value we have
$\lim_{\epsilon \rightarrow 0^+} \left[ \int_{-\infty}^{0-\epsilon} x^n \frac{1}{\pi \gamma \left[ 1 + \left(\frac{x-x_0}{\gamma} \right)^2 \right]}dx + \int_{0+\epsilon}^{\infty} x^n \frac{1}{\pi \gamma \left[ 1 + \left(\frac{x-x_0}{\gamma} \right)^2 \right]}dx \right].$
So my question reduces to whether the above expression is defined.

Wolfram Alpha tells me that
$\int x^n \frac{1}{\pi \gamma \left[ 1 + \left(\frac{x-x_0}{\gamma} \right)^2 \right]}dx = \frac{\gamma x^{n+1}}{2 \pi \sqrt{- \gamma^2} (n+1)} \left( \frac{_2F_1 \left(1, n+1;n+2;\frac{x}{m - \sqrt{- \gamma^2}} \right)}{m-\sqrt{- \gamma^2}} - \frac{_2F_1 \left(1, n+1;n+2;\frac{x}{m + \sqrt{- \gamma^2}} \right)}{m + \sqrt{- \gamma^2}} \right) + C$
where $_2F_1 (\cdot, \cdot ; \cdot ; \cdot )$ is a hypergeometric function.

The odd integral moments have principal values of zero while the even integral moments have infinite principal values.

By definition, the Cauchy Principal Value [w.r.t. infinity] of an integral over the real line is
$$\operatorname{pv}\int_\mathbb{R} g(x)\,\mathrm{d}x = \lim_{R\to\infty} \int_{-R}^R g(x)\,\mathrm{dx}.$$
When $g$ is antisymmetric -- that is, $g(x) = -g(x)$ for all $x,$ obviously the right hand side is always zero, whence its limit is zero.
The Cauchy density is proportional to $1/(1+x^2),$ which is symmetric about $0.$  Consequently, for any antisymmetric integrable function $h,$
$$g(x)  = \frac{h(x)}{1+x^2}$$
is also integrable and antisymmetric, whence
$$\operatorname{pv} \int_\mathbb R \frac{h(x)}{1+x^2}\,\mathrm{d}x = 0.$$
This is the case for $h(x) = x^{2k+1}$ when $k$ is a natural number.  It is not the case when $h(x) = x^{2k}.$  With $k=0,$ the integral converges, but for $k\ge 1$ use the simple fact that when $|x|\ge 1,$
$$\frac{x^{2k}}{1+x^2} \ge \frac{x^{2k}}{x^2+x^2} = \frac{1}{2}x^{2k-2} \ge \frac{1}{2}$$ to estimate
$$\begin{aligned}
\lim_{R\to\infty}\int_{-R}^R \frac{x^{2k}}{1+x^2}\,\mathrm{d}x &\ge \lim_{R\to\infty}\left(\int_{-R}^{-1} + \int_1^R\right) \frac{x^{2k}}{1+x^2}\,\mathrm{d}x\\
&\ge \lim_{R\to\infty}\left(\int_{-R}^{-1} + \int_1^R\right) \frac{1}{2}\,\mathrm{d}x\\
&= \lim_{R\to\infty} R - 1 =\infty.
\end{aligned}$$

A similar argument by symmetry handles the case of $negative$ integral $n,$ but now the principal value is computed by excising a small neighborhood $(-\epsilon,\epsilon)$ around zero and taking the limit (from above) as $\epsilon\to 0.$  Again there's cancellation for odd $n$ but divergence for even $n.$  Alternatively, change variables from $x$ to $1/x,$ which reduces the integral to the form explicitly considered here.
