Aug-cc-pVDZ was employed for the benchmarks, coupled with the corresponding auxiliary basis sets when necessary. States with 𝛿 TPA Q-CC3 < 3 a.u. are considered irrelevant for the purposes of the paper and they have not been included in the statistics, their values are reported in the SI for completeness. States for which TPA transitions are forbidden have not been reported. For such reason acetylene and streptocyanine-1 have been excluded from the statistical analysis, indeed the 1𝐴 𝑢 and 2𝐴 𝑢 states of the former are TPA forbidden, while the 1𝐵 2 state of the latter (the only state studied) have 𝛿 TPA value below assumed threshold. All CC calculations use the frozen-core approximation, with the length gauge for the interaction operator. Standard convergence criteria (10 -7 a.u.) was applied for the ground state energies and CC equations in Dalton which was used for computing Q-CC3, Q-CCSD, Q-CC2, and TD-DFT values. E-CCSD (RI-CD-EOM-CCSD) calculations were performed using Q-Chem (v6.0) and applying the following parameters: 10 -9 a.u. for ground-state energies, 10 -8 a.u. for CC equations, and 10 -7 a.u. for EOM Davidson convergence. The CD of two-electron integrals uses the standard threshold, 10 -3 a.u., as it was observed to be sufficient for the TPA strengths. DEA-CCSD calculations have been performed on systems with a significant double excitation character on Q-Chem using the same thresholds as its EOM-CCSD version (indeed the RI-CD aproximations are not yet implemented for this method), and employing the charge +2 for reproducing the double excitations. I-ADC2 calculations are achieved with Q-Chem using default convergence criteria. The coupled cluster implementations employed in Dalton and Q-chem is non-Hermitian, therefore, left and right second-order transition moments in Eq. 2 are generally different and in the worst cases they might even present different signs. Since the product of these two moments appears in Eq. 2, unphysical negative terms may arise, impacting the total TPA strength. To address this, cases where one or more products of left and right moments in Eq. 2 are negative and represent at least 20% of the total value of 𝛿 TPA are excluded from the analysis. The excluded cases are documented in Table in the SI. We stress that I-ADC2 and TD-DFT formulations, originating from Hermitian theories, are free of this difficulty. Finally, states expected to exhibit resonanceenhanced TPA effects are also excluded and detailed in Table in the SI. These (high-energy) states present an excitation energy twice that of a lower excited-state, potentially leading to divergent TPA strengths, which can be cured only through damped response theory, which was not considered here.