Limit of difference of pseudo-arcsinh and arcsinh $$\DeclareMathOperator{\erf}{erf}$$
$$\DeclareMathOperator{\arcsinh}{arcsinh}$$
Let $f$ a function defined by : $$f(x) = \int_0^x \frac{\erf(u)}{u}du$$Where $\erf$ is the error-function. This function look like $\arcsinh(x)$ :

In black there is $f(x)$ and in red this is $\arcsinh(x)$. Now, we can see (or at least conjecture) that $$\kappa := \lim_\limits{x \rightarrow +\infty} \left[ \int_0^x \frac{\erf(u)}{u}du - \arcsinh(x)\right]$$exists and is finite. That would imply that $f(x) \sim \kappa + \arcsinh(x)$ and more accuratly $f(x) \approx \kappa + \arcsinh(x)$. By using $\arcsinh(x) \sim \ln(2x)$ we would have $$f(x) \sim \kappa + \ln(2x)$$Now, heres the graph of $f(x)-\arcsinh(x)$ :

In black this is $f(x)-\arcsinh(x)$. We can easily see that it converges toward $\kappa$ (in red). Here, we have $\kappa \approx 0.288607829951$.
My question is, what is the exact value of $\kappa$ ?
 A: $$f(x) = \int_0^x \frac{\text{erf}(u)}{u}du=\frac{2 x }{\sqrt
   {\pi }}\,
   _2F_2\left(\frac{1}{2},\frac{1}{2};\frac{3}{2},\frac{3}{2};-x^2\right)$$Using series for large values of $x$
$$f(x)-\sinh ^{-1}(x)=\frac{\gamma }{2}-\frac{1}{4 x^2}+O\left(\frac{1}{x^3}\right)$$ which is not much.
In fact, $\gamma$ appears because the first terms of the expansion of $f(x)$ are
$$\log \left({x}\right)-\frac{1}{2}\psi \left(\frac{1}{2}\right)=\log (2 x)+\frac{\gamma }{2}$$
If you want much closer, instead of $\sinh ^{-1}(x)$ use
$\sinh ^{-1}\left(e^{\gamma /2} x\right)$ which would give
$$f(x)-\sinh ^{-1}\left(e^{\gamma /2} x\right)\sim -\frac{e^{-\gamma }}{4 x^2}$$
This would be a very idea for an approximation of the solution of $f(x)=k$.
The initial estimate would then be
$$x_0=e^{-\gamma /2} \sinh (k)$$
Trying for $k=3.456$, Newton iterates will be
$$\left(
\begin{array}{cc}
n & x_n \\
 0 & 11.86091721 \\
 1 & 11.87273377 \\
 2 & 11.87273966
\end{array}
\right)$$
Edit
We can improve the estimate of the solution of equation $f(x)=k$ at the price of a single evaluation of the hypergeometric function. Using the original Householder method, we should have
$$x_1=x_0+\frac{3 f(x_0) \left(f(x_0) f''(x_0)-2 f'(x_0)^2\right)}{f(x_0)^2 f'''(x_0)+6 f'(x_0)^3-6 f(x_0) f'(x_0) f''(x_0)}$$ with
$$f'(x_0)=\frac{\text{erf}(x_0)}{x_0} \qquad f''(x_0)=\frac{2 e^{-x_0^2}}{\sqrt{\pi } x_0}-\frac{f'(x_0)}{x_0}$$
$$f'''(x_0)=-\frac{4 e^{-x_0^2} \left(x_0^2+1\right)}{\sqrt{\pi } x_0^2}+2\frac{f'(x_0)}{x_0^2}$$
Applied to the working case, this would give
$$x_0=\color{red}{11.8}609172053431$$
$$x_1=\color{red}{ 11.87273965533}87$$ while the "exact" solution is
$$x =\color{red} { 11.8727396553392}$$ and we could continue for better and better using the first iteration of an higher order method.
A: $$\kappa=\underset{x\to \infty }{\text{lim}}\left(\int_0^x \frac{\text{erf}(u)}{u} \, du-\text{arcsinh } x\right)=\frac{\gamma }{2}\approx 0.2886078324507665$$
