# Two-Photon Absorption Strengths of Small Molecules : Reference CC3 Values and Benchmarks

Carmelo Naim, Robert Zaleśny, Denis Jacquemin

► To cite this version:

Carmelo Naim, Robert Zaleśny, Denis Jacquemin. Two-Photon Absorption Strengths of Small Molecules: Reference CC3 Values and Benchmarks. Journal of Chemical Theory and Computation, 2024, 20 (20), pp.9093-9106. 10.1021/acs.jctc.4c01090. http://dx.doi.org/10.1021/acs.jctc.6001090

HAL Id: hal-05221765

https://hal.science/hal-05221765v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two-Photon Absorption Strengths of Small Molecules: Reference CC3 Values and Benchmarks

Carmelo Naim,* Robert Zaleśny, and Denis Jacquemin*

<figure><img src="image_2.png" /></figure>

Cite This: J. Chem. Theory Comput. 2024, 20, 9093–9106

<figure><img src="image_3.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_4.png" /></figure>

Metrics & More

<figure><img src="image_5.png" /></figure>

Article Recommendations

<figure><img src="image_6.png" /></figure>

Supporting Information

**ABSTRACT:** We present a large dataset of highly accurate two-photon transition strengths ($\delta^{TPA}$) determined for standard small molecules. Our reference values have been calculated using the quadratic response implementation of the third-order coupled cluster method including iterative triples (Q-CC3). The aug-cc-pVTZ atomic basis set is used for molecules with up to five non-hydrogen atoms, while larger molecules are assessed with aug-cc-pVDZ; the differences due to the basis sets are discussed. This dataset, encompassing 82 singlet transitions of various characters (Rydberg, valence, and double excitations), enables a comprehensive benchmark of smaller basis sets and alternative wavefunction methods when Q-CC3 calculations become beyond reach as well as time-dependent density functional theory (TD-DFT) approaches. The evaluated wavefunction methods include quadratic response and equation-of-motion CCSD approximations, Q-CC2, and second-order algebraic diagrammatic construction in its intermediate state representation (I-ADC2). In the TD-DFT framework, a set of five commonly used exchange-correlation functionals are evaluated. This extensive analysis provides a quantitative assessment of these methods, revealing how different system sizes, response intensities, and types of transitions affect their performances.

<figure><img src="image_7.png" /></figure>

## 1. INTRODUCTION

Two-photon absorption (TPA) is a nonlinear optical process in which an atom or molecule is excited by the simultaneous absorption of two photons. This phenomenon was initially proposed by Göppert-Mayer in 1931¹ and experimentally confirmed only 30 years later by Kaiser and Garrett.² The TPA process can be quantified by the probability of absorption of two photons, which is usually reported as a cross section ($\sigma^{TPA}$). In contemporary applications, materials or dyes exhibiting substantial $\sigma^{TPA}$ are extensively used in various technological fields, including optical data storage, 3D microfabrication, optical power limiting, bioimaging, and photodynamic therapy.³⁻¹⁴ For example, bioimaging applications require the product of the two-photon absorption cross section and fluorescence quantum yield to be larger than 50 GM (1 GM = cm⁴·s·photon⁻¹). This threshold safeguards bright two-photon images at laser powers that will not cause photodamages.¹¹ In the case of organic dyes, $\sigma^{TPA}$ values can reach several thousand GM.¹⁵ During the past five decades, computational quantum chemistry has evolved into an indispensable tool for interpreting experimental TPA measurements.¹⁶,¹⁷ Indeed, theoretical calculations provide valuable insights into structure–property relationships,¹⁶,¹⁸⁻²² and furthermore, they can be routinely employed to design new molecules with significant $\sigma^{TPA}$, obviating the need for generally expensive and complex experiments.¹⁶,²³⁻²⁶ $\sigma^{TPA}$ for the absorption of two photons having the same frequency can be expressed as

$$ \sigma^{\text{TPA}} = \frac{\pi^3 \alpha a_0^5 (2\omega)^2}{c} g(2\omega) \delta^{\text{TPA}} \quad (1) $$

where $\omega$ is the frequency of the incident photons, $\alpha$ is the fine structure constant, $c$ is the speed of light, $a_0$ is the Bohr radius, $g$ is the broadening function (typically Lorentzian or Gaussian, depending on the broadening mechanism²⁷), and $\delta^{TPA}$ is the orientationally averaged two-photon transition strength, which can be determined using various electronic-structure theories (see below). Despite significant developments, an accurate theoretical estimation of $\sigma^{TPA}$ in molecules remains challenging due to the necessity of using advanced quantum chemistry methods, including electron correlation effects and relying on a large basis set.²⁸ Many efforts have been made to find efficient methods capable of computing TPA cross sections in a cost-effective manner. Validation of the results of new computational methodologies can be approached in two ways: by comparison with experimental data or by comparison with benchmark values. The former approach presents several challenges, as reproducing experimental $\sigma^{TPA}$ requires consideration of multiple factors such as solvent effects and

Received: August 20, 2024

Revised: September 20, 2024

Accepted: September 23, 2024

Published: October 7, 2024

<figure><img src="image_8.png" /></figure>
vibrational couplings.<sup>9,21,25,29–37</sup> Moreover, σ<sup>TPA</sup> depends on ω<sup>2</sup> and δ<sup>TPA</sup>, each introducing distinct sources of error in the computations, potentially leading to misleading interpretations.<sup>38–40</sup> The choice of an appropriate broadening function and the proper selection of constants in eq 1 are also subjects of debate.<sup>41,42</sup> For these reasons, direct comparison of δ<sup>TPA</sup> with more accurate (albeit expensive) wavefunction methods is generally preferred to establish the accuracy of computational methods.<sup>28,43–45</sup> In the case of linearly polarized light, for an excited state i, the general formulation of δ<sup>TPA</sup> standing for both Hermitian and non-Hermitian theories is

$$ \delta^{\text{TPA}} = \frac{1}{15} \sum_{\mu,\nu} [M_{i \leftarrow 0}^{\mu\mu} M_{0 \leftarrow i}^{\nu\nu} + M_{i \leftarrow 0}^{\mu\nu} M_{0 \leftarrow i}^{\mu\nu} + M_{i \leftarrow 0}^{\nu\mu} M_{0 \leftarrow i}^{\nu\mu}] \quad (2) $$

where μ, ν ∈ x, y, z and M<sub>i←0</sub><sup>μν</sup> and M<sub>0←i</sub><sup>μν</sup> represent 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).<sup>28</sup> For most organic molecules, coupled cluster (CC) methods can be viewed as the most reliable approaches for computing δ<sup>TPA</sup>. 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 the available CC methods, the approximate coupled cluster singles, doubles, and triples model (CC3)<sup>46</sup> in its quadratic-response implementation (Q-CC3)<sup>28</sup> 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 O(N<sup>7</sup>) (where N is the number of orbitals), limits its application to small systems, preventing the study of large molecular assemblies. This is a significant limitation since the molecules of interest for actual TPA applications are actually quite large π-conjugated molecules. Notably, only 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, and diacetylene,<sup>28</sup> whereas Nanda and collaborators computed the TPA of twisted ethene as well as the lowest excited A<sub>g</sub> states of butadiene and hexatriene.<sup>47</sup> Moreover, Siriatayant and Andruniów determined Q-CC3 TPA strengths for the first excited state of a model chromophore for retinal.<sup>45</sup> All of these reference data were computed using the Dalton code.<sup>48</sup>

When Q-CC3 becomes beyond computational reach, the coupled cluster singles and doubles method (Q-CCSD)<sup>49</sup> is considered as an excellent alternative, as Paterson and collaborators showed that Q-CCSD, with a scaling of O(N<sup>6</sup>), gives a maximum of 20% error on δ<sup>TPA</sup> computed at the Q-CC3 level.

Alternatively, Q-CC2<sup>50</sup> was exhaustively employed for TPA calculations, though it has been reported to exhibit a general overestimation trend (by factor between 1.4 and 2) compared to Q-CC3,<sup>28,41,44</sup> which might however become smaller when considering medium-sized molecular systems.<sup>45</sup> The Q-CC2 computational scaling is O(N<sup>5</sup>), yet the resolution of identity (RI) scheme<sup>51,52</sup> implemented for instance in Turbomole<sup>53,54</sup> helps improving its computational efficiency, making Q-CC2 applicable to quite large systems.<sup>55–57</sup>

An alternative to quadratic-response CC methods is provided by the expectation-value approach,<sup>58</sup> widely used for property computations, with the equation-of-motion CC

(E-CC) theories<sup>58</sup> or the algebraic diagrammatic construction in the intermediate state representation (I-ADC)<sup>43,59</sup> schemes. An efficient E-CCSD method<sup>44</sup> is implemented in Q-Chem<sup>60</sup> employing both the RI and Cholesky decomposition (CD)<sup>61–63</sup> schemes, which allows significant time savings and reduced memory requirements compared to standard Q-CCSD with a decrease in the performance of only about 5% on small molecules.<sup>44</sup> This method has been successfully applied as a reference for studying medium-sized systems, e.g., stilbene and phenanthrene,<sup>38,64</sup> as well as fluorescent protein chromophores.<sup>41,44</sup> 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 A<sub>g</sub> states of ethene, butadiene, and hexatriene.<sup>47</sup> Within the ADC family, the I-ADC2 and I-ADC(3/2) methods, which scale as O(N<sup>5</sup>) and O(N<sup>6</sup>), respectively, have been implemented in Q-Chem by Knippenberg and collaborators.<sup>43</sup> Using Q-CC3 references for hydrogen fluoride, water, and ethylene, they demonstrated that I-ADC2 performs slightly better than Q-CC2.<sup>43</sup> I-ADC(3/2) outperforms both second-order methods, but it results in an average percentage error 4 times larger than its Q-CCSD counterpart.<sup>43</sup> It is also worth mentioning that methods employing the polarization propagator approach, derived from a Hermitian theory, prevent the occurrence of nonphysical behavior in the calculation of properties. Although Hermitian CC formulations have been proposed,<sup>65</sup> to the best of our 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.<sup>66</sup>

The quadratic-response formalism can also be applied to compute TPA strengths using time-dependent density functional theory (TD-DFT),<sup>67,68</sup> and this approach is employed in various software packages, e.g., Dalton,<sup>48</sup> GAMESS,<sup>69</sup> Turbomole,<sup>53,54</sup> and VeloxChem.<sup>70</sup> TD-DFT can be routinely employed for TPA calculations in medium and large dyes.<sup>22,26,56,71–74</sup> 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.<sup>73</sup> CAM-B3LYP<sup>75</sup> is often considered suitable for δ<sup>TPA</sup>, as it reproduces the Q-CC2 trends,<sup>56,73</sup> although it generally underestimates Q-CC2's δ<sup>TPA</sup> by a factor of 3.<sup>40,41,71,72</sup> Recent studies suggest that novel highly parametrized XCFs, such as MN15,<sup>76</sup> improve CAM-B3LYP results for push–pull systems.<sup>40,74,77,78</sup> Besides, quadratic response has been recently implemented for large-scale DFT methods,<sup>39,79,80</sup> enabling a qualitative estimation of TPA strengths on extended systems (up to thousands of atoms).

In this work, with the aim of including TPA data to the QUEST database<sup>81–85</sup> of highly accurate reference values for excited-state energies and properties, we obtained a collection of δ<sup>TPA</sup> computed at the Q-CC3 level for a series of small benchmark molecules. With these data at hand, we benchmark several wavefunction methods which could be used when Q-CC3 becomes unreachable, namely, Q-CCSD, E-CCSD, Q-CC2, and I-ADC2, as well as TD-DFT with XCFs that are frequently used in the study of medium and large systems.
Scheme 1. Schematic Representation of the Molecules Considered in This Study

1,2,3 (non-H) atom molecules

<figure><img src="image_9.png" /></figure>

4 (non-H) atom molecules

<figure><img src="image_10.png" /></figure>

5 (non-H) atom molecules

<figure><img src="image_11.png" /></figure>

6 (non-H) atom molecules

<figure><img src="image_12.png" /></figure>

2. COMPUTATIONAL DETAILS

This study includes 25 organic molecules containing from one to six non-hydrogen atoms, which are shown in Scheme 1. The corresponding geometries, optimized at the CC3/aug-cc-pVTZ level, were extracted from prior works,⁸¹,⁸³ and are provided in the SI. We consider low-lying singlet excited states having valence or Rydberg character as well as several states with non-negligible contributions from double excitations. All molecules studied belong to a symmetry group higher than $C_1$. Therefore, each excited state is identified by its irreducible representation, and their ranking starts from the lowest excited state, excluding the ground state in the numbering.

