Matrix exponential, containing a thermal state

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$).