Here is my trial, which is partially successful but still not fully answering to your question.
Using some coordinate change, I derived that
$$ I_{n} := \int_{0}^{\infty} \mathrm{erfc}^{n}(x) \, dx = \frac{1}{\sqrt{n}} \left( \frac{2}{\sqrt{\pi}} \right)^{n-1} \int_{T^{n-1}} \int_{0}^{\infty} s^{n-1}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds d\sigma_{x}, \tag{*} $$
where $T^{n-1} = \{ x \in [0, \infty)^{n} : x_{1} + \cdots + x_{n} = \sqrt{n} \}$ is an $(n-1)$-simplex and $d\sigma_{x}$ is the surface measure on $T^{n-1}$.
Case $n = 1$. Since $T^{0} = \{1\}$ and $d\sigma_{x} = \delta(x-1)$ is a unit mass located at $x = 1$, the equation $\text{(*)}$ gives no new information.
Case $n = 2.$ It is easy to find that
$$ \int_{0}^{\infty} s e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2|x|(|x|+1)}$$
and by suitable parametrization of the line segment $T^{1}$, we get
\begin{align*}
I_{2}
&= \sqrt{\frac{2}{\pi}} \int_{T^{1}} \frac{1}{2|x|(|x|+1)} \, d\sigma_{x} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\sqrt{2}} \frac{\sqrt{2} \, dt}{2\sqrt{t^{2} + (\sqrt{2}-t)^{2}} ( \sqrt{t^{2} + (\sqrt{2}-t)^{2}} + 1)} \\
&= \sqrt{\frac{2}{\pi}} \int_{0}^{\frac{\pi}{4}} \frac{d\theta}{1+\cos\theta}
= \sqrt{\frac{2}{\pi}} \tan\left(\frac{\pi}{8}\right)
= \frac{2 - \sqrt{2}}{\sqrt{\pi}}.
\end{align*}
We can also write
$$ I_{2} = \boxed{ \displaystyle \frac{2}{\sqrt{\pi}} - \frac{2\sqrt{2}}{\pi^{3/2}} \arctan(\infty) }. $$
Case $n = 3$. We have
$$ \int_{0}^{\infty} s^{2}e^{-(|x|^{2}-1)s^{2}} \mathrm{erfc}(s) \, ds = \frac{1}{2\sqrt{\pi}(|x|^{2}-1)} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right). $$
Now we know that $T^{2}$ is a regular triangle with side length $\sqrt{6}$. Thus by introducing polar coordinate change, we get
\begin{align*}
I_{3}
&= \frac{2}{\pi^{3/2}\sqrt{3}} \int_{T^{2}} \frac{1}{|x|^{2}-1} \left( \frac{\arctan\sqrt{|x|^{2}-1}}{\sqrt{|x|^{2}-1}} - \frac{1}{|x|^{2}} \right) \, d\sigma_{x} \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \int_{0}^{\frac{\sec\theta}{\sqrt{2}}} \frac{1}{r} \left( \frac{\arctan r}{r} - \frac{1}{r^{2} + 1} \right) \, drd\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{0}^{\frac{\pi}{3}} \left( 1- \frac{\arctan\left(\sec\theta / \sqrt{2}\right)}{\sec\theta / \sqrt{2}} \right) \, d\theta \\
&= \frac{4\sqrt{3}}{\pi^{3/2}} \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}}.
\end{align*}
I haven't tried the last integral, but Mathematica suggests that
$$ \int_{\frac{1}{\sqrt{2}}}^{\sqrt{2}} \left( 1 - \frac{\arctan u}{u} \right) \frac{du}{u\sqrt{2u^{2} - 1}} = \frac{\sqrt{3}\pi}{4} -\sqrt{\frac{3}{2}} \arctan\left( \sqrt{2} \right).$$
This gives
$$ I_{3} = \boxed{ \displaystyle \frac{3}{\sqrt{\pi}} - \frac{6\sqrt{2}}{\pi^{3/2}} \arctan\left( \sqrt{2} \right) }, $$
which can be numerically checked.
I am trying to generalize this calculation for higher dimensions, but it seems somewhat daunting. For example, if $n = 4$ we have to evaluate
$$ \int_{0}^{\infty} \mathrm{erfc}^{4}(x) \, dx = \frac{1}{\pi^{3/2}} \int_{T^{3}} \frac{2|x|+1}{(|x|+1)^{2}|x|^{3}} \, d\sigma_{x}. $$
Nevertheless, using the formula
$$ \int_{T^{n-1}} f(|x|) \, d\sigma_{x} = n! \int_{0}^{\sqrt{\frac{1}{n-1}}} \int_{0}^{\sqrt{\frac{n}{n-2}}s_{1}} \cdots \int_{0}^{\sqrt{\frac{3}{1}}s_{n-2}} f\left( \sqrt{1+s_{1}^{2} + \cdots + s_{n-1}^{2}} \right) \, ds_{n-1} \cdots ds_{1} $$
together with aid of Mathematica, I was able to evaluate $I_{4}$ and it was
$$ I_{4} = \boxed{ \displaystyle \frac{4}{\sqrt{\pi}} - \frac{24\sqrt{2}}{\pi^{3/2}} \arctan \left( \frac{1}{\sqrt{8}} \right) }. $$
These series of observations lead us to believe that $I_{n}$ is of the form
$$ I_{n} = \frac{n}{\sqrt{\pi}} - \frac{n!\sqrt{2}}{\pi^{3/2}} \arctan \alpha_{n} $$
for some reasonably nice $\alpha_{n}$ (with $\alpha_{1} = 0$, $\alpha_{2} = \infty$, $\alpha_{3} = \sqrt{2}$ and $\alpha_{4} = \frac{1}{\sqrt{8}}$), but inverse symbolic calculations show that this seems no longer the case for $n \geq 5$.
Indeed, for $n = 5$ we have
$$ I_{5} = \frac{5}{\sqrt{\pi}} - \frac{80\sqrt{2}}{\pi^{3/2}} \arctan \sqrt{\frac{2}{3}} + \frac{240\sqrt{2}}{\pi^{5/2}} A\left( \sqrt{\frac{5}{3}}, \sqrt{\frac{3}{2}}, \sqrt{\frac{4}{5}} \right), $$
where $A(p, q, r)$ is the generalized Ahmed's integral. I have no idea whether it will reduce to a familiar closed form or not.
By the substitution $x\mapsto x^{-1}$ it is not hard to see that
$$I(\nu^{-1},1)=\int^\infty_0\frac{\arctan^2{x}\arctan(\nu x)}{x^2}{\rm d}x$$
First, start off by re-expressing the integral
\begin{align}
\int^\pi_0x^2\cos(nx)\cot\left(\frac{x}{2}\right)\ {\rm d}x
&=-{\rm Re}\int_C\frac{z+1}{z-1}z^{n-1}\ln^2{z}\ {\rm d}z\\
&=(-1)^{n}\int^1_0\frac{1-x}{1+x}x^{n-1}(\ln^2{x}-\pi^2)\ {\rm d}x-\int^1_0\frac{1+x}{1-x}x^{n-1}\ln^2{x}\ {\rm d}x\\
\end{align}
where $C$ is the arc joining $z=1$ to $z=-1$. (The first equality follows from $z=e^{ix}$, whereas the second follows from the fact that the integral along $[-1,\epsilon]\cup\epsilon\exp(i[\pi,0])\cup[\epsilon,1]\cup C$ is $0$ since $z=1$ is a removable singularity and the indent around $z=0$ vanishes.)
Next, note that
\begin{align}
I_\nu(\nu^{-1},1)
&=\int^\frac{\pi}{2}_0\frac{x^2}{\tan{x}(\cos^2{x}+\nu^2\sin^2{x})}{\rm d}x\\
&=\frac{1}{4}\int^\pi_0\frac{x^2}{\tan\left(\frac{x}{2}\right)(1+(1-\nu^2)\cos{x}+\nu^2)}{\rm d}x\\
&=\frac{1}{8\nu}\int^\pi_0\left(x^2\cot\left(\frac{x}{2}\right)+2\sum^\infty_{n=1}\left(\frac{\nu-1}{\nu+1}\right)^nx^2\cos(nx)\cot\left(\frac{x}{2}\right)\right){\rm d}x\\
&=\frac{\pi^2}{4\nu}\ln{2}-\frac{7\zeta(3)}{8\nu}-\frac{\xi}{4\nu}\left(\int^1_0\frac{1-x}{(1+x)(1+\xi x)}(\ln^2{x}-\pi^2)+\frac{1+x}{(1-x)(1-\xi x)}\ln^2{x}\ {\rm d}x\right)
\end{align}
Here $\xi=\dfrac{\nu-1}{\nu+1}$. Utilising the partial fraction decompositions
$$\frac{1-x}{(1+x)(1+\xi x)}=\frac{\nu+1}{1+x}-\frac{\nu}{1+\xi x}$$
$$\frac{1+x}{(1-x)(1-\xi x)}=\frac{\nu+1}{1-x}-\frac{\nu}{1-\xi x}$$
in tandem with the easily verifiable fact
$$\int^1_0\frac{\ln^2{x}}{1+\lambda x}{\rm d}x=-\frac{2{\rm Li}_3(-\lambda)}{\lambda}$$
yields
\begin{align}
I_\nu(\nu^{-1},1)
&=\frac{\pi^2}{4\nu}\ln{2}-\frac{7\zeta(3)}{8\nu}-\frac{\xi}{4\nu}\left(\frac{\pi^2\nu}{\xi}\ln(1+\xi)-\pi^2(\nu+1)\ln{2}+\frac{7\zeta(3)}{2}(\nu+1)-\frac{4\nu}{\xi}\chi_3(\xi)\right)\\
&=\chi_3\left(\frac{\nu-1}{\nu+1}\right)-\frac{7\zeta(3)}{8}+\frac{\pi^2}{4}\ln\left(1+\frac{1}{\nu}\right)
\end{align}
where $\displaystyle\chi_s(z)=\sum_{n\ \text{odd}}\frac{z^n}{n^s}=\frac{1}{2}\left({\rm Li}_s(z)-{\rm Li}_s(-z)\right)$ is the Legendre-chi function.
Integrating back,
\begin{align}
I_\nu(\nu^{-1},1)
&=\underbrace{\frac{\pi^2}{4}\ln\left(\frac{(1+\nu)^{1+\nu}}{\nu^\nu}\right)-\frac{7\zeta(3)}{8}\nu}_{\text{Let this be C}}+2\int^\frac{1-\nu}{1+\nu}_1\frac{\chi_3(v)}{(1+v)^2}{\rm d}v\\
&=C-\left.\frac{2\chi_3(v)}{1+v}\right|^\frac{1-\nu}{1+\nu}_1+2\int^\frac{1-\nu}{1+\nu}_1\frac{\chi_2(v)}{v(1+v)}{\rm d}v\\
&=C+(1-\nu)\chi_3\left(\frac{1-\nu}{1+\nu}\right)-\frac{7\zeta(3)}{8}-\left.\color{white}{\frac{1}{1}}2\chi_2(v)\ln(1+v)\right|^\frac{1-\nu}{1+\nu}_1+\int^\frac{1-\nu}{1+\nu}_1\frac{\ln(1+v)\ln\left(\frac{1+v}{1-v}\right)}{v}{\rm d}v\\
&=C+(1-\nu)\chi_3\left(\frac{1-\nu}{1+\nu}\right)-\frac{7\zeta(3)}{8}+2\chi_2\left(\frac{1-\nu}{1+\nu}\right)\ln\left(\frac{1+\nu}{2}\right)+\frac{\pi^2}{4}\ln{2}\\
&\ \ \ \ +\frac{1}{2}\int^\frac{1-\nu}{1+\nu}_1\frac{\ln^2(1+v)-\ln^2(1-v)+\ln^2\left(\frac{1-v}{1+v}\right)}{v}{\rm d}v
\end{align}
Repeatedly integrating by parts, it is not hard to derive that
\begin{align}
\frac{1}{2}\int^\frac{1-\nu}{1+\nu}_1\frac{\ln^2(1+v)}{v}{\rm d}v
=&\ -\frac{1}{6}\ln^3\left(\frac{1+\nu}{2}\right)+\frac{7\zeta(3)}{8}-{\rm Li}_3\left(\frac{1+\nu}{2}\right)+{\rm Li}_2\left(\frac{1+\nu}{2}\right)\ln\left(\frac{1+\nu}{2}\right)\\
&\ +\frac{1}{2}\ln\left(\frac{1-\nu}{2}\right)\ln^2\left(\frac{1+\nu}{2}\right)\\
-\frac{1}{2}\int^\frac{1-\nu}{1+\nu}_1\frac{\ln^2(1-v)}{v}{\rm d}v
=&\ {\rm Li}_3\left(\frac{2\nu}{1+\nu}\right)-{\rm Li}_2\left(\frac{2\nu}{1+\nu}\right)\ln\left(\frac{2\nu}{1+\nu}\right)-\frac{1}{2}\ln\left(\frac{1-\nu}{1+\nu}\right)\ln^2\left(\frac{2\nu}{1+\nu}\right)\\
\frac{1}{2}\int^\frac{1-\nu}{1+\nu}_1\frac{\ln^2\left(\frac{1-v}{1+v}\right)}{v}{\rm d}v
=&\ -2\chi_3(\nu)+2\chi_2(\nu)\ln{\nu}+\frac{1}{2}\ln^2\nu\ln\left(\frac{1-\nu}{1+\nu}\right)
\end{align}
Therefore, we get, for $I(\nu^{-1},1)$,
\begin{align}
I(\nu^{-1},1)
=&\color{brown}{\ (1-\nu)\chi_3\left(\frac{1-\nu}{1+\nu}\right)-2\chi_3(\nu)+{\rm Li}_3\left(\frac{2\nu}{1+\nu}\right)-{\rm Li}_3\left(\frac{1+\nu}{2}\right)-\frac{7\zeta(3)}{8}\nu}\\
&\ \color{brown}{+2\chi_2(\nu)\ln{\nu}+2\chi_2\left(\frac{1-\nu}{1+\nu}\right)\ln\left(\frac{1+\nu}{2}\right)-{\rm Li}_2\left(\frac{2\nu}{1+\nu}\right)\ln\left(\frac{2\nu}{1+\nu}\right)}\\
&\ \color{brown}{+{\rm Li}_2\left(\frac{1+\nu}{2}\right)\ln\left(\frac{1+\nu}{2}\right)+\frac{\pi^2}{4}\ln\left(\frac{(1+\nu)^{1+\nu}}{\nu^\nu}\right)+\frac{\pi^2}{4}\ln{2}}\\
&\ \color{brown}{+\frac{1}{2}\ln^2\nu\ln\left(\frac{1-\nu}{1+\nu}\right)+\frac{1}{2}\ln\left(\frac{1-\nu}{2}\right)\ln^2\left(\frac{1+\nu}{2}\right)-\frac{1}{6}\ln^3\left(\frac{1+\nu}{2}\right)}\\
&\ \color{brown}{-\frac{1}{2}\ln\left(\frac{1-\nu}{1+\nu}\right)\ln^2\left(\frac{2\nu}{1+\nu}\right)}
\end{align}
If no mistakes were made, this formula should hold for $0<\nu\le1$. Further simplifications to the formula may be possible through some polylogarithm identities.
Best Answer
Cool problem. Both proposed asymptotic forms are correct. I'm only going to do a first order calc for both, and for $J_n$ I'll say where the proof needs some work. For $I_n$ use the expression given by Arathoon's first note. Observe that $t\cot(t)$ is unimodal decreasing and at $t=0$ has the value 1. Thus by raising it to a high power we expect it to become more sharply peaked. In fact it is gaussian because $$ t\cot(t) = 1-t^2/3-t^4/45... \sim \exp(-t^2/3)(1-7/90t^4 + ...) $$ Therefore $$I_n = \int_0^{\pi/2} \frac{dt}{\sin^2t} \big(1-(t\cot(t))^n \big) \sim \int_0^{\pi/2} \frac{dt}{\sin^2t} \big(1-\exp{(-n\,t^2/3)} \big).$$ Because $n \to \infty$ and because of the $\csc^2t,$ it is seen that the most important region is near $t=0.$ The trick now is to recognize that $\tan^2t = t^2 + O(t^4).$ Replace $t^2$ in the exponential with $\tan^2t$ because it just so happens (Mathematica knows it too) that $$ \int_0^{\pi/2} \frac{dt}{\sin^2t} \big(1-\exp{(-a\,\tan^2t)} \big)=\sqrt{a \pi}.$$ With $a=n/3,$ we indeed have $I_n \sim \sqrt{n \pi /3}.$
For $J_n$ I used a technique call $\textit{depoissonization}.$ Make an exponential power series, $$\sum_{n=0}^\infty \frac{y^n}{n!}J_n = \sum_{n=0}^\infty\frac{y^n}{n!} \int_0^{\infty}(1-t\cot^{-1}(t))^n\, dt = e^y\int_0^{\infty}\exp{(-y\,t\,\cot^{-1}{t} )} \, dt. $$ Since $t\cot^{-1}t$ is monotonically increasing, all we need is the behavior near $0$ to asymptotically estimate the integral. In fact, $t\cot^{-1}t = \pi\,t/2 +O(t^2).$ Using the first term we therefore find $$J(y) := e^{-y}\sum_{n=0}^\infty \frac{y^n}{n!}J_n \sim \int_0^{\infty} \exp{(-\pi\,y\,t/2)} \, dt = \frac{2}{\pi\,y}.$$ By depoissonization, $J_n \sim J(n) =2/(\pi n)$ as long as the next term in the sequence is smaller than this first term. I don't intend to show that. One thing that makes this problem interesting is that in fact it can be shown that $I_n = \sum_{k=1}^n(-1)^{k+1}\binom{n}{k}J_k$ so that $I_n$ and $J_n$ are binomial transforms. Putting in the first order asymptotic for for either $I_n$ or $J_n$ will $\textit{not}$ give you the asymptotics of its transform.