The importance of double excitations in each state was assessed using the %$T_1$ value, which provides the percentage of contribution of all single excitations in an excited state. States with %$T_1$ smaller than 80% are considered to have a significant double-excitation character. We report %$T_1$ values computed at the CC3/aug-cc-pVTZ level for transitions belonging to molecules containing up to five non-hydrogen atoms and at the CC3/aug-cc-pVDZ level for transitions of larger molecules.

The various excited states determined at the reference Q-CC3/aug-cc-pVTZ and Q-CC3/aug-cc-pVDZ levels were identified through the literature.⁸¹,⁸⁴ For lower-order methods, assignments were made using standard criteria, i.e., spatial symmetries, energy differences, composition of the states in terms of orbitals, and oscillator strengths. In some instances, states with the same symmetry and close energies exhibit strongly mixed character, rendering their assignment challenging. Given that experimental measurements typically cannot distinguish states with trifling energy differences, in a few cases we aggregated such contributions and treated them as a single state in the statistical analysis. The few states for which unambiguous identification was beyond reach were excluded from the statistical analysis (see the SI). The molecules with non-Abelian group symmetries, namely, acetylene, benzene, and triazine, are treated by quantum chemistry software using the closest Abelian group symmetry. As a result, such states are split in two degenerate states, which are Abelians. Therefore,

only for such cases, for an appropriate comparison with experiments, the $\delta^{TPA}$ reported should be multiplied by 2. A basis set analysis is presented below for Q-CC3 calculations, comparing the $\delta^{TPA}$ obtained with the d-aug-cc-pVTZ (for a few cases), aug-cc-pVTZ, aug-cc-pVDZ, and 6-31+G(d) basis sets for molecules containing from one to five non-hydrogen atoms. For molecules with six non-hydrogen atoms, aug-cc-pVTZ becomes too computationally demanding, and a comparison between aug-cc-pVDZ and 6-31+G(d) is provided.

The aug-cc-pVDZ basis set was employed for the benchmarks, coupled with the corresponding auxiliary basis sets when necessary.⁸⁶ States with $\delta^{TPA}_{Q-CC3} < 3$ au are considered irrelevant for the purposes of this paper and were not 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 this reason, acetylene and streptocyanine-1 have been excluded from the statistical analysis; indeed the $1A_u$ and $2A_u$ states of the former are TPA-forbidden, while the $1B_2$ state of the latter (the only state studied) has a $\delta^{TPA}$ value below the selected threshold. All CC calculations uses the frozen-core approximation with the length gauge for the interaction operator. Standard convergence criteria ($10^{-7}$ au) were 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 thresholds: $10^{-9}$ au for ground-state energies, $10^{-8}$ au for CC equations, and $10^{-7}$ au for EOM Davidson convergence. The CD of two-electron integrals used the standard threshold, $10^{-3}$ au, as it was observed to be sufficient for the TPA strengths.³⁸,⁴⁴,⁶⁴ DEA-CCSD calculations were performed on systems with a significant double-excitation character using Q-Chem with the same thresholds as for its EOM-CCSD version (indeed, the RI and CD approximations are not yet implemented for this method) and employing the charge +2 for reproducing the double excitations.⁴⁷,⁸⁷ I-ADC2 calculations were performed with Q-Chem using default convergence criteria. The coupled
<table><caption>Table 1. Q-CC3 Excitation Energies (ΔE) and Two-Photon Transition Strengths (δTPA) of Molecules with One or Two Non-hydrogen Atoms Computed with the d-aug-cc-pVTZ (da-TZ) and aug-cc-pVTZ (a-TZ) Basis Sets; Values of %T₁ at the Q-CC3/a-TZ Level Are Also Reported</caption><thead><tr><th rowspan="2">molecule</th><th rowspan="2">state</th><th rowspan="2">type</th><th rowspan="2">%T<sub>1</sub></th><th colspan="2">ΔE [eV]</th><th colspan="2">δ<sup>TPA</sup> [au]</th></tr><tr><th>da-TZ</th><th>a-TZ</th><th>da-TZ</th><th>a-TZ</th></tr></thead><tbody><tr><td>hydrogen sulfide</td><td>1B<sub>1</sub></td><td>Ryd (n → 3s)</td><td>94.3</td><td>6.2</td><td>6.2</td><td>6.3</td><td>6.2</td></tr><tr><td rowspan="2">water</td><td>1B<sub>1</sub></td><td>Ryd (n → 3s)</td><td>93.4</td><td>7.6</td><td>7.6</td><td>4.8</td><td>4.8</td></tr><tr><td>1A<sub>2</sub></td><td>Ryd (n → 3p)</td><td>93.6</td><td>9.4</td><td>9.4</td><td>42.2</td><td>44.9</td></tr><tr><td rowspan="4">ethylene<br>formaldehyde</td><td>1B<sub>1g</sub></td><td>Ryd (π → 3p)</td><td>95.3</td><td>8.0</td><td>8.0</td><td>74.3</td><td>69.5</td></tr><tr><td>1B<sub>2</sub></td><td>Ryd (n → 3s)</td><td>91.7</td><td>7.2</td><td>7.2</td><td>80.8</td><td>81.9</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (n → 3p)</td><td>92.4</td><td>8.0</td><td>8.1</td><td>4.0</td><td>4.3</td></tr><tr><td>1A<sub>1</sub></td><td>Ryd (n → 3p)</td><td>91.9</td><td>8.1</td><td>8.2</td><td>90.9</td><td>107.5</td></tr><tr><td></td><td>2A<sub>2</sub></td><td>Ryd (n → 3p)</td><td>91.7</td><td>8.4</td><td>8.6</td><td>3.8</td><td>4.9</td></tr></tbody></table>

cluster implementations employed in Dalton and Q-Chem are non-Hermitian; therefore, the 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 the left and right moments in eq 2 are negative and represent at least 20% of the total value of δTPA were excluded from the analysis. The excluded cases are documented in Table S1. We stress that I-ADC2 and TD-DFT formulations, originating from Hermitian theories, are free of this difficulty. Finally, states expected to exhibit resonance-enhanced TPA effects were also excluded and are detailed in Table S2. 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.

Five XCFs were tested, with the aim of evaluating the impact of the exact exchange percentage: we tested two global hybrid functionals, namely, B3LYP⁸⁹ and BH&HLYP,⁹⁰ with 20% and 50% of exact exchange, respectively, and three range-separated functionals, namely, CAM-B3LYP,⁷⁵ which has exact exchange going from 19% to 65% and a range separation parameter ω of 0.33 Bohr⁻¹, and two different versions of the LC-BLYP functional⁹¹ for which the exact exchange ranges from 0 to 100% and the range separation parameter is 0.47 or 0.33 Bohr⁻¹.

The statistical indicators employed below are the mean average percentage error (MAPE), the mean signed error (MSE), the mean absolute error (MAE), the root-mean-square deviation (RMSD), and the maximum positive and negative errors (MAX⁺ and MAX⁻), given by eqs 3–8:

$$ \text{MAPE} = 100 \cdot \frac{1}{n} \sum_{i=1}^{n} \left| \frac{x_i - t_i}{t_i} \right| \quad (3) $$

$$ \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (x_i - t_i) \quad (4) $$

$$ \text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |x_i - t_i| \quad (5) $$

$$ \text{RMSD} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (x_i - t_i)^2} \quad (6) $$

$$ \text{MAX}_{+} = \max_{i} \{x_i - t_i | x_i - t_i > 0\} \quad (7) $$

$$ \text{MAX}_{-} = \max_{i} \{t_i - x_i | t_i - x_i > 0\} \quad (8) $$

in which $x_i$ is the value obtained with the low-level method, $t_i$ is the reference value, and $n$ is the number of values. The indicator of linear correlation, R², and the standard deviation of errors (SDE) are also reported below.

### 3. RESULTS AND DISCUSSION

This section is structured as follows: initially, the Q-CC3 results are presented, and a comprehensive analysis is conducted to benchmark the impact of the basis set on δTPA. Subsequently, the accuracy of the lower-order methods is assessed across the entire set of states. This assessment encompasses evolution of the method performance across diverse δTPA ranges, considering the molecular dimensions and the nature of the states, namely, valence, Rydberg, and states showing significant double-excitation character. In the following, we discuss the δTPA. The interested reader will find in Tables S21–S33 similar statistical analyses for the excitation energies. The quality of these energies for various levels of theory was discussed in previous papers on the QUEST database.⁸¹⁻⁸⁵

#### 3.1. Reference Values and Basis Set Analysis.
As a preliminary analysis, we assess the reliability of the aug-cc-pVTZ basis set relative to its doubly augmented counterpart, d-aug-cc-pVTZ, for selected small molecules: hydrogen sulfide, water, ethylene, and formaldehyde. In fact, previous studies have shown that for these small molecules, a doubly augmented basis set might be necessary to reach basis set convergence for δTPA. Due to its highly diffuse nature, the d-aug-cc-pVTZ basis set can introduce significant state mixing, hindering an unambiguous identification of the states and hence sound comparison with aug-cc-pVTZ results. We observe this issue for the following states: 1A₁ of water, 5A₉ of ethylene, and 2A₁ of formaldehyde. Therefore, we report in Table 1 only the results for states in which comparisons between the two basis sets are not plagued by state-mixing issues.

