Consider following integral
$$\mathcal{I}(x) = \int_{\mathbb{R}^n} \frac{ e^{ix \cdot \xi}}{1 + |\xi|^2}\;d^n\xi$$
Strictly speaking, for $n > 1$, such integral doesn't exist in Lebesgue sense because the integrand falls out too slowly for large $|x|$. When the integrand is not absolutely integrable, the usual way to define the inverse Fourier transform is consider following
limit:
$$
\mathcal{I}_{alt}(x) = \lim_{\Lambda\to\infty}\int_{|\xi|\le\Lambda}\frac{ e^{ix \cdot \xi}}{1 + |\xi|^2}\;d^n\xi
$$
If I didn't make any mistake, this regularization using a cutoff $\Lambda$ works fine for $n \le 4$ but fails when $n \ge 6$ (I have no idea for $n = 5$). Furthermore, it is hard to deduce the exact expression of $\mathcal{I}_{alt}(x)$ this way.
For practical application, we need to find an alternate regularization for the integral.
What regularized version to use depends on application. If rotation symmetry is desirable,
we can compute following integral
$$\mathcal{I}_\epsilon(x) = \int_{\mathbb{R}^n} \frac{ e^{ix \cdot \xi}}{1 + |\xi|^2} e^{-\epsilon(1+|\xi|^2)}\;d^n\xi$$
and take the limit $\epsilon \to 0_{+}$ instead.
I will skip further technical justification and just derive the final result.
Using the integral representation
$$\frac{1}{1+|\xi|^2} e^{-\epsilon(1+|\xi|^2)} = \int_\epsilon^\infty e^{-\lambda(1+|\xi|^2)} d\lambda$$
We can rewrite $\mathcal{I}_\epsilon(x)$ as
$$\begin{align}
\mathcal{I}_\epsilon(x)
&= \int_\epsilon^\infty \int_{\mathbb{R}^n} e^{-\lambda-\lambda |\xi|^2 + ix\cdot\xi} d^n\xi d\lambda
= \int_\epsilon^\infty e^{-\lambda - \frac{|x|^2}{4\lambda}}\left(\int_{\mathbb{R}^n} e^{-\lambda \left|\xi - i\frac{x}{2\lambda} \right|^2} d^n\xi \right) d\lambda\\
&= \int_\epsilon^\infty e^{-\lambda - \frac{|x|^2}{4\lambda}}\left(\frac{\pi}{\lambda}\right)^{\frac{n}{2}} d\lambda
\end{align}
$$
Assume $|x| \ne 0$, introduce variable $t$ such that $\displaystyle\;\lambda = \frac{|x|}{2} e^{t}$, above expression becomes
$$\mathcal{I}_\epsilon(x) = \pi\left(\frac{2\pi}{|x|}\right)^{\frac{n}{2}-1} \int_{\log\frac{2\epsilon}{|x|}}^\infty e^{ -|x|\cosh(t)} e^{-(\frac{n}{2}-1)t} dt\tag{*1}$$
Compare this with the integral representation of the
modified Bessel function of the second kind.
$$K_\alpha(z) = \int_0^\infty e^{-z\cosh(t)}\cosh(\alpha t) dt$$
and notice as $\epsilon \to 0_{+}$, $\log\frac{2\epsilon}{|x|} \to -\infty$
and the integral in $(*1)$ is well behaved in this limit, we obtain
$$
\lim_{\epsilon\to 0_{+}}\mathcal{I}_\epsilon(x) = 2\pi\left(\frac{2\pi}{|x|}\right)^{\frac{n}{2}-1} K_{\frac{n}{2}-1}(|x|)\tag{*2}
$$
As a double check, consider the case $n = 1$. Since $\displaystyle\;K_{-\frac12}(z) = \sqrt{\frac{\pi}{2z}}e^{-z}$, $(*2)$ reduces to
$$\int_{-\infty}^{\infty} \frac{e^{ix\xi}}{1+\xi^2} d\xi = 2\pi\frac{1}{\sqrt{\frac{2\pi}{|x|}}}\sqrt{\frac{\pi}{2|x|}}e^{-|x|} = \pi e^{-|x|}$$
reproducing what we know for the $1$-dimensional case.