I was able to compute a nontrivial lower bound for the Gaussian case (thanks to Ying Xiao who privately proposed the approach). For convenience I will consider the distribution of $\sigma_1\sim{\cal N}(0,\frac1d I)$. Now, I will ignore the nonnegativity constraints on $\lambda$ and consider the following Lagrangian $$ {\cal L}(\lambda,\mu) = \|\Sigma \lambda\|_2^2-2\mu (\mathbf{1}^{\top}\lambda-1).$$ Writing the first-order optimality conditions, we get \begin{eqnarray*} \lambda^{\ast} &=& \mu^{\ast}(\Sigma^{\top}\Sigma)^{-1}\mathbf{1}\\ \mu^{\ast} &=& \dfrac{1}{\mathbf{1}^{\top}(\Sigma^{\top}\Sigma)^{-1}\mathbf{1}}. \end{eqnarray*}
This way $\|\Sigma\lambda^{\ast}\|_2^2=1/[\mathbf{1}^{\top}(\Sigma^{\top}\Sigma)^{-1}\mathbf{1}]$. Therefore, the quantity of interest can be lower bounded by $$\mathbb{E}[\|\Sigma\lambda^{\ast}\|_2]=\mathbb{E}\Big[\frac{1}{\sqrt{\mathbf{1}^{\top}(\Sigma^{\top}\Sigma)^{-1}\mathbf{1}}}\Big].$$ The term inside the expectation squared is known to have a chi-squared distribution with $d-M+1$ degrees of freedom (I found out about this here), thus the expectation above is of order $1/\sqrt{d-M+1}$$\sqrt{\frac{d-M+1}{d}}$.
I believe something similar should be true for the (scaled) Boolean distribution. On the other hand, I don't know how tight is this lower bound, but it turns out it suffices for my problem.