For both hydrogen sulfide and water, the aug-cc-pVTZ basis set provides δTPA values matching the d-aug-cc-pVTZ ones, with differences ranging from 0 to 2.7 au. For ethylene, we observe a slight underestimation of 5 au for the 1B₁g state when using aug-cc-pVTZ. For formaldehyde, while the 1B₂, 2B₂, and 2A₂ states are well-described with aug-cc-pVTZ, with a maximum error of 1 au, the 1A₁ state shows a significant
<table><caption>Table 2. Q-CC3 Excitation Energies (ΔE) and Two-Photon Transition Strengths (δTPA) of Small Molecules Computed with the aug-cc-pVTZ (a-TZ), aug-cc-pVDZ (a-DZ), and 6-31+G(d) (p-DZ) Basis Sets; Values of %T₁ at the Q-CC3/a-TZ Level Are Also Reported</caption><thead><tr><th rowspan="2">molecule</th><th rowspan="2">state</th><th rowspan="2">type</th><th rowspan="2">%T1</th><th colspan="3">ΔE [eV]</th><th colspan="3">δ<sup>TPA</sup> [au]</th></tr><tr><th>a-TZ</th><th>a-DZ</th><th>p-DZ</th><th>a-TZ</th><th>a-DZ</th><th>p-DZ</th></tr></thead><tbody><tr><td rowspan="2">hydrogen sulfide</td><td>1B<sub>1</sub></td><td>Ryd (n → 3s)</td><td>94.3</td><td>6.2</td><td>6.1</td><td>6.0</td><td>6.2</td><td>6.9</td><td>6.5</td></tr><tr><td>1B<sub>1</sub></td><td>Ryd (n → 3s)</td><td>93.4</td><td>7.6</td><td>7.5</td><td>8.3</td><td>4.8</td><td>4.7</td><td>3.9</td></tr><tr><td rowspan="2">water</td><td>1A<sub>2</sub></td><td>Ryd (n → 3p)</td><td>93.6</td><td>9.4</td><td>9.3</td><td>10.7</td><td>44.9</td><td>46.5</td><td>23.0</td></tr><tr><td>1A<sub>1</sub></td><td>Ryd (n → 3s)</td><td>93.6</td><td>10.0</td><td>9.9</td><td>10.8</td><td>9.1</td><td>8.8</td><td>8.3</td></tr><tr><td rowspan="2">ethylene</td><td>1B<sub>1g</sub></td><td>Ryd (π → 3p)</td><td>95.3</td><td>8.0</td><td>8.0</td><td>8.3</td><td>69.5</td><td>67.2</td><td>76.4</td></tr><tr><td>5A<sub>g</sub></td><td>Val (double)</td><td>20.3</td><td>13.4</td><td>13.6</td><td>13.8</td><td>894.1</td><td>1636.7</td><td>1947.9</td></tr><tr><td rowspan="5">formaldehyde</td><td>1B<sub>2</sub></td><td>Ryd (n → 3s)</td><td>91.7</td><td>7.2</td><td>7.1</td><td>7.2</td><td>81.9</td><td>91.5</td><td>63.9</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (n → 3p)</td><td>92.4</td><td>8.1</td><td>8.0</td><td>8.0</td><td>4.3</td><td>5.2</td><td>3.9</td></tr><tr><td>1A<sub>1</sub></td><td>Ryd (n → 3p)</td><td>91.9</td><td>8.2</td><td>8.1</td><td>8.3</td><td>107.5</td><td>116.4</td><td>54.5</td></tr><tr><td>2A<sub>2</sub></td><td>Ryd (n → 3p)</td><td>91.7</td><td>8.6</td><td>8.7</td><td>8.7</td><td>4.9</td><td>5.6</td><td>3.8</td></tr><tr><td>2A<sub>1</sub></td><td>Val (π → π*)</td><td>90.4</td><td>9.5</td><td>9.6</td><td>9.8</td><td>15.4</td><td>10.6</td><td>5.9</td></tr><tr><td rowspan="2">cyclopropene</td><td>1B<sub>2</sub></td><td>Val (π → π*)</td><td>95.0</td><td>6.7</td><td>6.8</td><td>7.0</td><td>10.7</td><td>12.5</td><td>9.8</td></tr><tr><td>2B<sub>1</sub></td><td>Ryd (π → 3s)</td><td>95.2</td><td>6.9</td><td>6.9</td><td>7.2</td><td>88.0</td><td>94.8</td><td>59.5</td></tr><tr><td rowspan="4">acetone</td><td>1B<sub>2</sub></td><td>Ryd (n → 3s)</td><td>90.5</td><td>6.4</td><td>6.3</td><td>6.6</td><td>4.7</td><td>5.9</td><td>4.2</td></tr><tr><td>2A<sub>2</sub></td><td>Ryd (n → 3p)</td><td>90.8</td><td>7.4</td><td>7.4</td><td>7.8</td><td>51.1</td><td>54.3</td><td>63.7</td></tr><tr><td>1A<sub>1</sub></td><td>Ryd (n → 3p)</td><td>90.6</td><td>7.5</td><td>7.4</td><td>7.8</td><td>648.7</td><td>725.1</td><td>793.9</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (n → 3p)</td><td>91.2</td><td>7.6</td><td>7.6</td><td>7.9</td><td>148.5</td><td>164.0</td><td>181.5</td></tr><tr><td rowspan="4">acrolein</td><td>1A<sub>′</sub></td><td>Val (π → π*)</td><td>91.6</td><td>6.7</td><td>6.7</td><td>6.8</td><td>88.1</td><td>105.8</td><td>118.1</td></tr><tr><td>2A<sub>′</sub></td><td>Val (n → π*)</td><td>79.4</td><td>6.8</td><td>6.8</td><td>6.9</td><td>5.0</td><td>5.3</td><td>3.8</td></tr><tr><td>2A<sub>′</sub></td><td>Ryd (n → 3s)</td><td>89.9</td><td>7.2</td><td>7.0</td><td>7.2</td><td>99.8</td><td>104.0</td><td>91.2</td></tr><tr><td>5A<sub>′</sub></td><td>Val (π → π*)</td><td>75.6</td><td>8.1</td><td>8.1</td><td>8.2</td><td>172.9</td><td>168.1</td><td>142.5</td></tr><tr><td rowspan="4">butadiene</td><td>1B<sub>g</sub></td><td>Ryd (π → 3s)</td><td>94.1</td><td>6.3</td><td>6.3</td><td>6.5</td><td>153.0</td><td>158.8</td><td>127.5</td></tr><tr><td>1A<sub>g</sub></td><td>Val (π → π*)</td><td>75.1</td><td>6.7</td><td>6.7</td><td>6.7</td><td>372.4</td><td>389.4</td><td>348.1</td></tr><tr><td>2B<sub>g</sub></td><td>Ryd (π → 3d)</td><td>94.2</td><td>7.5</td><td>7.5</td><td>7.7</td><td>134.3</td><td>163.0</td><td>163.6</td></tr><tr><td>3B<sub>g</sub></td><td>Ryd (π → 3d)</td><td>94.5</td><td>7.6</td><td>7.6</td><td>8.0</td><td>64.8</td><td>72.2</td><td>48.1</td></tr><tr><td rowspan="2">glyoxal</td><td>2B<sub>g</sub></td><td>Val (n → π*)</td><td>83.9</td><td>6.6</td><td>6.6</td><td>6.8</td><td>9.9</td><td>10.1</td><td>7.0</td></tr><tr><td>1B<sub>2</sub></td><td>Val (π → π*)</td><td>85.4</td><td>4.3</td><td>4.3</td><td>4.4</td><td>27.8</td><td>28.8</td><td>23.2</td></tr><tr><td rowspan="4">methylenecyclopropene</td><td>1B<sub>1</sub></td><td>Ryd (π → 3s)</td><td>93.6</td><td>5.4</td><td>5.4</td><td>5.7</td><td>218.3</td><td>226.8</td><td>202.8</td></tr><tr><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>92.8</td><td>6.1</td><td>6.1</td><td>6.2</td><td>50.5</td><td>51.9</td><td>52.1</td></tr><tr><td>1A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>94.0</td><td>5.8</td><td>5.7</td><td>6.1</td><td>137.8</td><td>148.8</td><td>99.0</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (π → 3p)</td><td>94.2</td><td>6.6</td><td>6.5</td><td>7.1</td><td>12.4</td><td>18.6</td><td>8.1</td></tr><tr><td rowspan="4">cyclopentadiene</td><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>78.9</td><td>6.6</td><td>6.6</td><td>6.7</td><td>125.9</td><td>126.7</td><td>126.2</td></tr><tr><td>1A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>93.8</td><td>6.1</td><td>6.0</td><td>6.3</td><td>137.6</td><td>146.1</td><td>121.6</td></tr><tr><td>1B<sub>2</sub></td><td>Val (π → π*)</td><td>93.0</td><td>6.3</td><td>6.4</td><td>6.5</td><td>15.0</td><td>18.2</td><td>18.3</td></tr><tr><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>85.0</td><td>6.6</td><td>6.6</td><td>6.7</td><td>74.4</td><td>79.1</td><td>84.8</td></tr><tr><td rowspan="4">furan</td><td>2A<sub>2</sub></td><td>Ryd (π → 3p)</td><td>93.6</td><td>6.8</td><td>6.7</td><td>7.0</td><td>16.2</td><td>22.4</td><td>12.4</td></tr><tr><td>1A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>92.8</td><td>5.2</td><td>5.1</td><td>5.3</td><td>114.5</td><td>125.6</td><td>109.2</td></tr><tr><td>1B<sub>2</sub></td><td>Ryd (mixed)</td><td>92.4</td><td>6.0</td><td>5.9</td><td>6.0</td><td>15.5</td><td>16.5</td><td>8.8</td></tr><tr><td>2A<sub>2</sub></td><td>Ryd (π → 3p)</td><td>93.0</td><td>6.0</td><td>5.9</td><td>6.3</td><td>75.0</td><td>76.0</td><td>54.7</td></tr><tr><td rowspan="4">pyrrole</td><td>2B<sub>1</sub></td><td>Ryd (mixed)</td><td>92.3</td><td>6.1</td><td>6.0</td><td>6.2</td><td>36.7</td><td>32.1</td><td>22.1</td></tr><tr><td>1B<sub>2</sub></td><td>Val (π → π*)</td><td>92.5</td><td>6.3</td><td>6.3</td><td>6.3</td><td>25.1</td><td>29.1</td><td>30.1</td></tr><tr><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>86.3</td><td>6.3</td><td>6.3</td><td>6.4</td><td>86.5</td><td>94.3</td><td>106.0</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (π → 3p)</td><td>92.6</td><td>6.8</td><td>7.0</td><td>7.2</td><td>4.0</td><td>3.7</td><td>6.6</td></tr><tr><td rowspan="5">thiophene</td><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>87.6</td><td>5.7</td><td>5.7</td><td>5.8</td><td>37.6</td><td>43.2</td><td>59.2</td></tr><tr><td>1B<sub>1</sub></td><td>Ryd (π → 3p)</td><td>90.1</td><td>6.1</td><td>6.2</td><td>6.2</td><td>4.5</td><td>2.8</td><td>4.9</td></tr><tr><td>1A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>92.6</td><td>6.1</td><td>6.1</td><td>6.3<sup>a</sup></td><td>110.9</td><td>115.8</td><td>68.3<sup>a</sup></td></tr><tr><td>2A<sub>2</sub></td><td>Ryd (π → 3p)</td><td>91.8</td><td>6.3</td><td>6.3</td><td>6.3<sup>a</sup></td><td>11.5</td><td>9.8</td><td>44.7<sup>a</sup></td></tr><tr><td>2B<sub>1</sub></td><td>Ryd (π → 3s)</td><td>92.8</td><td>6.5</td><td>6.4</td><td>6.6</td><td>117.4</td><td>116.2</td><td>117.6</td></tr><tr><td rowspan="2">streptocyanine-3</td><td>2A<sub>1</sub></td><td>Val (π → π*)</td><td>86.5</td><td>7.4</td><td>7.4</td><td>7.5</td><td>30.5</td><td>37.0</td><td>39.6</td></tr><tr><td>1B<sub>2</sub></td><td>Val (π → π*)</td><td>87.2</td><td>4.8</td><td>4.8</td><td>4.8</td><td>22.8</td><td>23.0</td><td>31.1</td></tr></tbody></table>

overestimation of 17.5 au (18%), suggesting that d-aug-cc-pVTZ would be needed to obtain basis set convergence in that case.

Next, we made an extensive basis set analysis for the other tested bases, starting by analyzing compounds containing one

to five non-hydrogen atoms. Excitation energies and δTPA values for the excited states in which δTPA > 3 au at the Q-CC3/aug-cc-pVTZ level are listed in Table 2, while states with smaller response are given in Table S3 in the SI. The 50 states considered in this analysis encompass 18 valence states (16 π

<footnote>a Mixed states (summed in the statistical analysis).</footnote>
→ π* transitions and two n → π* transitions), among which six exhibit significant double-excitation character, each with different contributions of the single excitations in their composition: 5A<sub>g</sub> of ethylene (20.0%), 2A<sup>′</sup> (79.4%) and 5A<sup>′</sup> (75.0%) of acrolein, 1A<sub>g</sub> of butadiene (75.1%), and 1A<sub>1</sub> of cyclopentadiene (78.9%); the values in parentheses are the %T<sub>1</sub> values. The 5A<sub>g</sub> state of ethylene, characterized by a high degree of double excitation (80.0%), is anticipated to be largely affected not only by the basis set, as previously demonstrated for energies by one of us,<sup>82,92</sup> but more importantly by the level of theory. Indeed, CC3 is far from an ideal method for such genuine double excitation. The δ<sup>TPA</sup> values of that state are 894.1, 1636.7, and 1947.9 au with aug-cc-pVTZ, aug-cc-pVDZ, and 6-31+G(d), respectively. These substantial differences led us to conclude that small basis sets are inappropriate for accurately describing this state, prompting us to exclude this transition from further statistical analysis. The other four states display significantly smaller contributions of double excitations, and therefore, they are expected to be relatively less challenging, with Q-CC3 likely being well-suited.

The remaining 32 states have Rydberg character. Table 2 illustrates that for the thiophene molecule, the 1A<sub>2</sub> and 2A<sub>2</sub> states look strongly mixed with the 6-31+G(d) basis set, whereas the situation is clearer with larger basis sets. Indeed, the TPA values are 110.9 au (1A<sub>2</sub>) and 11.5 au (2A<sub>2</sub>) with aug-cc-pVTZ compared with 115.8 au (1A<sub>2</sub>) and 9.8 au (2A<sub>2</sub>) with aug-cc-pVDZ. The Q-CC3/6-31+G(d) method is unable to distinguish these two states, separated by 0.06 eV only, resulting in a splitting of the δ<sup>TPA</sup> values (68.3 and 44.7 au). As previously stated, we consequently summed these two contributions and considered them as a single state in this case.

The statistical analysis performed on the resulting 48 cases is presented in Table 3, using the aug-cc-pVTZ values as a

