& \vdots & \vdots & \vdots\\
0 & 0 & \ldots &\lambda_n & b_{n1}&\ldots&b_{nm}\\
b_{11}& b_{21} & \ldots &b_{n1} & d_{11} &\ldots & d_{1m}\\
\vdots& \vdots & \ldots &\vdots & \vdots &\ldots & \vdots\\
b_{1m} & b_{2m} & \ldots & b_{nm} & d_{m1} & \ldots & d_{mm}
\end{matrix}
\right]
\left[
\begin{matrix}
x_1\\
x_2\\
\vdots\\
x_n\\
y_1\\
\vdots\\
y_m
\end{matrix}
\right]
=
\epsilon
\left[
\begin{matrix}
x_1\\
x_2\\
\vdots\\
x_n\\
y_1\\
\vdots\\
y_m
\end{matrix}
\right]~,$$ in which the $b$ and $d$ coefficients represent the extensions to the cluster constraint, which depend on the variable $p$, the $\lambda_i$’s are the precalculated eigenvalues of the $\mathcal{T}_1$ on the cluster, and $\epsilon$ is the lowest eigenvalue of the total matrix. This leads to $n+m$ linear equations: $$\begin{aligned}
\lambda_i x_i + \sum_{j=1}^m b_{ij} y_j =& \epsilon x_i\qquad\text{for}\qquad i = 1\ldots n~,\\
\sum_{i = 1}^n b_{ij} x_i + \sum_j d_{ij} y_i =& \epsilon y_j\qquad\text{for}\qquad j = 1\ldots m~.\end{aligned}$$ From the first $n$ equations the $x_i$’s can be eliminated and expressed as a function of the $y$’s: $$x_i = \frac{1}{\epsilon - \lambda_i}\sum_{k=1}^m {b_{ik}}y_k~.$$ When this is substituted into the remaining $m$ equations we get the following self-consistent eigenvalue problem: $$\sum_{k=1}^m\Sigma(\epsilon)_{jk}~y_k = \epsilon y_j~,\quad\text{for}\quad j = 1\ldots m~,\quad\text{with}\quad \Sigma(\epsilon)_{jk} 