As @Maxim mentioned, the Laplas method can be applied to the integral, and the steepest descend is along the real axis, near $t=1$. The asymptotics can be obtained as the asymptotics of special function, because $I(x)=2K_{-n-1}(2x)$.
We can also get the asymptotics (few several terms) directly.
Let
$$I_n(x)=\int_{0}^{\infty}t^{n}e^{-x(t+\frac{1}{t})}dt$$
We can consider the generating function instead, namely
$$J(\alpha,x)=\sum_{n=0}^\infty\frac{\alpha ^n}{n!}I_n(x)=\int_{0}^{\infty}e^{\alpha t}e^{-x(t+\frac{1}{t})}dt=\int_{0}^{\infty}e^{\alpha t}e^{-xf(t)}dt$$
where $f(t)=\frac{1}{t}+t$.
Then $\,\,I_n(x)=\frac{\partial^n }{\partial\alpha^n}J(\alpha,x)\Big|_{\alpha=0}$
$f'(t)=0$ at $t=1$; therefore near this point
$$f(t)=f(1)+f'(1)(t-1)+f''(1)\frac{(t-1)^2}{2!}+f'''(1)\frac{(t-1)^3}{3!}+...$$
$$f(t)=2+(t-1)^2-(t-1)^3+(t-1)^4-=...$$
$$J(\alpha,x)=e^{-2x}\int_{0}^{\infty}e^{\alpha t}e^{-x\big((t-1)^2-(t-1)^3+(t-1)^4-+...\big)}dt=e^{-2x}\int_{-1}^{\infty}e^{\alpha (t+1)}e^{-x\big(t^2-t^3+t^4-+...\big)}dt$$
$$=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\sqrt x}^{\infty}e^{\frac{\alpha t}{\sqrt x}}e^{-t^2}e^{\frac{t^3}{\sqrt x}-\frac{t^4}{x}+ -...}dt=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{\infty}\,-\,\,\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{-\sqrt x}$$
The second integral contains an additional exponentially small factor ($\sim e^{-x}$), so we drop it and focus on the first term. Expanding the integrand near $t=0$ into the series
$$J(\alpha,x)=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{\infty}\Big(1+\frac{\alpha t}{\sqrt x}+\frac{\alpha^2 t^2}{2x}+...\Big)\Big(1+\frac{t^3}{\sqrt x}+\frac{t^6}{2x}-\frac{t^4}{x}+O\big(x^{-\frac{3}{2}}\big)\Big)e^{-t^2}dt$$
Using $\int_{-\infty}^{\infty}e^{-t^2}t^{2k}dt=(-1)^k\frac{\partial^k}{\partial s^k}\int_{-\infty}^{\infty}e^{-st^2}dt\Big|_{s=1}=(-1)^k\sqrt\pi\frac{\partial^k}{\partial s^k}\frac{1}{\sqrt s}\Big|_{s=1}$, integrating term by term and grouping the terms with the same power of $x$, we get
$$J(\alpha,x)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\,e^{\alpha}\Big(1+\frac{\alpha^2+3\alpha+\frac{3}{4}}{4x}+O\big(\frac{1}{x^2}\big)\Big)$$
$$I_n(x)=\frac{\partial^n }{\partial\alpha^n}J(\alpha,x)\Big|_{\alpha=0}$$
$$I(0)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{3}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$
$$I(1)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{15}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$
$$I(2)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{35}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$
$$....$$
$$I_n(x)=\sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{(2n+3)(2n+1)}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$
The main asymptotic term of $I_n(x)$ does not depend on $n$.
Best Answer
In your example, the only critical point of $p$ is $z = 0$. $\operatorname{Re} i p$ is constant on the lines $u v = 0$ and $v = \pm u$, dividing the plane into eight octants, and $\operatorname{Im} i p$ is constant on the bisectors of the octants. The goal is to prove that the integral is asymptotically equivalent to the integral over a small part of a bisector around $z = 0$, then any other critical points will be irrelevant.
For that, you need to choose the bisectors of those angles inside which the contour can be deformed to the real line while $\operatorname{Re} i p$ remains less than its value at $z = 0$. There is only one pair of angles that fits the requirement.