<table><caption>Table 3. Statistical Analysis of δ<sup>TPA</sup> Values for Small Molecules (See Table 2 for Details); Reference Values Were Computed at the Q-CC3/aug-cc-pVDZ Level</caption><thead><tr><th rowspan="2">indicator</th><th colspan="2">all</th><th colspan="2">Rydberg</th><th colspan="2">valence</th></tr><tr><th>a-DZ</th><th>p-DZ</th><th>a-DZ</th><th>p-DZ</th><th>a-DZ</th><th>p-DZ</th></tr></thead><tbody><tr><td>n</td><td>48</td><td>48</td><td>31</td><td>31</td><td>17</td><td>17</td></tr><tr><td>MAPE [%]</td><td>11.5</td><td>22.5</td><td>11.9</td><td>21.9</td><td>10.7</td><td>23.7</td></tr><tr><td>MSE [au]</td><td>-5.9</td><td>0.9</td><td>-7.0</td><td>2.6</td><td>-3.7</td><td>-2.1</td></tr><tr><td>MAE [au]</td><td>6.7</td><td>15.1</td><td>7.7</td><td>17.5</td><td>4.8</td><td>10.8</td></tr><tr><td>RMSD [au]</td><td>13.4</td><td>27.1</td><td>15.8</td><td>31.9</td><td>7.0</td><td>14.7</td></tr><tr><td>MAX<sub>+</sub> [au]</td><td>4.8</td><td>53.1</td><td>4.6</td><td>53.1</td><td>4.8</td><td>30.4</td></tr><tr><td>MAX<sub>-</sub> [au]</td><td>76.5</td><td>145.2</td><td>76.5</td><td>145.2</td><td>17.7</td><td>30.0</td></tr><tr><td>SDE [au]</td><td>12.2</td><td>27.3</td><td>14.4</td><td>32.3</td><td>6.2</td><td>15.0</td></tr><tr><td>R<sup>2</sup></td><td>1.00</td><td>0.96</td><td>1.00</td><td>0.97</td><td>1.00</td><td>0.98</td></tr></tbody></table>

reference. According to the reported statistical indicators, the aug-cc-pVDZ results exhibit a MAPE of 11.5% and a MAE of 6.7 au, indicating robust reproduction of the reference data. In contrast, when employing the smaller 6-31+G(d) basis set we observe a stronger degradation of the performances in both the MAPE (22.5%) and MAE (15.1 au). It is noteworthy that all error indicators of aug-cc-pVDZ results are approximately half of their 6-31+G(d) counterparts. Similar conclusions can be drawn when analyzing valence and Rydberg states separately, as also illustrated in Table 3.

For states with six non-hydrogen atoms, Q-CC3/aug-cc-pVTZ becomes too computationally demanding. Encouraged by the performances of aug-cc-pVDZ for smaller compounds,

we took aug-cc-pVDZ data as references to evaluate the performance of 6-31+G(d) when larger compounds are considered. The δ<sup>TPA</sup> values and energies for 34 excited states showing δ<sup>TPA</sup> > 3 au with aug-cc-pVDZ are reported in Table 4, while states with smaller response values are reported in Table S3 of the SI. This set encompasses 13 valence states (seven π → π* and six n → π* transitions). The 1A<sub>g</sub> state of hexatriene and the 1E<sub>2g</sub> state of benzene exhibit a substantial double character with %T<sub>1</sub> values of 66.5% and 67.8%, respectively. Within this set, there are also 20 Rydberg states and one mixed state, the 3A<sub>1</sub> state of pyridazine. We evaluated the performance of 6-31+G(d) with respect to aug-cc-pVDZ in 82 states: 48 states collected in Table 2 and 34 states in Table 4. The results are listed in Table 5. For the full set, we found that δ<sup>TPA</sup> with the 6-31+G(d) basis set significantly differs from the aug-cc-pVDZ values with a MAPE of 24.7% and a MAE of 15.1 au. These values are essentially unchanged with respect to those in Table 3. The large MAX<sub>+</sub> in Table 5 comes from the 1A<sub>g</sub> state of hexatriene, with an overestimation that might be ascribed to the worse description of double excitations with the compact basis. By excluding the six states where double excitations are not negligible (i.e., %T<sub>1</sub> < 80%), we observe a slight increase of the MAPE but a large improvement of the other indicators. Table 5 also indicates generally better performance of 6-31+G(d) for the valence states, although the indicators stay in the same range for the Rydberg transitions. In Table S20, we present a statistical analysis where the systems are divided into two groups based on their size: those with one to four non-hydrogen atoms and those with five or six non-hydrogen atoms. For all single excitations, we observe a slight increase in both MAPE and MAE as the system size increases, with MAPE rising from 22.0% to 25.8% and MAE increasing from 14.7 to 15.3 au. This is likely due to the fact that valence states are better described in smaller molecules compared to larger ones, with MAPE values of 21.3% and 33.1%, respectively. For Rydberg states, however, the difference between the two groups appears to be minimal. Further comment should be made on the size of δ<sup>TPA</sup> and its implications for practical applications. The studied molecules exhibit rather small values of δ<sup>TPA</sup> and of the corresponding σ<sup>TPA</sup>. For instance, one of the largest δ<sup>TPA</sup> values was found for the 1A<sub>g</sub> transition in hexatriene (1358.3 au), which translates to a σ<sup>TPA</sup> value of 17 GM (assuming Lorentzian broadening of 0.2 eV). This value is an order of magnitude smaller than the values reported for push-pull systems such as 4-dimethylamino-4'-nitrostilbene.<sup>93</sup> Despite the small values, studying such small responses is crucial for evaluating the performance of QM methods and understanding their predictive power and limitations.

## 3.2. Benchmark of Lower-Order Methods. 3.2.1. General Trends.
In Figure 1, we provide error distribution patterns for all methods included in this study. The values of MAE, MSE, and SDE are given in atomic units, as are all values in the following. Obviously, CCSD methods significantly outperform the other methods in all of the indicators considered independently if the Q or E formalism is employed. Q-CCSD is the most accurate method, with a MAE of 5.7, an MSE of 0.1, and an SDE of 15.6. Similarly, E-CCSD reproduces Q-CC3 references well, though with slightly deteriorated accuracy compared to the more costly Q-CCSD, with a MAE of 7.9, an MSE of 3.9, and an SDE of 17.7; these differences can be mainly ascribed to the valence states (see below). Conversely, Q-CC2 and I-ADC2 exhibit similar but
<table><caption>Table 4. Q-CC3 Excitation Energies (ΔE) and Two-Photon Transition Strengths (δTPA) of Molecules with Six Non-hydrogen Atoms Computed with the aug-cc-pVDZ (a-DZ) and 6-31 + G(d) (p-DZ) Basis Sets; Values of %T₁ at the Q-CC3/a-DZ Level Are Also Reportedᵃ</caption><thead><tr><th rowspan="2">molecule</th><th rowspan="2">state</th><th rowspan="2">type</th><th rowspan="2">%T<sub>1</sub></th><th colspan="2">ΔE [eV]</th><th colspan="2">δ<sup>TPA</sup> [au]</th></tr><tr><th>a-DZ</th><th>p-DZ</th><th>a-DZ</th><th>p-DZ</th></tr></thead><tbody><tr><td rowspan="2">benzene*</td><td>1E<sub>1g</sub></td><td>Ryd (π → 3s)</td><td>92.8</td><td>6.5</td><td>6.7</td><td>87.5</td><td>73.7</td></tr><tr><td>1E<sub>2g</sub></td><td>Val (π → π*)</td><td>67.8</td><td>8.4</td><td>8.5</td><td>30.7</td><td>37.0</td></tr><tr><td rowspan="2">hexatriene</td><td>1A<sub>g</sub></td><td>Val (π → π*)</td><td>66.5</td><td>5.8</td><td>5.8</td><td>1385.3</td><td>1078.6</td></tr><tr><td>1B<sub>g</sub></td><td>Ryd (π → 3s)</td><td>93.4</td><td>5.8</td><td>6.1</td><td>9.8</td><td>8.2</td></tr><tr><td rowspan="6">pyrazine</td><td>1B<sub>2g</sub></td><td>Val (n → π*)</td><td>85.3</td><td>5.8</td><td>5.9</td><td>7.1</td><td>8.5</td></tr><tr><td>1A<sub>g</sub></td><td>Ryd (n → 3s)</td><td>91.1</td><td>6.5</td><td>6.7</td><td>86.6</td><td>63.3</td></tr><tr><td>2B<sub>3g</sub></td><td>Ryd (π → 3s)</td><td>93.8</td><td>7.1</td><td>7.4</td><td>149.9</td><td>150.2</td></tr><tr><td>2A<sub>g</sub></td><td>Ryd (n.d.)</td><td>90.7</td><td>7.9</td><td>8.1</td><td>211.1</td><td>173.1</td></tr><tr><td>1B<sub>1g</sub></td><td>Ryd (n.d.)</td><td>91.1</td><td>7.9</td><td>8.2</td><td>20.6</td><td>16.4</td></tr><tr><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>85.8</td><td>5.3</td><td>5.4</td><td>3.2</td><td>3.4</td></tr><tr><td rowspan="6">pyridazine</td><td>1B<sub>2</sub></td><td>Ryd (n → 3s)</td><td>88.4</td><td>6.1</td><td>6.3</td><td>39.0</td><td>42.2</td></tr><tr><td>2A<sub>1</sub></td><td>Ryd (n → 3p)</td><td>88.6</td><td>6.7</td><td>6.9</td><td>11.5</td><td>8.5</td></tr><tr><td>2B<sub>2</sub></td><td>Val (π → π*)</td><td>90.5</td><td>6.8</td><td>7.0</td><td>4.9</td><td>3.6</td></tr><tr><td>3A<sub>2</sub></td><td>Val (n → π*)</td><td>88.6</td><td>7.4</td><td>7.6</td><td>16.2</td><td>19.8</td></tr><tr><td>3A<sub>1</sub></td><td>Mixed (π → π*)</td><td>89.2</td><td>7.5</td><td>7.7</td><td>488.8</td><td>448.6</td></tr><tr><td>1A<sub>1</sub></td><td>Ryd (n → 3s)</td><td>90.2</td><td>6.6</td><td>6.8</td><td>49.3</td><td>29.3</td></tr><tr><td rowspan="6">pyridine</td><td>2A<sub>1</sub></td><td>Val (π → π*)</td><td>91.7</td><td>6.7</td><td>6.9</td><td>8.2</td><td>23.1</td></tr><tr><td>2A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>93.2</td><td>6.7</td><td>7.0</td><td>90.1</td><td>82.2</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (n → 3p)</td><td>89.5</td><td>7.3</td><td>7.5</td><td>14.5</td><td>11.3</td></tr><tr><td>2B<sub>1</sub></td><td>Ryd (π → 3p)</td><td>93.5</td><td>7.3</td><td>7.5</td><td>17.2</td><td>18.2</td></tr><tr><td>3A<sub>1</sub></td><td>Ryd (n → 3s)</td><td>89.8</td><td>7.3</td><td>7.6</td><td>36.8</td><td>39.4</td></tr><tr><td>3B<sub>1</sub></td><td>Ryd (π → 3s)</td><td>91.9</td><td>7.4</td><td>7.7</td><td>28.2</td><td>22.4</td></tr><tr><td rowspan="4">pyrimidine</td><td>3A<sub>2</sub></td><td>Ryd (π → 3s)</td><td>93.0</td><td>7.6</td><td>7.8</td><td>38.5</td><td>35.7</td></tr><tr><td>4B<sub>1</sub></td><td>Val (n → π*)</td><td>89.7</td><td>8.0</td><td>8.1</td><td>6.0</td><td>6.6</td></tr><tr><td>2A<sub>2</sub></td><td>Val (n → π*)</td><td>86.8</td><td>6.0</td><td>6.1</td><td>3.2</td><td>4.0</td></tr><tr><td>2B<sub>2</sub></td><td>Ryd (n → 3s)</td><td>90.2</td><td>6.6</td><td>6.8</td><td>33.5</td><td>29.0</td></tr><tr><td rowspan="4">tetrazine</td><td>1A<sub>1</sub></td><td>Val (π → π*)</td><td>91.6</td><td>6.9</td><td>7.1</td><td>3.0</td><td>5.0</td></tr><tr><td>2A<sub>1</sub></td><td>Ryd (n → 3p)</td><td>90.5</td><td>7.3</td><td>7.6</td><td>10.5</td><td>13.0</td></tr><tr><td>3B<sub>1</sub></td><td>Ryd (π → 3s)</td><td>93.5</td><td>7.4</td><td>7.6</td><td>66.7</td><td>58.9</td></tr><tr><td>1B<sub>3g</sub></td><td>Val (n → π*)</td><td>81.3</td><td>5.5</td><td>5.6</td><td>56.0</td><td>87.3</td></tr><tr><td rowspan="3">triazine*</td><td>1B<sub>2g</sub></td><td>Ryd (n → 3s)</td><td>91.0</td><td>6.3</td><td>6.5</td><td>33.2</td><td>24.8</td></tr><tr><td>1E'</td><td>Ryd (n → 3s)</td><td>90.8</td><td>7.2</td><td>7.4</td><td>36.6</td><td>28.5</td></tr><tr><td>2E"</td><td>Val (n → π*)</td><td>82.6</td><td>7.8</td><td>7.9</td><td>11.7</td><td>12.9</td></tr><tr><td></td><td>2E'</td><td>Val (π → π*)</td><td>89.5</td><td>8.0</td><td>8.1</td><td>5.0</td><td>1.3</td></tr></tbody></table>

