## 2 Theoretical Formulation
### Atomic Structure Calculation  
The non-relativistic two-electron ionic Hamiltonian, \(H_{A}\), in atomic units, is given by:  
\[\hat{H}_{A}=\hat{h}(\mathbf{r}_{1})+\hat{h}(\mathbf{r}_{2})+\frac{1}{|\mathbf{ r}_{1}-\mathbf{r}_{2}|}, \tag{1}\]  
where \(\hat{h}(\mathbf{r}_{i})=-\
abla_{i}^{2}/2-Z/r_{i}\) and \(\mathbf{r}_{i}\), \(i=1,2\), denotes each electron's coordinate. \(-\
abla_{i}^{2}/2\) represents the electronic kinetic energy operator, \(-Z/r_{i}\) the nucleus-\(i^{th}\) electron electrostatic (Coulombic) potential (where it is assumed that the nucleus defines the origin of the working coordinate system) and \(1/|\mathbf{r}_{1}-\mathbf{r}_{2}|\) represents the corresponding inter-electronic electrostatic potential. \(Z\) is the atomic number, which is 3 for Li\({}^{+}\), 10 for Ne\({}^{8+}\) and 18 for Ar\({}^{16+}\). First, we numerically solve the Schrodinger equation (SE), \((\hat{h}(\mathbf{r})-\epsilon)\phi_{elm_{i}}(\mathbf{r})=0\), for the mono-electronic systems Li\({}^{2+}\), Ne\({}^{9+}\) and Ar\({}^{17+}\). To this end, we adopt a separation of variables approach where the one-electron orbitals are expressed as \(\phi_{elm_{i}}(\mathbf{r})=[P_{e}(r)/\rho]Y_{lm_{i}}(\theta,\phi)\) where \(Y_{lm_{i}}(\theta,\phi)\), are the spherical harmonic functions. Projection of \(\phi_{elm_{i}}(\mathbf{r})\) onto the one-electron SE, followed by angular integration leads to the one-dimensional radial differential equation for the unknown radial orbitals, \(P_{e}(r)\):  
\[\left[-\frac{1}{2}\frac{d^{2}}{dr^{2}}+\frac{l(l+1)}{2r^{2}}-\frac{Z}{r}- \epsilon\right]\;P_{e}(r)=0. \tag{2}\]  
The present method assumes that the radial configuration space of the electron is limited to a finite radius, \(R\), with the boundary conditions chosen as: \(P_{e}(0)=P_{e}(R)=0\). The immediate consequence of this assumption is to allow for a finite matrix representation of the physical Hamiltonian with a discretized eigenergy spectrum including both the bound (\(\epsilon<0\)) and the continuum eigenstates (\(\epsilon>0\)). Moreover, we also see that both the eigenenergies and eigenstates are dependent on the particular angular momentum number, \(l\). In the following, we adopt the following discretized notation for the eigenvalues and the eigenfunction: \(\varepsilon\rightarrow\varepsilon_{nl}\) and \(P_{\varepsilon}(r)\to P_{nl}\). The positive-energy eigenstates, \(\varepsilon_{nl}\geq 0\), exhibit oscillatory asymptotic behaviour, similar to that expected from continuum-like eigenstates. On the other hand, the negative-energy eigenstates, \(\varepsilon_{nl}<0\), have an exponentially decaying asymptotic behaviour associated with the bound spectrum of the finite Hamiltonian. The numerical solution then proceeds by expanding the radial orbitals, \(P_{nl}(r)=\sum_{i}c_{i}^{(nl)}B_{i}^{(k_{b})}(r)\), on a nonorthogonal set of B-spline polynomials of order \(k_{b}\) and total number \(n_{b}\) defined in the finite interval \([0,\,R]\)[24]. The particular choice of a B-spline basis, as opposed to a Gaussian or Slater-type basis, is dictated by their superior ability to represent the continuum solutions with great accuracy. This expansion leads to a diagonalization matrix problem, where the solution provides the unknown coefficients, \(c_{i}^{(nl)}\)[25]. Having calculated the radial orbitals, \(P_{nl}(r)\), for each partial wave \(l=0,1,2,...\), we proceed with the solution of the SE for the two-electron ions Li\({}^{+}\), Ne\({}^{8+}\) and Ar\({}^{16+}\):  
\[\hat{H}_{A}\Psi_{EL}(\mathbf{r}_{1},\mathbf{r}_{2})=E_{L}\Psi_{EL}(\mathbf{r}_ {1},\mathbf{r}_{2}). \tag{3}\]  
The calculational method proceeds along similar lines to the one-electron case.