Closed, sum-free form for the $n$-th derivative of $\operatorname{arcsinh}(\frac1x)$ in $x=1$ During research involving the Born–Jordan quantization I came across the expression
$$
\frac{d^k}{dx^k}\operatorname{arcsinh}\Big(\frac1x\Big)\Big|_{x=1}\tag1
$$
for $k\in\mathbb N_0$. It is not too hard to write this expression as a sum 
$$
(1)=\frac{\sqrt{2}}{2^k}\sum_{j=0}^{k-1} a_j^k\tag{2a}
$$
for any $k\in\mathbb N$ where $(a_j^k)_{j\in\mathbb Z,k\in\mathbb N}$ is a recursive sequence of integers given by
$$
a_j^k:=\begin{cases} a_0^1=-1&\\a_j^k=0&\text{if }j<0\text{ or }j\geq k\\ a_j^{k+1}=a_j^k(2j-k)+a_{j-1}^k(2j-3k-1)&\text{else}\end{cases},\tag{2b}
$$
(basically a modified version of Pascal's triangle). Unfortunately, I so far was not able to find a closed (sum-free) form of $(1)$ / $(2a)$ for arbitrary $k\in\mathbb N_0$.
This recursive sequence is nice for explicit calculations (especially speeds up things for larger $n$) - but I'm interested how $(1)$, or rather of the arising matrix elements
$$
M_{nn}:=2\sum_{k=0}^n\begin{pmatrix}n\\k\end{pmatrix}\frac{2^k}{k!}\Big(\frac{d^k}{dx^k}\operatorname{arcsinh}\Big(\frac1x\Big)\Big|_{x=1}\Big)\tag3
$$
behave for large $n$.

Explicit calculations (up until $n=200$) suggest that $M_{nn}\overset{n\to\infty}\longrightarrow0$ with $M_{nn}=\mathcal O(\frac1n)$ for $n\to\infty$. What could be an approach to potentially prove this? What is more realistic: trying to find some bound for $(3)$ or trying to find a closed, sum-free form for $(1)$, respectively $(2a)$ / $(2b)$?

Being fairly new here I hope this "question" (or rather problem) is suitable for mathoverflow. If it is not, feel free to tell me so I can outsource it to math.stackexchange. Thanks in advance for any answer or comment!
 A: Here is a proof that $M_{nn}=O(n^{-1})$, as the OP conjectured. The function
$$f(z):=\operatorname{arcsinh}(1/z),\qquad |z|<1,$$
is holomorphic in the unit disk, hence by Cauchy's formula
$$ \frac{f^{(k)}(1)}{k!}=\frac{1}{2\pi i}\int_{|z-1|=r}\frac{f(z)}{(z-1)^{k+1}}\,dz,\qquad 0<r<1.$$
Here and later circles are positively oriented. It follows, by applying the binomial formula, that
$$\frac{1}{2}M_{nn}=\frac{1}{2\pi i}\int_{|z-1|=r}\frac{f(z)}{z-1}\left(\frac{z+1}{z-1}\right)^n\,dz,\qquad 0<r<1.$$
We make the change of variables $z:=(1+w)/(1-w)$, keep track of the image of the original circle (which is another circle), and deform it to a circle centered at the origin. This way we obtain
$$\frac{1}{2}M_{nn}=\frac{1}{2\pi i}\int_{|w|=\rho}\frac{g(w)}{1-w}\cdot\frac{dw}{w^{n+1}},\qquad 0<\rho<1,$$
where $g(w)$ is the following holomorphic function in the unit disk:
$$g(w):=\operatorname{arcsinh}\left(\frac{1-w}{1+w}\right),\qquad |w|<1.$$
This way we see, again by Cauchy's formula, that $M_{nn}$ is twice the sum of the first $n$ Taylor coefficients of $g(w)$ at the origin. (Added: This can also be derived without complex analysis, as in Max Alekseyev's response.) Using the identities
$$g'(w)=\frac{-\sqrt{2}}{(1+w)(1+w^2)^{1/2}},\qquad|w|<1,$$
$$\frac{1}{(1+w^2)^{1/2}}=\sum_{k=0}^\infty\frac{(-1)^k}{4^k}\binom{2k}{k}w^{2k},\qquad|w|<1,$$
it is straightforward to derive that
$$\frac{1}{2}M_{nn}=\operatorname{arcsinh}(1)-\sqrt{2}\sum_{2k\leq n}\frac{(-1)^k}{4^k}\binom{2k}{k}\sum_{2k\leq m\leq n-1}\frac{(-1)^m}{m+1}.$$
The inner sum is positive and it is of size $O(k^{-1})$, hence in the outer sum the $k$-th term has sign $(-1)^k$ and it is of size $O(k^{-3/2})$. It follows that $M_{nn}$ converges. By applying a more careful asymptotic analysis in $k$ and keeping track of the error term $O(n^{-1})$ coming from the inner sum, the convergence is seen to occur with speed $O(n^{-1})$. Finally, by Abel's theorem, 
$$\lim_{n\to\infty} M_{nn}=2g(1)=0,\qquad\text{hence in fact}\qquad M_{nn}=O(n^{-1}).$$
A: The numbers $c_k = \frac{1}{k!} \left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1}$ are the coefficients in the expansion:
$$\operatorname{arcsinh}\frac1x = c_0 + c_1(x-1) + c_2(x-1)^2 + \dots.$$
It follows that $2^kc_k$ is the coefficient of $t^k$ in $\operatorname{arcsinh}\frac1{1+2t}$.
Now,
$$M_{nn} = 2\sum_{k=0}^n \binom{n}{k} 2^kc_k = 2\cdot [t^n]\ (1+t)^n \operatorname{arcsinh}\frac1{1+2t}.$$
Using Lagrange inversion, we get the generating function for $M_{nn}$:
$$\sum_{n\geq 0} M_{nn} t^n = \frac{2}{1-t} \operatorname{arcsinh}\frac{1-t}{1+t}.$$
From here the asymptotic for $M_{nn}$ can be obtained using the standard tools (e.g., see the answer from GH from MO).

