Overall, this study highlights the necessity for improved cost-efficient approximations for computing TPA crosssections in large systems. Indeed, the current low-scaling (less then O (𝑁 6 )) state-of-the-art methods appears to have variable accuracies and it is not clear what would be the accuracy in complex calculations, indicating a need for further advancements in wavefunction methods and/or XCF. In this regard, the recent development of quadratic response for BSE/GW presents a promising avenue for future exploration.