As a further analysis we investigated the results for valence states with large 𝛿 TPA (larger than 20 a.u. at the CC3/aug-cc-pVDZ level), since they are the more interesting for actual applications. We underline that this subset includes only 9 transitions, so these results should be interpreted with caution. The results of MAPE and 𝑅 2 are reported in Fig. , while the other indicators can be found in Tables S11-S15 in the SI. Q-CCSD is again the most adequate method with a MAPE of 24.4% and a 𝑅 2 of 0.93. Interestingly, for this subset, Q-CC2 and E-CCSD show similar performances: MAPE of 34.3% (44.8%), and 𝑅 2 of 0.88 (0.86) for CC2 (E-CCSD). We underline that the performances of Q-CCSD, E-CCSD and Q-CC2 are mainly driven by the 1𝐵 3𝑔 state of tetrazine. In fact its Q-CC3 𝛿 TPA of 56.0 a.u., is largely overestimated by both Q-CCSD (86.6 a.u.) and E-CCSD (118.3 a.u.), while CC2 provides a rather accurate estimate of 47.7 a.u. Indeed if we exclude this molecule the MAPE becomes: 21.1% for Q-CCSD, 37.4% for E-CCSD and 37.1% for Q-CC2. I-ADC2 still performs quite poorly with a MAPE of 96.7% and MAE of 43.1, but has a relatively acceptable 𝑅 2 of 0.76. TD-DFT methods perform similarly as in the whole set in terms of MAPE but they show extremely poor 𝑅 2 values. B3LYP is the "best" XCF, with a MAPE of 44.5%, and a low 𝑅 2 of 0.34. The three range-separated functionals show similar results: MAPEs of 47.4%, 47.0%, and 48.3%; and 𝑅 2 values of 0.17, 0.21, and 0.16 for CAM-B3LYP, LC-BLYP33, and LC-BLYP47, respectively. BH&HLYP performs the worst among the TD-DFT methods, with a MAPE of 63.2%, and an extremely low 𝑅 2 of 0.04. This of course raises questions for actual TD-DFT applications for TPA simulations. The performances of the wavefunction methods Q-CCSD, E-CCSD, Q-CC2, I-ADC2, and DEA-CCSD (which raw values are reported in Table in the SI) have been evaluated for the six electronic states with significant contributions from double excitations. TD-DFT methods were excluded due to their inability to account for double excitations, and any ac-curate reproduction of 𝛿 TPA by TD-DFT would likely be due to significant error cancellation. Table reports the absolute percentage error (APE) for these methods compared to reference Q-CC3 calculations. Values with APEs exceeding 200% are not detailled. The DEA-CCSD method provides notably accurate results for the 1𝐴 1 state of cyclopentadiene and the 1𝐴 𝑔 state of butadiene, with APEs of 15% and 23%, respectively. These results significantly improve upon those provided by both Q-CCSD (96% and 118%) and E-CCSD (129% and 173%). In contrast, Q-CC2 and I-ADC2 show poor performances, with errors exceeding 200%. For the 1𝐴 𝑔 state of hexatriene, DEA-CCSD is the only method that yields a reasonably accurate result (APE of 62%), whereas all other methods provide very poor values. On the other hand, for the two states of acrolein DEA-CCSD appears quite unreliable, indeed we have an APE larger then 200%. For the 2𝐴" state, Q-CCSD is the best method with an APE of 65%, followed by E-CCSD (79%), I-ADC2 (96%), and Q-CC2 (115%). For the 5𝐴' state, the second order methods surprisngly outperform the CCSD-based methods, with Q-CC2 and I-ADC2 providing better results than both Q-CCSD and E-CCSD. All methods perform reasonably well for the 1𝐸 2𝑔 state of benzene, with E-CCSD showing the lowest APE (4%), and DEA-CCSD has the highest APE (63%).