To get an explicit formula for $\left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1}$ (and thus for $c_k$), we notice that 
$$\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x = -  \left(\frac{d}{dx}\right)^{k-1} (x^2+x^4)^{-\frac12}.$$
To expand the last expression, one can use Faà di Bruno's formula:
$$\left(\frac{d}{dx}\right)^{k-1} (x^2+x^4)^{-\frac12}$$
$$ = \sum_{i=1}^{k-1} \binom{-\frac{1}{2}}{i} i! (x^2+x^4)^{-\frac12-i} B_{k-1,i}(2x+4x^3,2+12x^2,24x,24,0,0,\dots,0),$$
where $B_{k-1,j}$ are Bell polynomials.
Evaluating at $x=1$, for $k>0$, we get 
$$\left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1} =
-\sum_{i=1}^{k-1} \binom{-\frac{1}{2}}{i} i! 2^{-\frac12-i} B_{k-1,i}(6,14,24,24,0,0,\dots,0)$$
$$ = -\frac{(k-1)!}{\sqrt{2}}\sum_{j_1+2j_2+3j_3+4j_4=k-1} \frac{(-1)^{j_1+j_2+j_3+j_4}(2(j_1+j_2+j_3+j_4))!}{(j_1+j_2+j_3+j_4)!j_1!j_2!j_3!j_4!} 2^{-3(j_1+j_2+j_3+j_4)} 6^{j_1}7^{j_2}4^{j_3}$$
$$=-\frac{(k-1)!}{\sqrt{2}}\sum_{j_1+2j_2+3j_3+4j_4=k-1} \frac{(-1)^{j_1+j_2+j_3+j_4} (2(j_1+j_2+j_3+j_4))!}{(j_1+j_2+j_3+j_4)!j_1!j_2!j_3!j_4!} 2^{-2j_1-3j_2-j_3-3j_4} 3^{j_1}7^{j_2}.$$
A: Let 
$$
f(x):=\operatorname{arcsinh}\Big(\frac1x\Big). 
$$
Then 
\begin{equation}
 -f''(x)=\Big(\frac1x+\Re\frac1{x+i}\Big)f'(x). 
\end{equation}
So, by Leibniz's formula, for $k=0,1,\dots$ 
\begin{equation}
 -f^{(k+2)}(x)=\sum_{j=0}^n\binom kj (-1)^j j!\Big(\frac1{x^{j+1}}+\Re\frac1{(x+i)^{j+1}}\Big)f^{(k-j+1)}(x); 
\end{equation}
so, for $a_j:=f^{(j)}(1)$ we have the one-index recurrence 
\begin{equation}
  a_{k+2}=\sum_{j=0}^n\binom kj (-1)^{j+1} j!\big(1+2^{-(j+1)/2}\Re e^{i(j+1)\pi/4}\big)a_{k-j+1}, 
\end{equation}
from which I think it should be not too hard to get asymptotics of $a_n$ and then maybe of $M_{nn}$. 