<table><caption>Table 5. Statistical analysis of δTPA computed at 6-31+G(d) for the states in Table 2 and 4. Reference values are computed at the Q-CC3/aug-cc-pVDZ level</caption><thead><tr><th>indicator</th><th>all</th><th>all singly excited</th><th>Rydberg</th><th>valence</th></tr></thead><tbody><tr><td>n</td><td>82</td><td>76</td><td>51</td><td>24</td></tr><tr><td>MAPE [%]</td><td>24.7</td><td>25.3</td><td>23.4</td><td>30.2</td></tr><tr><td>MSE [au]</td><td>-9.2</td><td>-5.1</td><td>-8.6</td><td>3.8</td></tr><tr><td>MAE [au]</td><td>15.1</td><td>11.3</td><td>13.4</td><td>5.6</td></tr><tr><td>RMSD [au]</td><td>38.5</td><td>18.1</td><td>20.4</td><td>9.0</td></tr><tr><td>MAX<sub>+</sub> [au]</td><td>306.7</td><td>62.0</td><td>62.0</td><td>5.6</td></tr><tr><td>MAX<sub>-</sub> [au]</td><td>68.7</td><td>68.7</td><td>68.7</td><td>31.2</td></tr><tr><td>SDE [au]</td><td>37.6</td><td>17.4</td><td>18.7</td><td>8.3</td></tr><tr><td>R<sup>2</sup></td><td>0.97</td><td>0.97</td><td>0.98</td><td>0.96</td></tr></tbody></table>

significantly inferior performances compared with CCSD, with MAEs of 21.5 (22.8) for Q-CC2 (I-ADC2), over 4 times larger than the Q-CCSD value. Moreover, they demonstrate a quite large overestimation of δTPA, resulting in MSEs of 20.2 (18.8) and large error spreads with SDEs of 49.7 (34.4) for Q-CC2

(I-ADC2). Among the tested TD-DFT methods, we observe, as expected, significantly different results based on the type of functional and the amount of exact exchange included. The most accurate methods are the BH&HLYP and CAM-B3LYP functionals, which yield similar results with respective MAEs of 15.8 and 14.6, MSEs of -2.9 and -3.9 (indicating a slight underestimation), and a significant spread of the errors with SDEs of 30.1 and 32.2. Less accurate are the two LC functionals, which show similar trends. LC-BLYP33 exhibits a MAE of 19.1, an MSE of -6.8, and a large SDE of 42.2. LC-BLYP47, while having similar values for MAE (21.9) and SDE (47.4), exhibits an MSE of -14.5, more than twice as large as its ω = 0.33 Bohr⁻¹ counterpart. B3LYP is the worst-performing method among those tested for all indicators, with a MAE of 25.5, an MSE of 16.9, and an SDE of 54.1. Interestingly, all TD-DFT methods (excluding B3LYP) outperform or slightly underperform the Q-CC2 and I-ADC2 wavefunction methods for all indicators for a smaller computational cost. A similar trend is observed for the excitation energies (see Tables S22–S24), with functionals like

<footnote>ᵃMolecules labeled with * have non-Abelian group symmetry, and their δTPA values must be multiplied by two when compared with experiments (see section 2 for more details).</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 1.** Error distributions, MAEs, MSEs, and SDEs in atomic units with respect to reference Q-CC3/aug-cc-pVDZ $\delta^{TPA}$ for states without a large contribution of double excitations.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 2.** MAPE and $R^2$ of $\delta^{TPA}$ for all of the wavefunction and TD-DFT methods for the sets with all singly excited, Rydberg, and valence states. References are the $\delta^{TPA}_{Q-CC3}/\text{aug-cc-pVDZ}$ calculations with a threshold of 3 au. The values of the indicators are reported in Tables S8–S10.</figcaption></figure>

CAM-B3LYP and LC-BLYP33 slightly outperforming both CC2 and ADC2, which tend to underestimate the reference values. This indicates that benchmarking TD-DFT on the basis of Q-CC2 results might not always be a fully trustworthy procedure.

In Figure 2, the MAPE and $R^2$ indicators are reported, including separate analyses for valence and Rydberg transitions. Looking at the MAPE for all of the singly excited states, one can see excellent results for Q-CCSD (15.2%), and E-CCSD (20.2%). The other wavefunction methods deliver worse performances, resulting in MAPEs of 43.9% for Q-CC2 and 57.1% for I-ADC2. TD-DFT again outperforms the second-order methods but remains less accurate than the CCSD variants, following the order CAM-B3LYP (34.6%), BH&HLYP (39.3%), LC-BLYP33 (42.7%), LC-BLYP47 (43.2%), and B3LYP (55.6%). Looking at $R^2$, we observe that all of the methods have excellent $R^2$ values for the whole set of molecules, ranging from 0.99 for Q-CCSD to 0.87 for

LC-BLYP33. Interestingly, the methods that underperform on absolute and percentage errors, namely, Q-CC2, I-ADC2, and B3LYP, do provide $R^2$ close to 1.

3.2.2. *Subgroups.* To understand whether the size of the systems affects the accuracy of the methods, we performed separate statistical analyses in two groups as already done above. We collected the related indicators in Tables S16–S18. From this analysis, we observe small variations in the trends for each method when moving from smaller to larger molecules. Interestingly, we observe how the performance of the Q-CCSD, E-CCSD, Q-CC2, and I-ADC2 methods decreases. For example, the MAPE of Q-CCSD passes from 6.5% to 18.4%, that of Q-CC2 from 40.4% to 45.1%, and that of I-ADC2 from 46.1% to 60.8% when going from smaller to larger molecules. Similar trends are found for excitation energies at the CCSD level, while CC2 and ADC2 show a slight increase in performance for larger molecules. On the other hand, the performances of TD-DFT methods for both energies and $\delta^{TPA}$
<figure><img src="image_15.png" /><figcaption>**Figure 3.** MAPEs and $R^2$ values of $\delta^{\text{TPA}}$ for all of the wavefunction and TD-DFT methods for the sets of all the singly excited states for different $\delta^{\text{TPA}}$ thresholds (in au). References are at the Q-CC3/aug-cc-pVDZ level. The values of all the indicators are listed in Tables S11–S13.</figcaption></figure>

tend to improve, with for instance the MAPE of CAM-B3LYP going from 38.7% to 32.8%. An exception to this trend is the LC-BLYP47 functional, which shows MAPEs of 36.8% for small molecules and 45.6% for the large ones with a similar trend for energies (Tables S32 and S33).

A further analysis can be made by collecting data according to the different magnitudes of $\delta^{\text{TPA}}$, in particular considering the intervals $t_1 = 3 \le \delta^{\text{TPA}}_{Q-CC3} < 20$, $t_2 = 20 \le \delta^{\text{TPA}}_{Q-CC3} < 80$, and $t_3 = 80 \le \delta^{\text{TPA}}_{Q-CC3}$. In Figure 3 we report the MAPEs and $R^2$ values computed for these three groups, while the rest of the indicators are reported in Tables S11–S13.

Looking at both MAPE and $R^2$, we can observe how the general trends of Figure 2 are replicated for all of the methods at different thresholds. However, it is noteworthy that the MAPE ($R^2$) decreases (increases) with increasing $\delta^{\text{TPA}}_{Q-CC3}$ threshold. Indeed, looking at the Q-CCSD (E-CCSD) results, we can observe that it goes from 25.9% (29.6%) for the $t_1$ interval to 10.4% (16.9%) for the $t_2$ interval and 5.4% (8.1%) for the $t_3$ interval. Similar decreases are obtained for the other methods (see Tables S11 and S13). Analogous trends can be found for $R^2$. For the states in $t_1$, Q-CCSD and E-CCSD are the only ones showing $R^2$ larger than 0.5 (0.84 and 0.82, respectively), while in $t_2$ all methods have a correlation larger than 0.4 (in which Q-CCSD excels with an $R^2$ of 0.85). On the other hand, in the $t_3$ interval all methods show correlation ($R^2 > 0.91$). It is therefore good news that as the molecules become more interesting for practical applications, all methods become more reliable.

We also performed separate statistical analyses on valence and Rydberg states to gain deeper insights into the performance of such methods. We report the resulting indicators in Figure 2 and Tables S8–S10. The $3A_1$ state of pyridazine has been excluded from such statistics due to its mixed character. The TPA response of this state is largely underestimated by several methods, with differences (in au) of 117.0 for Q-CCSD, 117.9 for E-CCSD, 111.1 for I-ADC2, 234.8 for CAM-B3LYP, 302.8 for LC-BLYP33, and 252 for LC-BLYP47. In contrast, CC2 performs very well with an underestimation of only 4.2. For B3LYP, this state is also responsible for the large MAX<sub>s</sub> error of 349.2. For BH&HLYP, this state might be resonance-enhanced (see Table S2), and it has therefore been excluded. Focusing on the Rydberg transitions reveals an improvement in the MAPE for all methods for both $\delta^{\text{TPA}}$ and excitation energies (see Tables S22 and S24), indicating that they are globally easier to describe than valence states. Q-CCSD and E-CCSD emerge as the top performers with very small MAPEs of 10.5% and 12.1%, respectively, closely followed by CAM-B3LYP (25.2%) and BH&HLYP (26.2%). LC-based methods exhibit MAPEs of

32.8% (LC-BLYP33) and 36.0% (LC-BLYP47), while Q-CC2, I-ADC2, and B3LYP consistently show the poorest performances for Rydberg states, with MAPEs in the 40–50% range. Nevertheless all methods demonstrate excellent correlation, with $R^2$ values exceeding 0.95. For valence states we observe a general decline in accuracy compared to Rydberg transitions across all methods. Despite this, Q-CCSD and E-CCSD remain the most satisfying approaches, with MAPEs of 24.8% and 36.6%, respectively. Interestingly, unlike in the case of Rydberg states, Q-CC2 follows closely with a MAPE of 54.2%. As expected, these trends are consistent with those observed for excitation energies (see Table S23), given that CC2 and ADC2 provide more accurate energies for valence transitions than for Rydberg transitions.<sup>85,94</sup> The performances of all TD-DFT methods are comparable, ordered as CAM-B3LYP (55.0%), LC-BLYP47 (57.4%), BH&HLYP (62.4%), LC-BLYP33 (65.0%), and B3LYP (67.7%). I-ADC2 has the worst MAPE at 98.9%. When examining the determination's coefficients, wavefunction methods demonstrate excellent results: 0.98 for both Q-CCSD and E-CCSD, 0.95 for Q-CC2, and 0.90 for I-ADC2. In contrast, TD-DFT values show limited capacity to reproduce trends, with $R^2$ values ranging from 0.66 for B3LYP to 0.49 for BH&HLYP only, despite relatively small absolute and percentage errors<sup>85,94</sup>.

As a further analysis, we investigated the results for valence states with large $\delta^{\text{TPA}}$ (larger than 20 au at the CC3/aug-cc-pVDZ level), since they are more interesting for actual applications. We underline that this subset includes only nine transitions, so these results should be interpreted with caution. The results of MAPE and $R^2$ are reported in Figure 4, while the other indicators can be found in Tables S11–S15. Q-CCSD is again the most adequate method, with a MAPE of 24.4% and

<figure><img src="image_16.png" /><figcaption>**Figure 4.** MAPEs and $R^2$ values of $\delta^{\text{TPA}}$ for all of the wavefunction and TD-DFT methods for the set of valence states showing $20 \le \delta^{\text{TPA}}_{Q-CC3}$. References are obtained at the Q-CC3/aug-cc-pVDZ level. The values of the indicators are reported in Tables S11–S15. The legend is reported in Figures 2 and 3.</figcaption></figure>
an R² of 0.93. Interestingly, for this subset, Q-CC2 and E-CCSD show similar performances: MAPE of 34.3% (44.8%) and R² 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 1B₃g state of tetrazine. In fact its Q-CC3 δTPA of 56.0 au is significantly overestimated by both Q-CCSD (86.6 au) and E-CCSD (118.3 au), while CC2 provides a rather accurate estimate of 47.7 au. Indeed, when we exclude this transition, 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%, but it has a relatively acceptable R² of 0.76. TD-DFT methods perform similarly as in the whole set in terms of MAPE but show extremely poor R² values. This discrepancy can be explained by examining Figure S4: TD-DFT methods tend to underestimate δTPA for small values compared to Q-CC3 for valence states while overestimating δTPA for large response values. This results in an increase in the correlation when the entire set is considered. B3LYP is the “best” XCF, with a MAPE of 44.5% and a low R² of 0.34. The three range-separated functionals show similar results: MAPEs of 47.4%, 47.0%, and 48.3% and R² 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 R² of 0.04. This of course raises questions about actual TD-DFT applications for TPA simulations.

