\delta_{\nu\nu'} - \frac{1}{2} \sum_{p_1h_1} D_{p_1h_1}X_{p_1h_1}^{\nu',{J}^{\pi}} 
                                                       X_{p_1h_1}^{\nu ,{J}^{\pi} \ast}  
] 
D_{ph}^{1/2}D_{p'h}^{1/2}Y_{ph}^{\nu ,{J}^{\pi}} Y_{p'h}^{\nu',{J}^{\pi}\ast } 
, 
\label{E:OBDMp} 
\\ 
\langle 0 | a_{h}^{\dagger} a_{h'} | 0 \rangle 
&=& 
\delta_{hh'} -  
\hat{j}_h^{-2}    \sum_{{J}^{\pi} \nu \nu' p  }
\hat{J}^2          [ 
\delta_{\nu\nu'} 
- \frac{1}{2} \sum_{p_1h_1} D_{p_1h_1}X_{p_1h_1}^{\nu',{J}^{\pi}} 
                                                       X_{p_1h_1}^{\nu ,{J}^{\pi} \ast}  
] 
D_{ph}^{1/2}D_{ph'}^{1/2}Y_{ph}^{\nu ,{J}^{\pi}} Y_{ph'}^{\nu',{J}^{\pi}\ast } 
, 
\label{E:OBDMh} \end{aligned}$$

where $\hat{j}\equiv \sqrt{(2j+1)}$. These equations, which are exact up to terms $O(|Y|^4)$, are not linear in the OBDM elements. Therefore, for given $X$ and $Y$ amplitudes (at a given step of the ERPA iteration), they have to be solved iteratively as well: one begins with an initial guess for the occupation probabilities on the r.h.s., calculates the OBDM on the l.h.s., diagonalizes it, thus obtaining a new set of occupation probabilities, and so on, until convergence is reached. Let us call this the OBDM iteration.

This double iterative procedure of solving the ERPA equations can be outlined as follows: In order to begin the ERPA iteration, we need an initial choice for the ground state. We start with the HF ground state. We build the ERPA matrix for each ${J}^{\pi}$ that we wish to take into account by using Eqs. (\[E:Amat\]) and (\[E:Bmat\]) and assuming the HF occupation probabilities 