## 2 Theoretical Formulation
### Two- and Three-Photon Ionization Cross Section Formulation  
According to LOPT, the \(N\)-photon partial-wave ionization cross section, following the absorption of \(N\) photons of energy \(\omega\), from a system in its ground state \(|g\rangle\) (of energy \(E_{g}\)) to a final continuum state of energy \(E\) and angular momentum \(L\), is given by [10; 27]:  
\[\sigma_{L}^{(N)}(\omega)=2\pi(2\pi\alpha)^{N}\omega^{N}|M_{E_{L}}^{(N)}|^{2} \;\delta(E-E_{g}-N\omega), \tag{6}\]  
where \(\alpha\) is the fine structure constant. The total \(N\)-photon cross section is given by \(\sigma_{N}(\omega)=\sum_{L}\sigma_{L}^{(N)}(\omega)\). The \(N\)-photon transition amplitude, \(M_{E_{L}E_{L}}^{(N)}\), in the case of two-photon absorption, \(N=2\), reduces to:  
\[M_{E_{L}}^{(2)}=\sum_{E_{p}^{\prime}}\frac{D_{g;E_{p}^{\prime}}D_{E_{p}^{ \prime};E_{L}}}{E_{g}+\omega-E_{p}^{\prime}},\qquad L=S,D, \tag{7}\]  
while for three-photon absorption, \(N=3\), it is given by:  
\[M_{E_{p}}^{(3)}=\sum_{E_{p}^{\prime}}\frac{D_{g;E_{p}^{\prime}}} {E_{g}+\omega-E_{p}^{\prime}}\left[\sum_{E_{p}^{\prime}}\frac{D_{E_{p}^{\prime };E_{p}^{\prime}}D_{E_{p}^{\prime};E_{p}}}{E_{g}+2\omega-E_{g}^{\prime}}+\sum_ {E_{p}^{\prime}}\frac{D_{E_{p}^{\prime};E_{D}^{\prime}}D_{E_{p}^{\prime};E_{p} }}{E_{g}+2\omega-E_{D}^{\prime}}\right], \tag{8}\] \[M_{E_{p}}^{(3)}=\sum_{E_{p}^{\prime}}\sum_{E_{p}^{\prime}}\sum_{E _{p}^{\prime}}\frac{D_{E_{g}E_{p}^{\prime}}D_{E_{p}^{\prime};E_{p}^{\prime}}D _{E_{p}^{\prime};E_{p}}}{(E_{g}+\omega-E_{p}^{\prime})(E_{g}+2\omega-E_{p}^{ \prime})}, \tag{9}\]  
where \(D_{E_{L};E_{L}^{\prime}}\equiv\langle\Psi_{EL}|\hat{D}|\Psi_{EL^{\prime}}\rangle\) is the two-electron dipole matrix element and \(\Psi_{EL}=\Psi_{EL}(\mathbf{r}_{1},\mathbf{r}_{2})\) is the CI two-electron eigenstate given by Equation (4). The integrals over the angular symmetries include both the bound and continuum states. The length form of the interaction operator is \(\hat{D}=-\hat{e}_{p}\cdot(\mathbf{r}_{1}+\mathbf{r}_{2})\) and the velocity form is \(\hat{D}=\hat{e}_{p}\cdot(\
abla_{1}+\
abla_{2})/\omega\), with \(\hat{e}_{p}\) representing the polarization unit vector of the radiation. The detailed expressions of the dipole matrix elements in terms of the calculated CI two-electron wavefunctions can be found in Ref. [23].  
At this point, it is appropriate to comment on the adoption of the dipole approximation used in evaluating transition amplitudes. In the present context of the X-ray regime, the use of the dipole approximation is justified by considering the scaling of the dipole transitions of the one-electron (hydrogenic) systems (see Equation (2)). The mean distance of the electron from the ionic core scales with the nuclear charge, \(Z\), as \(\langle r\rangle\sim 1/Z\), while the ionization potential scale as \(IP(Z)\sim Z^{2}\). In the present work, we assume transitions with X-ray photons below the first ionization threshold, \(\omega<|E(1s)-E(1s^{2})|\). Since generally \(|E(1s)-E(1s^{2})|<IP(Z)\), we thus have \(\omega<Z^{2}\). The dipole approximation is based on the validity of \(|\mathbf{k}_{\gamma}\cdot\mathbf{r}|<<1\), where \(\mathbf{k}_{\gamma}\) is the wavevector of the X-ray photon and \(\mathbf{r}\) the electron's position. From the above considerations, we then have \(|\mathbf{k}_{\gamma}\cdot\mathbf{r}|\simeq k_{\gamma}\langle r\rangle=\omega \langle r\rangle/c\sim Z^{2}/(Zc)=Z/c<<1\), for all \(Z=3,10,18\) that we consider here (\(c\simeq 137.036\)).  
The angular momentum of the final states, following photoabsorption (\(L=0,2\) for two-photon absorption and \(L=1,3\) for three-photon absorption), is dictated by the selection rules implicit in the dipole matrix element. For the present case of linearly polarized light, we have \(\Delta M_{L}=0\) and \(\Delta L=\pm 1\)[26]. Since we start from the ion's ground state (\(M_{L}=0\)), we only need to keep the transitions between the \(M_{L}=0\) states.