In view of the high activity in experimental non-linear X-ray spectroscopies, the paucity of robust ab initio methods capable of describing the underlying non-linear properties limits the potential impact of new experimental techniques, especially given the above-mentioned contributions from competing processes in the experimental signal. Being the simplest higherorder process from the theoretical point of view, X2PA can be used as a platform for future developments of theoretical methods for modeling higher-order multiphoton processes, such as X-ray sum-frequency generation and X-ray second-harmonic generation. 2PA cross sections can be computed using different approaches. Formally, 2PA moments can be obtained from the first residues of quadratic response functions (first electricdipole hyperpolarizability tensor) to build the 2PA strength tensor and cross sections. Alternatively, elements of the 2PA strength tensor can be evaluated as first residues of the cubic response functions (second electricdipole hyperpolarizability tensor). In practice, these relationships are exploited in the resonance convergent complex polarization propagator (CPP) approach and in damped response theory. In both these approaches, a phenomenological imaginary inverse lifetime (damping) parameter is introduced to compute the 2PA strengths from the imaginary part of the cubic response function. The CPP approach has been exploited within the framework of time-dependent DFT (TDDFT) for computing both valence 2PA 62 and X2PA spectra; Ref. 63 by Fahleson et al. being the only published theoretical work on X2PA until now. As described below, we here explore a computational strategy based on the many-body EOM-CC framework.