Matrix exponential, containing a thermal state This question was originally posted on MSE, and I'm cross posting it here. 
Define an infinite matrix $$ M = 
\begin{bmatrix}
0 & -1 & 0 & 0 & \cdots \\
1 & 0 & -2 & 0 & \cdots \\
0 & 2 & 0 & -3 & \cdots \\
0 & 0 & 3 & 0 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{bmatrix}$$
Numerically, I've found that the first column of $\exp(M)$ is given by $\alpha(1,e^{-\lambda},e^{-2\lambda},e^{-3\lambda},\dots)^T$, where $\lambda \approx 0.27$ and $\alpha = \sqrt{1-e^{-2\lambda}} \approx 0.65$. 

Question: How can we prove analytically that the first column of $\exp(M)$ has the stated form?


Context: This question comes from a quantum mechanical model of two harmonic oscillators coupled by a Hamiltonian of the form $\hat{H} \propto \hat{a}_1^\dagger \hat{a}_2^\dagger - \hat{a}_1 \hat{a}_2$, where $\hat{a}_i$ is the lowering operator for oscillator $i$.  The matrix $M$ is essentially this Hamiltonian on the subspace spanned by $\{\left|nn\right>\}$, and so $\exp(M)$ is a time evolution operator.  The exponentially decaying coefficients in the first column indicates that either of the individual harmonic oscillators (after tracing out the other one) is in a thermal state.
I've tried several tricks with the BCH formula, but was discouraged by the fact that the commutators are not very cooperative.  
 A: It is actually easier to compute the first column of $e^{tM}$ for each $t \in \mathbb{R}$ rather than computing only the first column of $e^M$.
Indeed, let $u(t) = e^{tM}e_1$ for each $t \in \mathbb{R}$, where $e_1$ denotes the first canonical unit vector. Then $u(t)$ is the first column of $e^{tM}$.
Claim: For every $k \in \mathbb{N} := \{1,2,\dots\}$ and every $t \in \mathbb{R}$ we have
\begin{align*}
  u_k(t) = \frac{1}{\cosh(t)} \tanh(t)^{k-1}. \tag{*}
\end{align*}
Proof. Let $v_k(t)$ denote the right handside of $(*)$ for each $k$ and each $t$. A short computation shows that $\dot v_k(t) = (k-1)v_{k-1}(t) - k v_{k+1}(t) = (Mv(t))_k$ for each $k \in \mathbb{N}$ and each $t \in \mathbb{R}$ (where we assign an arbitary value to $v_0(t)$ to make sure that the equality also makes sense for $k = 1$). Hence, $v(t)$ solves the differential equation $\dot v(t) = Mv(t)$. Moreover, we have $v(0) = e_1$, so indeed $v(t) = e^{tM}e_1 = u(t)$ for all $t \in \mathbb{R}$.
Remarks:
(1) In the above "proof" I was a bit sloppy about properties of "unique solutions of differential equations". More precisely: (i) I only computed the derivative of $v(t)$ componentwise, which does not show that the mapping $\mathbb{R} \ni t \mapsto v(t) \in \ell^2$ is differentiable with respect to the norm in $\ell^2$. (ii) The "matrix" $M$ is not a continuous linear mapping, but an unbounded linear operator on $\ell^2$, so to make precise assertions about "unique solvability" of initial value problems, one needs to employ the theory of $C_0$-semigroups. 
However, given the context of the question, I don't think that this is the most relevant point here, so I guess there's no need to go into detail here.
(2) Concerning the question how to find the solution presented in $(*)$: The numerical computation of the OP suggests that $u(t)$ is of the form $u_k(t) = g(t) f(t)^{k-1}$ for two functions $f$ and $g$. The differential equation for $u(t)$ yields the differential equation
\begin{align*}
  \dot f(t) = 1 - \frac{k}{k-1} f(t)^2 - \frac{1}{k-1} \frac{\dot g(t)}{g(t)} f(t) \tag{A}
\end{align*}
for $k \ge 2$ and the differential equation
\begin{align*}
  -f(t) = \frac{\dot g(t)}{g(t)} \tag{B}
\end{align*}
for $k = 1$ (provided that $g(t)$ is non-zero). Plugging $(B)$ into $(A)$, we obtain $\dot f(t) = 1 - f(t)^2$, which is solved by $\tanh(t)$. Then we can plug this into $(B)$ again and compute $g(t) = \frac{1}{\cosh(t)}$.
(3) By setting $t = 1$ we obtain $u_k(1) = \frac{1}{\cosh(1)} \tanh(1)^{k-1}$, which is exactly the solution suggested by the OP's numerical experiment (where $\alpha = \frac{1}{\cosh(1)} \approx 0.6481$ and $\lambda = -\ln(\tanh(1)) \approx 0.2723$).
A: Using the numerical findings in Jochen Glueck's answer I'd like to add the following comment.
Let P denote the lower triangular Pascal-matrix; $V(x)$ a vector containing the consecutive powers of $x$ :$V(x)=[1,x,x^2,...]$ and  $\,^dV(x)$ its arrangement as diagonal matrix
Then let Q denote the similarity scaling of P : 
  $$Q = \,^dV(\sinh(1)) \cdot P \cdot \,^dV(\sinh(1))^{-1}
  $$
Then your matrix-exponential E can be described as LU-composition of Q with its inverted transpose
$$ E = \,^dV(1/\cosh(1)) \cdot Q \cdot (Q^{-1})^\tau \cdot \,^dV(1/\cosh(1))  \cdot \cosh(1)
$$
(It is easy to formulate this in Pari/GP and to see how this approximates well the matrix E when computed by the eponential-series with your initial matrix as argument)