The performances of the wavefunction methods Q-CCSD, E-CCSD, Q-CC2, I-ADC2, and DEA-CCSD (for which raw values are reported in Table S6) 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 accurate reproduction of δTPA by TD-DFT would likely be due to significant error cancellation. Table 6 reports the absolute percentage error (APE) for these methods compared to the reference Q-CC3 calculations. Values with APEs exceeding 200% are not detailed.

<table><caption>**Table 6.** Absolute Percentage Errors (APEs) for the δTPA of Molecules with Double Excitations<sup>a</sup></caption><thead><tr><th>molecule</th><th>state</th><th>Q-CCSD</th><th>E-CCSD</th><th>DEA-CCSD</th><th>Q-CC2</th><th>I-ADC2</th></tr></thead><tbody><tr><td rowspan="2">acrolein</td><td>2A"</td><td>65</td><td>79</td><td>>200</td><td>115</td><td>96</td></tr><tr><td>5A'</td><td>159</td><td>>200</td><td>>200</td><td>38</td><td>43</td></tr><tr><td>butadiene</td><td>1A<sub>g</sub></td><td>118</td><td>173</td><td>23</td><td>>200</td><td>>200</td></tr><tr><td>CycloP</td><td>1A<sub>1</sub></td><td>96</td><td>129</td><td>15</td><td>>200</td><td>>200</td></tr><tr><td>benzene</td><td>1E<sub>2g</sub></td><td>15</td><td>4</td><td>63</td><td>28</td><td>48</td></tr><tr><td>hexatriene</td><td>1A<sub>g</sub></td><td>>200</td><td>>200</td><td>62</td><td>>200</td><td>>200</td></tr></tbody></table>

The DEA-CCSD method provides accurate results for the 1A₁ state of cyclopentadiene and the 1A₈ 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 1A₈ 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 to be quite unreliable; indeed, we have an APE

larger than 200%. For the 2A" 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 5A' state, the second-order methods surprisingly 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 1E₂g state of benzene, with E-CCSD showing the lowest APE (4%) and DEA-CCSD having the highest APE (63%).

In summary, Q-CCSD generally provides acceptable results for states with limited double-excitation character but fails for the 1A₈ state of hexatriene that has a %T₁ of 66.5%. DEA-CCSD improves upon E-CCSD and Q-CCSD for certain systems but can be unreliable for others. Such behavior is already expected, as DEA-CCSD is supposed to be effective for excited states in which mainly excitations from HOMO to LUMO are involved,⁴⁷ while in the considered states of acrolein and benzene the main contributions come from lower occupied and higher unoccupied orbitals.⁸²

## CONCLUSIONS

This study evaluated the performance of several computational methods and basis sets for calculating the vertical two-photon transition strengths (δTPA) of the singlet excited states in small molecules using Q-CC3 results as benchmarks. The analysis underscored the importance of employing a not too compact basis set. Specifically, benchmarks on small molecules using the aug-cc-pVTZ basis set indicated that the aug-cc-pVDZ basis set introduces an error of approximately 10%. Using aug-cc-pVDZ as a reference, we observed that 6-31+G(d) reduces the performances up to 24%, and thus, we advise against the use of such small basis sets for performing reference calculations.

Among the lower-order methods evaluated, Q-CCSD demonstrated superior performance and reliability, showing no significant differences of accuracy between Rydberg and valence transitions. Consequently, Q-CCSD is recommended when feasible. Conversely, I-ADC2 consistently performed the worst across all subsets and cannot be advised from the present work, though its performance for large dyes could possibly be more satisfying. For Rydberg states, E-CCSD emerged as a reliable and cost-effective alternative to Q-CCSD, while Q-CC2 underperformed. Indeed for such cases certain TD-DFT XCFs, such as BH&HLYP and CAM-B3LYP, should be preferred when CCSD calculations are impractical. For valence states, the performance of all methods diminished, particularly for states with large δTPA, which are significant for practical applications. Here E-CCSD showed some limitations, and Q-CC2 appeared to be a valuable choice considering its lower computational scaling. Although TD-DFT methods exhibited excellent percentage errors compared to the reference, their very small correlation suggests that the relative error pattern might significantly change for larger molecules.

Overall, this study highlights the necessity for improved cost-efficient approximations for computing TPA cross sections in large systems. Indeed, the current low-scaling (less than O(N⁶)) state-of-the-art methods appear 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 XCFs. In this regard, the recent development of quadratic response for BSE/GW presents a promising avenue for future exploration.⁶⁶

<footnote><sup>a</sup>Reference values were obtained at the Q-CC3/aug-cc-pVDZ level. CycloP stands for Cyclopentadiene.</footnote>
## ASSOCIATED CONTENT

**Supporting Information**

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acs.jctc.4c01090.

Transitions which have not be included in the statistical analysis, $\delta^{TPA}$ for lower-scaling methods, additional statistical analyses, and the Cartesian coordinates of all molecules (PDF)

## AUTHOR INFORMATION

**Corresponding Authors**

Denis Jacquemin – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; Institut Universitaire de France (IUF), F-75005 Paris, France; orcid.org/0000-0002-4217-0708; Email: Denis.Jacquemin@univ-nantes.fr

Carmelo Naim – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France; Email: carmelo.naim@univ-nantes.fr

**Author**

Robert Zaleśny – Faculty of Chemistry, Wroclaw University of Science and Technology, PL-50370 Wroclaw, Poland; orcid.org/0000-0001-8998-3725

Complete contact information is available at:
https://pubs.acs.org/10.1021/acs.jctc.4c01090

**Notes**

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

This work was supported by the French Agence Nationale de la Recherche (ANR) under Contract ANR-21-CE07-0058-2 (CONDOR). This research used resources of the GLiCID Computing Facility (Ligerien Group for Intensive Distributed Computing, 10.60487/glicid, Pays de la Loire, France). R.Z. acknowledges financial support from National Science Centre (Grant 2018/30/E/ST4/00457).

## REFERENCES

(1) Göppert-Mayer, M. Über Elementarakte Mit Zwei Quantensprüngen. *Ann. Phys.* **1931**, *401*, 273–294.

(2) Kaiser, W.; Garrett, C. G. B. Two-Photon Excitation in $CaF_2$:Eu$^{2+}$. *Phys. Rev. Lett.* **1961**, *7*, 229–231.

(3) Denk, W.; Strickler, J. H.; Webb, W. W. Two-photon Laser Scanning Fluorescence Microscopy. *Science* **1990**, *248*, 73–76.

(4) He, G. S.; Bhawalkar, J. D.; Zhao, C. F.; Prasad, P. N. Optical Limiting Effect in a Two-photon Absorption Dye Doped Solid Matrix. *Appl. Phys. Lett.* **1995**, *67*, 2433–2435.

(5) Belfield, K. D.; Schafer, K. J.; Liu, Y.; Liu, J.; Ren, X.; Stryland, E. V. Multiphoton-Absorbing Organic Materials for Microfabrication, Emerging Optical Applications and Non-destructive Three-dimensional Imaging. *J. Phys. Org. Chem.* **2000**, *13*, 837–849.

(6) Gao, D.; Agayan, R. R.; Xu, H.; Philbert, M. A.; Kopelman, R. Nanoparticles for Two-photon Photodynamic Therapy in Living Cells. *Nano Lett.* **2006**, *6*, 2383–2386.

(7) Dvornikov, A. S.; Walker, E. P.; Rentzepis, P. M. Two-photon Three-dimensional Optical Storage Memory. *J. Phys. Chem. A* **2009**, *113*, 13633–13644.

(8) Wang, B.-G.; König, K.; Halbhuber, K.-J. Two-Photon Microscopy of Deep Intravital Tissues and its Merits in Clinical Research. *J. Microsc.* **2010**, *238*, 1–20.

(9) Drobizhev, M.; Makarov, N. S.; Tillo, S. E.; Hughes, T. E.; Rebane, A. Two-photon Absorption Properties of Fluorescent Proteins. *Nat. Methods* **2011**, *8*, 393–399.

(10) Streets, A. M.; Li, A.; Chen, T.; Huang, Y. Imaging Without Fluorescence: Nonlinear Optical Microscopy for Quantitative Cellular Imaging. *Anal. Chem.* **2014**, *86*, 8506–8513.

(11) Kim, H. M.; Cho, B. R. Small-Molecule Two-photon Probes for Bioimaging Applications. *Chem. Rev.* **2015**, *115*, 5014–5055.

(12) Shen, Y.; Shuhendler, A. J.; Ye, D.; Xu, J.-J.; Chen, H.-Y. Two-photon Excitation Nanoparticles for Photodynamic Therapy. *Chem. Soc. Rev.* **2016**, *45*, 6725–6741.

(13) Olesiak-Banska, J.; Waszkielewicz, M.; Obstarczyk, P.; Samoc, M. Two-Photon Absorption and Photoluminescence of Colloidal Gold Nanoparticles and Nanoclusters. *Chem. Soc. Rev.* **2019**, *48*, 4087–4117.

(14) Pascal, S.; David, S.; Andraud, C.; Maury, O. Near-Infrared Dyes for Two-Photon Absorption in the Short-Wavelength Infrared: Strategies Towards Optical Power Limiting. *Chem. Soc. Rev.* **2021**, *50*, 6613–6658.

(15) He, G. S.; Tan, L.-S.; Zheng, Q.; Prasad, P. N. Multiphoton Absorbing Materials: Molecular Designs, Characterizations, and Applications. *Chem. Rev.* **2008**, *108*, 1245–1330.

(16) Terenziani, F.; Katan, C.; Badaeva, E.; Tretiak, S.; Blanchard-Desce, M. Enhanced Two-photon Absorption of Organic Chromophores: Theoretical and Experimental Assessments. *Adv. Mater.* **2008**, *20*, 4641–4678.

(17) Salem, M. A.; Gedik, M.; Brown, A. Two Photon Absorption in Biological Molecules. In *Handbook of Computational Chemistry; Leszczynski, J., Kaczmarek-Kedziera, A., Puzyn, T., Papadopoulos, M. G., Reis, H., Shukla, M. K., Eds.; Springer: Dordrecht, The Netherlands, 2016; pp 1875–1893. DOI: 10.1007/978-3-319-27282-5_47*.

(18) Meyers, F.; Marder, S.; Pierce, B.; Brédas, J. Electric Field Modulated Nonlinear Optical Properties of Donor-acceptor Polyneses: Sum-over-states Investigation of the Relationship Between Molecular Polarizabilities ($\alpha$, $\beta$, and $\gamma$) and Bond Length Alternation. *J. Am. Chem. Soc.* **1994**, *116*, 10703–10714.

(19) Bartkowiak, W.; Zaleśny, R.; Leszczynski, J. Relation Between Bond-length Alternation and Two-photon Absorption of a Push-Pull Conjugated Molecules: a Quantum-Chemical Study. *Chem. Phys.* **2003**, *287*, 103–112.

(20) Alam, M. M.; Chattopadhyaya, M.; Chakrabarti, S.; Ruud, K. High-polarity Solvents Decreasing the Two-photon Transition Probability of Through-Space Charge-Transfer Systems – A Surprising In Silico Observation. *J. Phys. Chem. Lett.* **2012**, *3*, 961–966.

(21) Wielgus, M.; Zaleśny, R.; Murugan, N. A.; Kongsted, J.; Ågren, H.; Samoc, M.; Bartkowiak, W. Two-photon Solvatochromism II: Experimental and Theoretical Study of Solvent Effects on the Two-Photon Absorption Spectrum of Reichardt’s Dye. *ChemPhysChem* **2013**, *14*, 3731–3739.

(22) Beerepoot, M. T. P.; Alam, M. M.; Bednarska, J.; Bartkowiak, W.; Ruud, K.; Zaleśny, R. Benchmarking the Performance of Exchange-Correlation Functionals for Predicting Two-Photon Absorption Strengths. *J. Chem. Theory Comput.* **2018**, *14*, 3677–3685.

