Bilinear transformation of continuous time state space system
I've had the same question last week, but I've managed to find how to derive it (getting rid of those $z$ terms is indeed tricky). I will give here detailed demonstration of how to arrive to the result given in 1 (with, in your notation, $\alpha = 2 \lambda$).

So we define our new discrete-time function transfer as

$$
\begin{array}{rcl}
H_d(z) &amp;=&amp; H(\frac{\alpha\left(z-1\right)}{z+1}) \\
&amp;=&amp; \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&amp;=&amp; \frac{z+1}{\alpha} \mathbf{C} \left[ \left(z-1\right)\mathbf{I}- \frac{z+1}{\alpha}\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&amp;=&amp; \frac{z+1}{\alpha} \mathbf{C} \left[ z \left(\mathbf{I} - \frac{1}{\alpha} \mathbf{A} \right) - \left(\mathbf{I} + \frac{1}{\alpha} \mathbf{A}\right) \right]^{-1}\mathbf{B} + \mathbf{D} \\
\end{array}
$$

For sake of notation, let $\mathbf{P} = \mathbf{I} - \frac{1}{\alpha} \mathbf{A}$ and $\mathbf{Q} = \mathbf{I} + \frac{1}{\alpha} \mathbf{A}$ (important remark useful later on: $\mathbf{P} + \mathbf{Q} = 2 \mathbf{I}$). Then

$$
\begin{array}{rcl}
H_d(z) &amp;=&amp; \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\
&amp;=&amp; \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\
&amp;=&amp; \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \left( \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} \right) + \mathbf{D} \\
&amp;=&amp; \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D}
\end{array}
$$
where we defined
$ \mathbf{A_d} = \mathbf{P}^{-1} \mathbf{Q} $
and
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}}  \mathbf{P}^{-1} \mathbf{B} $.
Furthermore, we have

$$
\begin{array}{rcl}
\frac{z+1}{\sqrt{2 \alpha}} \mathbf{C}
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \mathbf{P}^{-1} \left( \mathbf{P} + \mathbf{Q} \right) \right] \\
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right]
\end{array}
$$

which gives us

$$
\begin{array}{rcl}
H_d(z)
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&amp;=&amp; \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} + \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\
&amp;=&amp; \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} 
\end{array}
$$
with
$ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} $
and
$ \mathbf{D_d} = D + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} $.
Replacing $\mathbf{P}$ and $\mathbf{P}$ gives the results in 1.

PS: the same idea can be used to prove the discrete-time state-space representation found using  Generalized Bilinear Transform ($s \leftarrow \alpha \frac{z - 1}{\beta z + \left(1-\beta\right)}$) or First-order Holder methods.