We provide a (partial) solution in these sections
- closed-form expressions
- numerical results
- asymptotic behaviour
- discussion
Closed-form expressions
Since the original integral is evaluated by Mathematica only for the first two values we shall look here for equivalent expressions for the integral and see if Mathematica performs better with these.
Equivalent expressions
First, for easier reading, we give the formulas in Latex, then we write down the Mathematica expressions.
Equivalent expressions for $f(n)$ are (the proof of the equivalences is left to the reader. I confine myself to brief hints, the complete proofs can be provided on request).
$$f1(n)={\frac{1}{\sqrt{\pi }}\int_{-\infty }^{\infty } u \exp \left(-u^2\right) \left(\frac{1}{2}\text{erfc}(u)\right)^{n-1} \, du}$$
Noticing that $\text{erfc} = 1 - \text{erf}$ and expanding the power into a binomial sum gives
$$f2(n) = - \frac{2^{1-n}}{\sqrt{\pi }} \sum _{m=1}^{\left\lfloor n/2\right\rfloor } g(2 m-1) \binom{n-1}{2m-1}$$
where
$$g(k)={\int_{-\infty }^{\infty } u \exp \left(-u^2\right) \text{erf}(u)^{k} \, du}$$
Notice that $g(k) = 0$ for even $k$. We have already taken this into account in $f2$ where only odd terms contribute.
The following forms are equivalent to $g$ for odd $k$ but they will be considered also for even $k$ in the following
$$g1(k)={2\int_{0 }^{\infty } u \exp \left(-u^2\right) \text{erf}(u)^{k} \, du}$$
Letting $u^2 \to t$ gives
$$g2(k)={\int_{0 }^{\infty } \exp \left(-t \right) \text{erf}(\sqrt{t})^{k} \, du}$$
The next form is a multiple integral over the $k$-dimensional hypercube obtained by inserting the integral representation of $\text{erf}$ into $g2$ and performing the $t$-Integration:
$$g3(k)=\left(\frac{2}{\sqrt{\pi }}\right)^k \Gamma \left(\frac{k}{2}+1\right) \int _0^1 ... \int _0^1\frac{1}{\left(\sum _{i=1}^k y(i)^2+1\right){}^{\frac{k}{2}+1}}dy(1) ... dy(k)$$
The corresponding Mathematica expressions (with the same names) are:
f1[n_] :=
1/Sqrt[\[Pi]]
Integrate[
u Exp[-u^2] (1/2 Erfc[u])^(n - 1), {u, -\[Infinity], \[Infinity]}]
f2[n_] :=
2^(1 - n)/Sqrt[\[Pi]]
Sum[Binomial[n - 1, k] (-1)^k g[k], {k, 1, n - 1,2}]
g[k_] := Integrate[u Exp[-u^2] Erf[u]^k, {u, -\[Infinity], \[Infinity]}]
g1[k_] := 2 Integrate[u Exp[-u^2] Erf[u]^k, {u, 0, \[Infinity]}]
g2[k_] := Integrate[Exp[-t] Erf[Sqrt[t]]^k, {t, 0, \[Infinity]}]
g3[k_] :=
(2/Sqrt[\[Pi]])^k Gamma[1 + k/2] Integrate[1/(1 + Sum[y[i]^2, {i, 1, k}])^(1 + k/2), Sequence @@ Table[{y[i], 0, 1}, {i, 1, k}]]
Results
We have shown that the calculation of $f$ can be traced back to evaluating the function $g$.
It turns out that $g3$ gives closed-form results for the highest $k$:
Table[{k, g3[k]}, {k, 1, 3}]
$$\left(
\begin{array}{cc}
1 & \frac{1}{\sqrt{2}} \\
2 & \frac{2 \sqrt{2} \tan ^{-1}\left(\frac{1}{\sqrt{2}}\right)}{\pi } \\
3 & \frac{3 \sqrt{2} \cot ^{-1}\left(2 \sqrt{2}\right)}{\pi } \\
\end{array}
\right)$$
Unfortunately, this "highest $k$ is still very modest. Already $g3(4)$ remains unevaluated.
Looking for a rule in the first 3 terms I came up with
gguess[k_] := k Sqrt[2] ArcTan[1/Sqrt[k^2-1]])/\[Pi]}
but it doesn't work for $k=4$ as can be seen by numerical comparison.
Finally, using $f2$ with $g3$ (called $f23$) we obtain the closed-form expressions for $f(n)$
{#, f23[#]} & /@ {1, 2, 3, 4, 5} // Simplify
$$\left(
\begin{array}{cc}
1 & 0 \\
2 & -\frac{1}{2 \sqrt{2 \pi }} \\
3 & -\frac{1}{2 \sqrt{2 \pi }} \\
4 & -\frac{3 \cos ^{-1}\left(-1/3\right)}{4 \sqrt{2} \pi ^{3/2}} \\
5 & -\frac{3 \sqrt{2}}{4 \pi ^{3/2}}\cos ^{-1}\left(\sqrt{\frac{2}{3}}-\frac{1}{6}\right)\\
\end{array}
\right)$$
numerically
% // N
(* Out[797]= {{1., 0.}, {2., -0.199471}, {3., -0.199471}, {4., -0.18197}, {5., -0.164468}} *)
Numerical results
Define the numerical analogue to our functions as
fn[n_] := NIntegrate[s^(n - 1) InverseErf[1 - 2 s], {s, 0, 1}]
f1n[n_] := 1/Sqrt[\[Pi]]
NIntegrate[u Exp[-u^2] (1/2 Erfc[u])^(n - 1), {u, -\[Infinity], \[Infinity]}]
As expected $f1n$ is faster than $fn$;
AbsoluteTiming[fn[10^4]]
(* Out[430]= {1.88997, -0.00027235} *)
AbsoluteTiming[f1n[10^4]]
(* Out[431]= {0.022427, -0.00027235} *)
The first few discrete values are (excluding the trivial case $n = 1$), using
tfn[nn_] := Table[{n, fn[n] // Chop}, {n, 2, nn}]
tfn[10]
(* Out[676]=
{
{2, -0.199471}, {3, -0.199471}, {4, -0.18197}, {5, -0.164468}, {6,-0.149342}, {7, -0.136591}, {8, -0.12583}, {9, -0.116674}, {10, -0.108806}
}
Notice the agreement with the numerical values of the closed-form expressions for $n = 2..4$ which is also a verification of these.
Plotting it
ListPlot[tfn[20],
PlotLabel -> "Numerical values of f(n) for n = 2 .. 20",
AxesLabel -> {"n", "f(n)"}]
Asymptotic behaviour
The asymptotics for large $n$ is most conveniently calculated in the original form of $f$.
The integrand of $f$ is given by
fi[n_, s_] := s^(n - 1) InverseErf[1 - 2 s]
The "exotic" function InverseErf[] does not seem to have an official name (please correct my if I am wrong), hence for brevity we shall call it fre(x) and define
fre[x_] := InverseErf[x]
The function $\text{fre}$ is, of course, just $\text{erf}$ rotated by $\pi /4$ and looks like this
Plot[fre[x], {x, -1, 1}, PlotLabel -> "The function fre(x) = InverseErf[x]",
AxesLabel -> {"x", "fre(x)"}, PlotRange -> All]
The complete integrand for some small values of $n$ is shown here
Plot[{fi[2, s], fi[3, s], fi[5, s]}, {s, 0, 1},
PlotLabel -> "The integrand for some n\nn = 2 (blue), 3 (yellow), 5 (green)",
AxesLabel -> {"s", "fi(n,s)"}, PlotRange -> {-2, .2}]
We see that, as expected, for increasing $n$ the appreciable contribution to the integral concentrates around $s = 1$
Close to $s = 1$ the integrand is
Simplify[Series[s^(n - 1) fre[1 - 2 x], {x, 1, 1}] // Normal, 0 < x < 1]
(* Out[733]=
-s^(-1 + n) Sqrt[-Log[1 - x] -
1/2 Log[2 \[Pi] (Log[1/(2 \[Pi])] - 2 Log[1 - x])]]
*)
First let us look at the plot of $f$ for very large $n$ and try to find analytic expression suggested by the series expansion
Plot[{1, - f1n[n] n/Sqrt[(-Log[1/n] - Log[-Log[1/n]])], -f1n[n] n/Sqrt[
ProductLog[n]], - f1n[n] n/Sqrt[Log[n] ]}, {n, 2, 10^6},
PlotRange -> {0.85, 1.05}, AxesLabel -> {"n", "f/fappr"},
PlotLabel ->
"Behavior of f for large n\nrelative to various approximations\nfappr ~ \
asympt(ProductLog)/n (yellow)\nfappr ~ \!\(\*SqrtBox[\(ProductLog[n]\)]\)/n \
(green)\nfappr ~ \!\(\*SqrtBox[\(Log[n]\)]\)/n (red)"]
We see that, up to a factor close to unity the asymptotics is best described by
$$f(n)\sim \frac{1}{n}\sqrt{W(n)}$$
where $W(n)$ is the LambertW function (Lambert 1758, http://mathworld.wolfram.com/LambertW-Function.html ) which is the inverse of $n = W \exp (W)$.
An analytic derivation of the asymptotics will be given later.
Discussion
Relation to Laplace transform
The function $g2$ is in fact a Laplace transform
$$glp(n) = \mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^n\right](s)$$
at $s\to 1$.
The function $g$ in the form $g3$ is easily generalized to the Laplace form by replacing the $1$ before the sum by $s$, and it can be calculated for the powers 1, 2, and 3. Adding the recent discovery of Mariusz Iwaniuk for $g(4)$ (which, however, contains an unevaluated integral) we have
$$\left(
\begin{array}{cc}
1 & \frac{1}{s \sqrt{s+1}} \\
2 & \frac{4 \tan ^{-1}\left(\frac{1}{\sqrt{s+1}}\right)}{\pi s \sqrt{s+1}} = \frac{4}{\pi s \sqrt{s+1}}\sin ^{-1}\left(\frac{1}{\sqrt{s+2}}\right) \\
3 & \frac{6 \tan ^{-1}\left(\frac{1}{\sqrt{s+1} \sqrt{s+3}}\right)}{\pi s \sqrt{s+1}} = \frac{6}{\pi s \sqrt{s+1}}\sin ^{-1}\left(\frac{1}{s+2}\right)\\
4 & \frac{4 \sqrt{\frac{1}{s+1}}}{s}-\frac{96 \int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{a^2+2}\right)}{\left(a^2+1\right) \sqrt{a^2+2}} \, da}{\pi ^2 s \sqrt{s+1}}\\
\end{array}
\right)$$
In the published tables of Laplace transformations (see e.g. http://authors.library.caltech.edu/43489/1/Volume%201.pdf 4.12 (4)) I have found only the case $n = 1$, hence, as a byproduct of this study, we can add two closed-form expressions and one with an integral left to the tables.
EDIT 06.02.2017
Mariusz Iwaniuk pointed out in a comment that
$$g2(4,s) = \frac{4 \sqrt{\frac{1}{s+1}}}{s}-\frac{96 \int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{a^2+2}\right)}{\left(a^2+1\right) \sqrt{a^2+2}} \, da}{\pi ^2 s \sqrt{s+1}}$$
and that the remaining integral is related to Ahmed's integral (https://arxiv.org/pdf/1411.5169.pdf).
The numerical check $g2(4,s=1) = g2n(4)$ is ok.
The integral we need is not Ahmed's integral but a generalization of it. The upper integration limit is not unity but $\sqrt{1+s}$.
EDIT 07.02.2017
In an attempt to "crack" the remaining integral I have defined some integrals appearing naturally in this context, and have derived relations between these (they include the relation of my comment).
The integrals are
$$M(s)=\int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{y^2+2}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy\tag{7.1}$$
$$B(s)= \int_0^1 \frac{\cot ^{-1}\left(\frac{\sqrt{y^2+2}}{\sqrt{s+1}}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy\tag{7.2}$$
$$C(s) = \frac{1}{2} \int_0^s \frac{\cot ^{-1}\left(\sqrt{z+3}\right)}{\sqrt{z+1} (z+2) \sqrt{z+3}} \, dz\tag{7.3}$$
Note added on 15.02.17:
While trying to understand the recent EDIT of Mariusz Iwaniuk I found that the integral C can be transformed into a simpler form by the substitution $z\to \frac{1}{\sin (w)}-2$
$$C_1(s) = \frac{1}{2} \int_{\sin ^{-1}\left(\frac{1}{s+2}\right)}^{\frac{\pi }{6}} \cot ^{-1}\left(\sqrt{\csc (w)+1}\right) \, dw\tag{7.3.1}$$
The relations are
$$M(s)=\frac{1}{4} \pi \tan ^{-1}\left(\sqrt{s+1}\right)-B(s) \tag {7.4}$$
$$B(s)=\frac{1}{4} \pi \tan ^{-1}\left(\sqrt{s+1}\right)-\frac{\pi ^2}{32} -C(s)\tag{7.5}$$
and, combining these two gives the interesting relation
$$M(s)=\frac{\pi ^2}{32} + C(s)\tag{7.6}$$
It turns out that, for the case of the OP, i.e. $s=1$, the main part of $M(1)$ is given by the explicit fraction, and only about 13% are due to integral $C(1)$.
In fact, numerically for $s=1$ the equation reads (in the above order)
0.3540215 = 0.3084251+0.0455963, and the relation is 0.0455963/0.3540215 = 0.128795.
We can find a good approximate expression for $C$. In fact, since $\cot^{-1}$ is a decreasing function we can replace it in the integrand by its maximum value $ \cot^{-1}\left(\sqrt{3}\right)$ to obtain an upper bound (i.e. the function in question is smalller). Now the integral can be performed, and the result for $M$ becomes
$$M(s) \lesssim \frac{\pi ^2}{288}+\frac{1}{6} \pi \tan ^{-1}\left(\sqrt{\frac{s+1}{s+3}}\right)$$
Notice: I have the derivation and some corollaries ready in draft state, but I have to beg for your patience as it takes some time for me to type it in.
EDIT 08.02.17: Generating function and difference-differential equation
Remember that we wish to calculate the Laplace-transformation of the n-th power of the Erf-function of a square root
$$glp(n,s)=\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^n\right](s)$$
Here's a new idea: instead of looking into single powers we try to "untie the Gordian knot" and cover them all at once.
This is generally done using generating functions. It is convenient here to use the exponential generating function
$$h(z, s)=\int_0^{\infty } \exp (-s t) \exp \left(-z\; \text{erf}\left(\sqrt{t}\right)\right) \, dt\tag{8.1}$$
After some manipulations including partial integration combined with differentiation I came up with this equation for the generating function
$$ddeq=s (s+1) \frac{\partial h(z,s)}{\partial s}+\left(\frac{3 s}{2}+1\right) h(z,s) -\frac{z^2 h(z,s+2)}{\pi }=\frac{1}{2}\tag{8.2}$$
This is a linear inhomogeneous partial difference-differential equation of first order.
This type of equation is new to me. Hence I don't know yet how to approach a solution.
Some first observations
1) the "generating" parameter $z$ appears only in one position as a factor
2) if $z = 0$, the eqation reduces to a simple ODE which can readily be solved with a simple result
DSolve[ddeq /. z -> 0, h[0, s], s]
(* {{h[0, s] -> 1/s + C[1]/(s Sqrt[1 + s])}} *)
The first part is simply $glp(n=0,s)$, as it should be. It is interesting that the next term $glp(n=1,s)$ also appears automatically.
EDIT 18.10.17: Derivation of the difference-differential equation
It's my pleasure to comply with Mariusz Iwaniuk's recent request and present the details of the derivation of equ. (8.2).
The text needs some formatting, sorry.
The derivation consists of a chain of partial integrations and derivatives.
§1)
We consider the expression
D[Exp[-s t] Exp[-z Erf[Sqrt[t]]], t] // Expand
-E^(-s t - z Erf[Sqrt[t]]) s - (E^(-t - s t - z Erf[Sqrt[t]]) z)/(
Sqrt[\[Pi]] Sqrt[t])
Integrating over t gives
Exp[-s t] Exp[-z Erf[Sqrt[t]]] /. t -> 0
1
Limit[Exp[-s t] Exp[-z Erf[Sqrt[t]]], t -> \[Infinity], Assumptions -> s > 0]
0
Hence
1 == s Integrate[E^(-s t - z Erf[Sqrt[t]]) , {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]]
Integrate[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]
(* 1 *) 1 ==
s h[z, s] +
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]
Differentiating the equation with respect to s gives
0 = D[s h[z, s], s] +
z/Sqrt[\[Pi]]
Integrate[
D[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], s], {t, 0, \[Infinity]}]
but
D[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], s]
-E^(-t - s t - z Erf[Sqrt[t]]) Sqrt[t]
inserting this gives
(* 2 *) 0 =
D[s h[z, s], s] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}]
§2 2nd partial integration
Consider
(notice that we have s+1 now)
D[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t] // Expand
E^(-(1 + s) t - z Erf[Sqrt[t]])/(2 Sqrt[t]) -
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t] -
E^(-(1 + s) t - z Erf[Sqrt[t]]) s Sqrt[t] - (
E^(-t - (1 + s) t - z Erf[Sqrt[t]]) z)/Sqrt[\[Pi]]
Integrating this expression gives zero since
Limit[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t -> 0]
0
Limit[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t -> \[Infinity],
Assumptions -> s > 0]
0
0 == Integrate[E^(-(1 + s) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t, 0, \[Infinity]}] - (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 2) t - z Erf[Sqrt[t]]) , {t, 0, \[Infinity]}]
or, noticing that the last integral is h(z,s+2), i.e. h(z,s) with s shifted by two units,
(* 3 *) Integrate[E^(-(1 + s) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t,
0, \[Infinity]}] == (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]] h[z, s + 2]
Now the left hand side of (* 3 ) can be replaced using ( 1 *)
(* 1 rep. *) 1 ==
s h[z, s] +
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]
or
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t, 0, \[Infinity]}] == Sqrt[\[Pi]]/(2 z) (1 - s h[z, s])
giving
Sqrt[\[Pi]]/(
2 z) (1 - s h[z, s]) == (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]] h[z, s + 2]
Now on the r.h.s. of (* 3 *) we use (*2 *)
(* 2 rappel *) 0 =
D[s h[z, s], s] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}]
or
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] ==
Sqrt[\[Pi]]/z D[s h[z, s], s]
to get from (* 3 *)
(* 3a *) Sqrt[\[Pi]]/(
2 z) (1 - s h[z, s]) == (1 + s) Sqrt[\[Pi]]/z D[s h[z, s], s] +
z/Sqrt[\[Pi]] h[z, s + 2]
simplifying and sorting
1/2 (1 - s h[z, s]) == (1 + s) D[s h[z, s], s] + z^2/\[Pi] h[z, s + 2]
1/2 (1 - s h[z, s]) == (1 + s) (h[z, s] + s D[h[z, s], s]) +
z^2/\[Pi] h[z, s + 2]
1/2 (1 - s h[z, s]) == (1 + s) h[z, s] + (1 + s) s D[h[z, s], s] +
z^2/\[Pi] h[z, s + 2]
1/2 - s /2 h[z, s] - (1 + s) h[z, s] -
z^2/\[Pi] h[z, s + 2] == (1 + s) s D[h[z, s], s]
gives finally the equation
(* 4 *) ddeq = (1 + s) s D[h[z, s], s] + (1 + (3 s )/2) h[z, s] -
z^2/\[Pi] h[z, s + 2] == 1/2;
End of derivation of (8.2).
n = 2wil be: $-\frac{1}{2 \sqrt{2 \pi }}$ $\endgroup$