(23) Albota, M.; Beljonne, D.; Brédas, J.-L.; Ehrlich, J. E.; Fu, J.-Y.; Heikal, A. A.; Hess, S. E.; Kogej, T.; Levin, M. D.; Marder, S. R.; McCord-Maughon, D.; Perry, J. W.; Röckel, H.; Rumi, M.; Subramaniam, G.; Webb, W. W.; Wu, X.-L.; Xu, C. Design of Organic Molecules With Large Two-photon Absorption Cross Sections. *Science* **1998**, *281*, 1653–1656.

(24) Zojer, E.; Beljonne, D.; Kogej, T.; Vogel, H.; Marder, S. R.; Perry, J. W.; Brédas, J. L. Tuning the Two-Photon Absorption Response of Quadrupolar Organic Molecules. *J. Chem. Phys.* **2002**, *116*, 3646–3658.

(25) Rossano-Tapia, M.; Olsen, J. M. H.; Brown, A. Two-photon Absorption Cross-sections in Fluorescent Proteins Containing Non-Canonical Chromophores Using Polarizable QM/MM. *Front. Mol. Biosci.* **2020**, *7*, 111.

(26) Knysh, I.; Jassar, M. B.; Ośmiałowski, B.; Zaleśny, R.; Jacquemin, D. In Silico Screening of Two-photon Absorption
Properties of a Large Set of Bis-Difluoroborate Dyes. *ChemPhotoChem* **2022**, *6*, e202200137.

(27) Myers, A. B. Molecular Electronic Spectral Broadening in Liquids and Glasses. *Annu. Rev. Phys. Chem.* **1998**, *49*, 267–295.

(28) Paterson, M. J.; Christiansen, O.; Pawlowski, F.; Jørgensen, P.; Hättig, C.; Helgaker, T.; Salek, P. Benchmarking Two-Photon Absorption with CC3 Quadratic Response Theory, and Comparison with Density-functional Response Theory. *J. Chem. Phys.* **2006**, *124*, 054322.

(29) Macak, P.; Luo, Y.; Norman, P.; Ågren, H. Electronic and Vibronic Contributions to Two-Photon Absorption of Molecules with Multi-branched Structures. *J. Chem. Phys.* **2000**, *113*, 7055–7061.

(30) Painelli, A.; Del Freo, L.; Terenziani, F. Vibronic Contributions to Resonant NLO Responses: Two-Photon Absorption in Push-pull Chromophores. *Chem. Phys. Lett.* **2001**, *346*, 470–478.

(31) Lin, N.; Luo, Y.; Ruud, K.; Zhao, X.; Santoro, F.; Rizzo, A. Differences in Two-Photon and One-Photon Absorption Profiles Induced by Vibronic Coupling: The Case of Dioxaborine Heterocyclic Dye. *ChemPhysChem* **2011**, *12*, 3392–3403.

(32) Kamarchik, E.; Krylov, A. I. Non-Condon Effects in the One- and Two-photon Absorption Spectra of the Green Fluorescent Protein. *J. Phys. Chem. Lett.* **2011**, *2*, 488–492.

(33) Silverstein, D. W.; Jensen, L. Vibronic Coupling Simulations for Linear and Nonlinear Optical Processes: Theory. *J. Chem. Phys.* **2012**, *136*, 064111.

(34) Silverstein, D. W.; Jensen, L. Vibronic Coupling Simulations for Linear and Nonlinear Optical Processes: Simulation Results. *J. Chem. Phys.* **2012**, *136*, 064110.

(35) Hu, H.; Przhonska, O. V.; Terenziani, F.; Painelli, A.; Fishman, D.; Ensley, T. R.; Reichert, M.; Webster, S.; Bricks, J. L.; Kachkovski, A. D.; Hagan, D. J.; Van Stryland, E. W. Two-Photon Absorption Spectra of a Near-Infrared 2-Azaazulene Polymethine Dye: Solvation and Ground-state Symmetry Breaking. *Phys. Chem. Chem. Phys.* **2013**, *15*, 7666–7678.

(36) Rossano-Tapia, M.; Brown, A. Quantum Mechanical/molecular Mechanical Studies of Photophysical Properties of Fluorescent Proteins. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2022**, *12*, e1557.

(37) Petrusevich, E. F.; Reis, H.; Ośmiałowski, B.; Jacquemin, D.; Luis, J. M.; Zaleśny, R. One- and Two-photon Absorption Spectra of Organoboron Complexes: Vibronic and Environmental Effects. *Phys. Chem. Chem. Phys.* **2024**, *26*, 13239–13250.

(38) de Wergifosse, M.; Houk, A. L.; Krylov, A. I.; Elles, C. G. Two-photon Absorption Spectroscopy of Trans-stilbene, Cis-stilbene, and Phenanthrene: Theory and Experiment. *J. Chem. Phys.* **2017**, *146*, 144305.

(39) de Wergifosse, M.; Beaujean, P.; Grimme, S. Ultrafast Evaluation of Two-Photon Absorption with Simplified Time-Dependent Density Functional Theory. *J. Phys. Chem. A* **2022**, *126*, 7534–7547.

(40) Elayan, I. A.; Rib, L.; A. Mendes, R.; Brown, A. Beyond Explored Functionals: A Computational Journey of Two-photon Absorption. *J. Chem. Theory Comput.* **2024**, *20*, 3879–3893.

(41) Beerepoot, M. T. P.; Friese, D. H.; List, N. H.; Kongsted, J.; Ruud, K. Benchmarking Two-Photon Absorption Cross Sections: Performance of CC2 and CAM-B3LYP. *Phys. Chem. Chem. Phys.* **2015**, *17*, 19306–19314.

(42) Bednarska, J.; Zaleśny, R.; Tian, G.; Murugan, N. A.; Ågren, H.; Bartkowiak, W. Nonempirical Simulations of Inhomogeneous Broadening of Electronic Transitions in Solution: Predicting Band Shapes in One- and Two-Photon Absorption Spectra of Chalcones. *Molecules* **2017**, *22*, 1643.

(43) Knippenberg, S.; Rehn, D. R.; Wormit, M.; Starcke, J. H.; Rusakova, I. L.; Trofimov, A. B.; Dreuw, A. calculations of Nonlinear Response Properties Using the Intermediate State Representation and the Algebraic-diagrammatic Construction Polarization Propagator Approach: Two-photon Absorption Spectra. *J. Chem. Phys.* **2012**, *136*, 064107.

(44) Nanda, K. D.; Krylov, A. I. Two-Photon Absorption Cross Sections within Equation-of-Motion Coupled-cluster Formalism

Using Resolution-of-the-identity and Cholesky Decomposition Representations: Theory, Implementation, and Benchmarks. *J. Chem. Phys.* **2015**, *142*, 064118.

(45) Sirimatayanant, S.; Andruniów, T. benchmarking Two-Photon Absorption Strengths of Rhodopsin Chromophore Models with CC3 and CCSD Methodologies: An Assessment of Popular Density Functional Approximations. *J. Chem. Phys.* **2023**, *158*, 094106.

(46) Koch, H.; Christiansen, O.; Jørgensen, P.; Sanchez De Merás, A. M.; Helgaker, T. The CC3Model: An Iterative Coupled Cluster Approach Including Connected Triples. *J. Chem. Phys.* **1997**, *106*, 1808–1818.

(47) Nanda, K. D.; Gulania, S.; Krylov, A. I. theory, Implementation, and Disappointing Results for Two-Photon Absorption Cross Sections within the Doubly Electron-attached Equation-of-Motion Coupled-cluster Framework. *J. Chem. Phys.* **2023**, *158*, 054102.

(48) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansík, B.; Jensen, H. J. A.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. The Dalton Quantum Chemistry Program System. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2014**, *4*, 269–284.

(49) Purvis, G. D.; Bartlett, R. J. A Full Coupled-cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. *J. Chem. Phys.* **1982**, *76*, 1910–1918.

(50) Christiansen, O.; Koch, H.; Jørgensen, P. The Second-order Approximate Coupled Cluster Singles and Doubles Model CC2. *Chem. Phys. Lett.* **1995**, *243*, 409–418.

(51) Whitten, J. L. coulombic Potential Energy Integrals and Approximations. *J. Chem. Phys.* **1973**, *58*, 4496–4501.

(52) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in ab Initio Theory. An Application in MP2 Energy Calculations. *Chem. Phys. Lett.* **1993**, *208*, 359–363.

(53) *Turbomole V7.2* 2017, a Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH 1989–2007, *Turbomole GmbH* since 2007; available from https://www.turbomole.org/.

(54) Balasubramani, S. G.; Chen, G. P.; Coriani, S.; Diedenhofen, M.; Frank, M. S.; Franzke, Y. J.; Furche, F.; Grotjahn, R.; Harding, M. E.; Hättig, C.; Hellweg, A.; Helmich-Paris, B.; Holzer, C.; Huniar, U.; Kaupp, M.; Marefat Khah, A.; Karbalaei Khani, S.; Müller, T.; Mack, F.; Nguyen, B. D.; Parker, S. M.; Perlt, E.; Rappoport, D.; Reiter, K.; Roy, S.; Rückert, M.; Schmitz, G.; Sierka, M.; Tapavicza, E.; Tew, D. P.; van Wüllen, C.; Voora, V. K.; Weigend, F.; Wodyński, A.; Yu, J. M. *Turbomole*: Modular Program Suite for ab Initio Quantum-chemical and Condensed-matter Simulations. *J. Chem. Phys.* **2020**, *152*, 184107.

(55) Friese, D. H.; Hättig, C.; Ruud, K. Calculation of Two-Photon Absorption Strengths with the Approximate Coupled Cluster Singles and Doubles Model CC2 Using the Resolution-of-Identity Approximation. *Phys. Chem. Chem. Phys.* **2012**, *14*, 1175–1184.

(56) Grabarek, D.; Andruniów, T. Assessment of Functionals for TDDFT Calculations of One- and Two-photon Absorption Properties of Neutral and Anionic Fluorescent Proteins Chromophores. *J. Chem. Theory Comput.* **2019**, *15*, 490–508.

(57) Petrusevich, E. F.; Glodek, M.; Antoniak, M. A.; Muziol, T.; Plażuk, D.; Siomra, A.; Nyk, M.; Ośmiałowski, B.; Zaleśny, R.
Difluoroborate-based Bichromophores: Symmetry Relaxation and Two-photon Absorption. *Spectrochim. Acta, Part A* **2023**, 295, 122600.

(58) Bartlett, R. J. The Coupled-cluster Revolution. *Mol. Phys.* **2010**, *108*, 2905–2920.

(59) Dreuw, A.; Wormit, M. The Algebraic Diagrammatic Construction Scheme for the Polarization Propagator for the Calculation of Excited States. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2015**, *5*, 82–95.

