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) &=& H(\frac{\alpha\left(z-1\right)}{z+1}) \\ &=& \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \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} \\ &=& \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) &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\ &=& \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} \\ &=& \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
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} &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\ &=& \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] \\ &=& \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) &=& \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} \\ &=& \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} \\ &=& \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} \end{array} $$ with $ \mathbf{C_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} $ and $ \mathbf{D_d} = \mathbf{D} + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} $. Replacing $\mathbf{P}$ and $\alpha = 2 \lambda$ 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.
