I'm interested in the following problem which arises from some "random matrix theory" calculations. Let $\phi,s_1,s_2, p > 0$ with $p \in [0,1]$, and set $p_1=p$, $p_2=1-p$, and $q_k := p_k s_k$ for $k=1,2$. Let $\Sigma$ be a $d \times d$ positive-definite matrix which has a limiting spectral density $\nu$ which is compactly supported on $\mathbb R_+$. Consider the equations \begin{align} \eta &= r(\theta):= \overline{\operatorname{tr}} \Sigma(\Sigma + \theta I_d)^{-1},\\ \frac{1}{\theta} &= \frac{q_1}{\lambda + s_1\phi\theta\eta} + \frac{q_2}{\lambda + s_1 \phi\theta\eta}. \end{align} Let $\theta(\lambda)$ be the unique nonnegative solution.

Question. In terms of $\phi,p_1,p_2,s_1,s_2,\nu$, what is an analytic / algebraic description of $\theta(0^+):= \lim_{\lambda \to 0^+}\theta(\lambda)$ ?

My hope is that there is a deep mathematical result which allows one to provide such a description (by which I mean an explicit formula, some fixed-point equation, etc.)

Example. Suppose $\Sigma = \sigma^2 I_d$. Then, $r(\theta) = \sigma^2 / (\sigma^2 + \theta)$. If we further assume $s_1=s_2 = 1$, then the equations become $$ \frac{1}{\theta}= \frac{1}{\lambda + \phi \theta \sigma^2 / (\sigma^2 + \theta)}. $$ This is a quadratic equation with unique nonnegative solution given by $$ \theta(\lambda) = \frac{\lambda + (\phi-1)\sigma^2 + \sqrt{(\lambda + (\phi-1)\sigma^2)^2+4\lambda \sigma^2}}{2}. \tag{1} $$ It follows that $\theta(0^+) = \sigma^2\max(\phi-1,0)$.

Context: The expression $\theta(\lambda)$ given in (1) above is related to the Stieltjes transform of the celebrated Marchenko-Pastur law.