(60) Epifanovsky, E.; Gilbert, A. T. B.; Feng, X.; Lee, J.; Mao, Y.; Mardirossian, N.; Pokhilko, P.; White, A. F.; Coons, M. P.; Dempwolff, A. L.; Gan, Z.; Hait, D.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Kussmann, J.; Lange, A. W.; Lao, K. U.; Levine, D. S.; Liu, J.; McKenzie, S. C.; Morrison, A. F.; Nanda, K. D.; Plasser, F.; Rehn, D. R.; Vidal, M. L.; You, Z.-Q.; Zhu, Y.; Alam, B.; Albrecht, B. J.; Aldossary, A.; Alguire, E.; Andersen, J. H.; Athavale, V.; Barton, D.; Begam, K.; Behn, A.; Bellonzi, N.; Bernard, Y. A.; Berquist, E. J.; Burton, H. G. A.; Carreras, A.; Carter-Fenk, K.; Chakraborty, R.; Chien, A. D.; Closser, K. D.; Cofer-Shabica, V.; Dasgupta, S.; de Wergifosse, M.; Deng, J.; Diedenhofen, M.; Do, H.; Ehlert, S.; Fang, P.-T.; Fatehi, S.; Feng, Q.; Friedhoff, T.; Gayvert, J.; Ge, Q.; Gidofalvi, G.; Goldey, M.; Gomes, J.; Gonzalez-Espinoza, C. E.; Gulania, S.; Gunina, A. O.; Hanson-Heine, M. W. D.; Harbach, P. H. P.; Hauser, A.; Herbst, M. F.; Hernandez Vera, M.; Hodecker, M.; Holden, Z. C.; Houck, S.; Huang, X.; Hui, K.; Huynh, B. C.; Ivanov, M.; Jasz, A.; Ji, H.; Jiang, H.; Kaduk, B.; Kahler, S.; Khistyaev, K.; Kim, J.; Kis, G.; Klunzinger, P.; Kozcor-Benda, Z.; Koh, J. H.; Kosenkov, D.; Koulias, L.; Kowalczyk, T.; Krauter, C. M.; Kue, K.; Kunitsa, A.; Kus, T.; Ladjanszki, I.; Landau, A.; Lawler, K. V.; Lefrancois, D.; Lehtola, S.; Li, R. R.; Li, Y.-P.; Liang, J.; Liebenthal, M.; Lin, H.-H.; Lin, Y.-S.; Liu, F.; Liu, K.-Y.; Loipersberger, M.; Luenser, A.; Manjanath, A.; Manohar, P.; Mansoor, E.; Manzer, S. F.; Mao, S.-P.; Marenich, A. V.; Markovich, T.; Mason, S.; Maurer, S. A.; McLaughlin, P. F.; Menger, M. F. S. J.; Mewes, J.-M.; Mewes, S. A.; Morgante, P.; Mullinax, J. W.; Oosterbaan, K. J.; Paran, G.; Paul, A. C.; Paul, S. K.; Pavosevic, F.; Pei, Z.; Prager, S.; Proynov, E. I.; Rak, A.; Ramos-Cordoba, E.; Rana, B.; Rask, A. E.; Rettig, A.; Richard, R. M.; Rob, F.; Rossomme, E.; Scheele, T.; Scheurer, M.; Schneider, M.; Sergueev, N.; Sharada, S. M.; Skomorowski, W.; Small, D. W.; Stein, C. J.; Su, Y.-C.; Sundstrom, E. J.; Tao, Z.; Thirman, J.; Tornai, G. J.; Tsuchimochi, T.; Tubman, N. M.; Veccham, S. P.; Vydrov, O.; Wenzel, J.; Witte, J.; Yamada, A.; Yao, K.; Yeganeh, S.; Yost, S. R.; Zech, A.; Zhang, I. Y.; Zhang, X.; Zhang, Y.; Zuev, D.; Aspuru-Guzik, A.; Bell, A. T.; Besley, N. A.; Bravaya, K. B.; Brooks, B. R.; Casanova, D.; Chai, J.-D.; Coriani, S.; Cramer, C. J.; Cserey, G.; DePrince, A. E.; DiStasio, R. A.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Goddard, W. A.; Hammes-Schiffer, S.; Head-Gordon, T.; Hehre, W. J.; Hsu, C.-P.; Jagau, T.-C.; Jung, Y.; Klamt, A.; Kong, J.; Lambrecht, D. S.; Liang, W.; Mayhall, N. J.; McCurdy, C. W.; Neaton, J. B.; Ochsenfeld, C.; Parkhill, J. A.; Peverati, R.; Rassolov, V. A.; Shao, Y.; Slipchenko, L. V.; Stauch, T.; Steele, R. P.; Subotnik, J. E.; Thom, A. J. W.; Tkatchenko, A.; Truhlar, D. G.; Van Voorhis, T.; Wesolowski, T. A.; Whaley, K. B.; Woodcock, H. L.; Zimmerman, P. M.; Faraji, S.; Gill, P. M. W.; Head-Gordon, M.; Herbert, J. M.; Krylov, A. I. Software for the Frontiers of Quantum Chemistry: An Overview of Developments in the Q-Chem 5 Package. *J. Chem. Phys.* **2021**, *155*, 084801.

(61) Koch, H.; Sánchez de Merás, A.; Pedersen, T. B. Reduced Scaling in Electronic Structure Calculations Using Cholesky Decompositions. *J. Chem. Phys.* **2003**, *118*, 9481–9484.

(62) Aquilante, F.; Gagliardi, L.; Pedersen, T. B.; Lindh, R. Atomic Cholesky Decompositions: A Route to Unbiased Auxiliary Basis Sets for Density Fitting Approximation With Tunable Accuracy and Efficiency. *J. Chem. Phys.* **2009**, *130*, 154107.

(63) Weigend, F.; Kattannek, M.; Ahlrichs, R. approximated Electron Repulsion Integrals: Cholesky Decomposition Versus Resolution of the Identity Methods. *J. Chem. Phys.* **2009**, *130*, 164106.

(64) de Wergifosse, M.; Elles, C. G.; Krylov, A. I. Two-Photon Absorption Spectroscopy of Stilbene and Phenanthrene: Excited-State Analysis and Comparison with Ethylene and Toluene. *J. Chem. Phys.* **2017**, *146*, 174102.

(65) Tucholska, A. M.; Modrzejewski, M.; Moszynski, R. transition Properties From The Hermitian Formulation of The Coupled Cluster Polarization Propagator. *J. Chem. Phys.* **2014**, *141*, 124109.

(66) Rauwolf, N.; Klopper, W.; Holzer, C. Non-linear Light–Matter Interactions from the Bethe–Salpeter Equation. *J. Chem. Phys.* **2024**, *160*, 061101.

(67) Salek, P.; Vahtras, O.; Guo, J.; Luo, Y.; Helgaker, T.; Ågren, H. Calculations of Two-photon Absorption Cross Sections by Means of Density-functional Theory. *Chem. Phys. Lett.* **2003**, *374*, 446–452.

(68) Parker, S. M.; Rappoport, D.; Furche, F. Quadratic Response Properties From TDDFT: Trials and Tribulations. *J. Chem. Theory Comput.* **2018**, *14*, 807–819.

(69) Barca, G. M. J.; Bertoni, C.; Carrington, L.; Datta, D.; De Silva, N.; Deustua, J. E.; Fedorov, D. G.; Gour, J. R.; Gunina, A. O.; Guidez, E.; Harville, T.; Irle, S.; Ivanic, J.; Kowalski, K.; Leang, S. S.; Li, H.; Li, W.; Lutz, J. J.; Magoulas, I.; Mato, J.; Mironov, V.; Nakata, H.; Pham, B. Q.; Piecuch, P.; Poole, D.; Pruitt, S. R.; Rendell, A. P.; Roskop, L. B.; Ruedenberg, K.; Sattasathuchana, T.; Schmidt, M. W.; Shen, J.; Slipchenko, L.; Sosonkina, M.; Sundriyal, V.; Tiwari, A.; Galvez Vallejo, J. L.; Westheimer, B.; Wloch, M.; Xu, P.; Zahariev, F.; Gordon, M. S. Recent Developments in the General Atomic and Molecular Electronic Structure System. *J. Chem. Phys.* **2020**, *152*, 154102.

(70) Rinkevicius, Z.; Li, X.; Vahtras, O.; Ahmadzadeh, K.; Brand, M.; Ringholm, M.; List, N. H.; Scheurer, M.; Scott, M.; Dreuw, A.; Norman, P. Veloxchem: A Python-driven Density-functional Theory Program for Spectroscopy Simulations in High-performance Computing Environments. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2020**, *10*, e1457.

(71) Beerepoot, M. T. P.; Friese, D. H.; Ruud, K. Intermolecular Charge Transfer Enhances Two-Photon Absorption in Yellow Fluorescent Protein. *Phys. Chem. Chem. Phys.* **2014**, *16*, 5958–5964.

(72) Salem, M. A.; Brown, A. Two-Photon Absorption in Fluorescent Protein Chromophores: TDDFT and CC2 Results. *J. Chem. Theory Comput.* **2014**, *10*, 3260–3269.

(73) Choluj, M.; Alam, M. M.; Beerepoot, M. T. P.; Sitkiewicz, S. P.; Matito, E.; Ruud, K.; Zaleśny, R. Choosing Bad versus Worse: Predictions of Two-Photon-Absorption Strengths Based on Popular Density Functional Approximations. *J. Chem. Theory Comput.* **2022**, *18*, 1046–1060.

(74) Ahmadzadeh, K.; Li, X.; Rinkevicius, Z.; Norman, P.; Zaleśny, R. Toward Accurate Two-photon Absorption Spectrum Simulations: Exploring the Landscape Beyond the Generalized Gradient Approximation. *J. Phys. Chem. Lett.* **2024**, *15*, 969–974.

(75) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). *Chem. Phys. Lett.* **2004**, *393*, 51–57.

(76) Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. MN15: A Kohn–Sham Global–Hybrid Exchange–Correlation Density Functional with Broad Accuracy for Multi-reference and Single-reference Systems and Noncovalent Interactions. *Chem. Sci.* **2016**, *7*, S032–S051.

(77) Grotjahn, R.; Furche, F. Gauge-Invariant Excited-State Linear and Quadratic Response Properties within the meta-Generalized Gradient Approximation. *J. Chem. Theory Comput.* **2023**, *19*, 4897–4911.

(78) Grotjahn, R.; Furche, F. Comment On: "Toward Accurate Two-Photon Absorption Spectrum Simulations: Exploring the Landscape Beyond the Generalized Gradient Approximation. *J. Phys. Chem. Lett.* **2024**, *15*, 6237–6240.

(79) Rossano-Tapia, M.; Brown, A. Determination of Two-Photon-Absorption Cross Sections using Time-Dependent Density Functional Theory Tight Binding: Application to Fluorescent Protein Chromophores. *J. Chem. Theory Comput.* **2019**, *15*, 3153–3161.

(80) de Wergifosse, M.; Grimme, S. the Exact Integral Simplified Time-Dependent Density Functional Theory (XSTD-DFT). *J. Chem. Phys.* **2024**, *160*, 204110.

(81) Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. A Mountaineering Strategy to Excited States:
Highly Accurate Reference Energies and Benchmarks. *J. Chem. Theory Comput.* **2018**, *14*, 4360–4379.

(82) Loos, P.-F.; Boggio-Pasqua, M.; Scemama, A.; Caffarel, M.; Jacquemin, D. Reference Energies for Double Excitations. *J. Chem. Theory Comput.* **2019**, *15*, 1939–1956.

(83) Loos, P.-F.; Lipparini, F.; Boggio-Pasqua, M.; Scemama, A.; Jacquemin, D. A Mountaineering Strategy to Excited States: Highly Accurate Energies and Benchmarks for Medium Sized Molecules. *J. Chem. Theory Comput.* **2020**, *16*, 1711–1741.

(84) Loos, P.-F.; Scemama, A.; Jacquemin, D. The Quest for Highly Accurate Excitation Energies: A Computational Perspective. *J. Phys. Chem. Lett.* **2020**, *11*, 2374–2383.

(85) Véril, M.; Scemama, A.; Caffarel, M.; Lipparini, F.; Boggio-Pasqua, M.; Jacquemin, D.; Loos, P.-F. QUESTDB: A Database of Highly Accurate Excitation Energies for the Electronic Structure Community. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2021**, *11*, e1517.

(86) Weigend, F.; Köhn, A.; Hättig, C. efficient Use of the Correlation Consistent Basis Sets in Resolution of the Identity MP2 Calculations. *J. Chem. Phys.* **2002**, *116*, 3175–3183.

(87) Gulania, S.; Kjønstad, E. F.; Stanton, J. F.; Koch, H.; Krylov, A. I. Equation-of-Motion Coupled-Cluster Method with Double Electron-Attaching Operators: Theory, Implementation, and Benchmarks. *J. Chem. Phys.* **2021**, *154*, 114115.

(88) Kristensen, K.; Kauczor, J.; Thorvaldsen, A. J.; Jørgensen, P.; Kjærgaard, T.; Rizzo, A. Damped Response Theory Description of Two-Photon Absorption. *J. Chem. Phys.* **2011**, *134*, 214104.

(89) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. *J. Phys. Chem.* **1994**, *98*, 11623–11627.

(90) Becke, A. D. a New Mixing of Hartree-Fock and Local Density-Functional Theories. *J. Chem. Phys.* **1993**, *98*, 1372–1377.

(91) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A Long-range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals. *J. Chem. Phys.* **2001**, *115*, 3540–3544.

(92) Kossoski, F.; Boggio-Pasqua, M.; Loos, P.-F.; Jacquemin, D. Reference Energies for Double Excitations: Improvement and Extension. *J. Chem. Theory Comput.* **2024**, *20*, 5655–5678.

(93) Wielgus, M.; Bartkowiak, W.; Samoc, M. Solvent Effects on Two-photon Absorption Cross Section of 4-Dimethylamino-4'-Nitrostilbene (DANS). *Chem. Phys. Lett.* **2012**, *554*, 113–116.

(94) Tajti, A.; Szalay, P. G. Investigation of the Impact of Different Terms in the Second Order Hamiltonian on Excitation Energies of Valence and Rydberg States. *J. Chem. Theory Comput.* **2016**, *12*, 5477–5482.