with 𝜇, 𝜈 ∈ 𝑥, 𝑦, 𝑧 and 𝑀 𝜇𝜈 𝑖←0 and 𝑀 𝜇𝜈 0←𝑖 representing the left and right second-order transition moments, which are equal in Hermitian theories, see Eq. S1 in the Supporting Information (SI) for their expressions. For most organic molecules, the coupled cluster (CC) methods can be viewed as the most reliable approaches for computing 𝛿 TPA . In addition, the hierarchy of CC methods provides a systematic and well defined path for improving the accuracy by incorporating higher excitation orders in the wavefunction expansion. Among these methods, the approximate coupled cluster singles, doubles, and triples model (CC3), in its quadratic response implementation (Q-CC3), stands out as the most reliable for TPA calculations when a monodeterminantal description of the electronic state is appropriate. However, its computational cost, scaling as 𝑂(𝑁 7 ) (where 𝑁 is the number of orbitals), limits its application to small systems, preventing the study of large molecular assemblies. This is a significant limit since the molecules of interest for actual TPA applications are actually quite large 𝜋-conjugated molecules. Notably, just a few Q-CC3 reference TPA strengths are available in the literature. Paterson and collaborators studied a series of tiny molecules namely, water, ethylene, formaldehyde, diacetylene, whereas Nanda and collaborators computed the TPA of twisted ethene, as well as the lowest excited 𝐴 𝑔 states of butadiene, and hexatriene. Moreover, Siriatayant and Andruniów determined Q-CC3 TPA strengths for the first excited state of a model chromophore for retinal. All of these reference data have been computed using the Dalton code. When Q-CC3 becomes beyond computational reach coupled cluster singles and doubles method (Q-CCSD) is considered an excellent alternative, as Paterson and collaborators showed that Q-CCSD, with a scaling of O (𝑁 6 ) gives a maximum of 20% error on 𝛿 TPA computed at Q-CC3 level.