Alternatively Q-CC2 was exhaustively employed for TPA calculations, though it has been reported to exhibit a general overestimation trend (by factor between 1.4 and 2) as compared to Q-CC3, which might however become smaller when considering medium-sized molecular systems. Its computational scaling is O (𝑁 ), yet the resolution of identity (RI) scheme implemented for instance in Turbomole helps improving its computational efficiency, making Q-CC2 applicable to quite large systems. An alternative to quadratic response CC methods is provided by the expectation-value approach, widely used for property computations, with the equation-of-motion CC theories (E-CC) or the algebraic diagrammatic construction in the intermediate state representation (I-ADC) schemes. An efficient E-CCSD method is implemented in Q-Chem employing both the RI and Cholesky decomposition (CD) schemes, which allows significant time saving and reduced memory requirements compared to standard Q-CCSD, with a decrease of the performances of about 5% only on small molecules. This method has been successfully applied as a reference for studying medium-sized systems, e.g., stilbene, phenanthrene, as well as fluorescent protein chromophores. Recently, Nanda and Krylov, developed the so-called Double Electron Attached CCSD (DEA-CCSD) method, with the aim of improving the accuracy of E-CCSD for states having a double excitation character, such as the 𝐴 𝑔 states of ethene, butadiene, and hexatriene. Within the ADC family the I-ADC2 and I-ADC(3/2) methods, respectively scaling O (𝑁 5 ) and O (𝑁 7 ), have been implemented in Q-Chem by Knippengerg and collaborators. Using Q-CC3 references for hydrogen fluoride, water, and ethylene, they demonstrate that I-ADC2 performs slightly better than Q-CC2. I-ADC(3/2) outperforms both second-order methods, but it results in an average percentage error four times larger than its Q-CCSD counterpart. It is also worth mentioning that methods employing the polarization propagator approach, deriving from a Hermitian theory, prevent the occurrence of nonphysical behavior in the calculation of properties. Although Hermitian coupled-cluster (CC) formulations have been proposed, as far as the authors knowledge no implementation of quadratic response theory allowing for computing two-photon transition strengths is currently available. As a further alternative method, let us underline the recent implementation of the computation of nonlinear optical properties with the Bethe-Salpeter equation approach by Rauwolf and co-workers. The quadratic response formalism can also be applied to compute TPA strengths using Time-Dependent Density Functional Theory (TD-DFT), such approach is currently employed in various software packages, e.g., Dalton, GAMESS, Turbomole, and VeloxChem. TD-DFT can be routinely employed for absolute TPA calculations in medium and large dyes. As always, the accuracy of TD-DFT estimates strongly depends on the selected exchange correlation functional (XCF). It has indeed been observed how the choice of a proper XCF is particularly challenging for computing the TPA, as standard XCF considerably underperforms with respect to CC methods. CAM-B3LYP is often considered suitable for 𝛿 TPA , as it reproduces the Q-CC2 trends, although it generally underestimates Q-CC2's 𝛿 TPA by a factor of 3. Recent studies suggest that novel highly parameterized XCF, such as MN15, improve CAM-B3LYP results for push-pull systems. Besides quadratic response has been recently implemented for large scale DFT methods, enabling a qualitative estimation of TPA strengths on extensive systems (up to thousands of atoms).