# Two-photon Absorption Strengths of Small- and Medium-Sized 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- and Medium-Sized Molecules: Reference CC3 Values and Benchmarks. 2024. [https://hal.science/hal-04839952](https://hal.science/hal-04839952)

HAL Id: hal-04839952

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

Preprint submitted on 17 Dec 2024

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- and
Medium-Sized Molecules: Reference CC3 Values and
Benchmarks

Carmelo Naim,*,† Robert Zaleśny,‡ and Denis Jacquemin*,†,¶

†Nantes Université, CNRS, CEISAM UMR 6230, F-44000, Nantes, France.

‡Faculty of Chemistry, Wrocław University of Science and Technology, Wyb. Wyspiańskiego 27, PL-50370 Wrocław, Poland.

¶Institut Universitaire de France (IUF), F-75005, Paris, France.

E-mail: carmelo.naim@univ-nantes.fr; denis.jacquemin@univ-nantes.fr

Abstract

We present the largest dataset of highly-accurate vertical and degenerate two-photon transition strengths ($\delta^{TPA}$) for standard small- and medium-sized organic molecules, calculated using the quadratic response implementation of the third-order coupled cluster method that includes iterative triples (Q-CC3). The aug-cc-pVTZ basis set was used for all small molecules, while medium-sized molecules were assessed with aug-cc-pVDZ and 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 both small basis sets and, alternative wavefunction methods when Q-CC3 calculations become beyond reach. These methods include quadratic (Q) response and equation of motion CCSD approximations, Q-CC2, second-order algebraic diagrammatic construction in its intermediate state representation (I-ADC2), as well as time-dependent density functional theory (TD-DFT) with a set of 5 commonly used exchange-correlation functionals. This extensive analysis provides a quantitative assessment of these methods, revealing how different system sizes, response intensities, and types of transitions affect their performances.

TOC Graphic

<figure><img src="image_2.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 collaborators.² The TPA process can be quantified by the probability of absorption of two photons, which is usually reported as a cross-section ($\sigma^{\text{TPA}}$). In contemporary applications, materials or dyes exhibiting substantial $\sigma^{\text{TPA}}$ are extensively used in various technological fields, including optical data storage, 3D microfabrication, optical power limit, bioimaging, and photodynamic therapy.³⁻¹⁴ For example, bioimaging applications require the product of two-photon absorption cross-section and fluorescence quantum yield larger than 50 GM (1 GM = cm⁴·s·photon⁻¹).¹¹ In the case of organic dyes $\sigma^{\text{TPA}}$ values can reach several thousands of 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/properties relationships¹⁶,¹⁸⁻²² and furthermore, they can be routinely employed to design new molecules with significant $\sigma^{\text{TPA}}$, obviating the need for generally expensive and complex experiments.¹⁶,²³⁻²⁶ The $\sigma^{\text{TPA}}$ for the absorption of two photons having the same frequency can be expressed as:

$$ \sigma^{\text{TPA}} = \frac{4\pi^3 \alpha a_0^5 \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^{\text{TPA}}$ is the orientationally-averaged two-photon transition strength, which can be determined based on electronic-structure theories (see below). Despite significant developments, an accurate theoretical estimation of $\sigma^{\text{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 in a cost-effective manner. Validation of the results of new computational methodologies can be approached in two ways: by comparisons with experimental data, or by comparison with benchmark values. The former approach presents several challenges, as reproducing experimental $\sigma^{\text{TPA}}$ requires consideration of multiple factors such as solvent effects and vibrational couplings.⁹,²¹,²⁵,²⁹⁻³⁷ Moreover, $\sigma^{\text{TPA}}$ depends on $\omega^2$ and $\delta^{\text{TPA}}$, each introducing distinct sources of error in the computations, potentially leading to misleading interpretations.³⁸ The choice of an appropriate broadening function and the proper selection of constants in Eq. 1 are also subjects of debate.³⁹,⁴⁰ For these reasons, direct comparison of $\delta^{\text{TPA}}$ with more accurate (albeit expensive) wavefunction methods is generally preferred to establish the accuracy of computational methods.²⁸,⁴¹⁻⁴³ In the case of a linearly polarized light, for an excited state $i$, the general formulation of $\delta^{\text{TPA}}$

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}^{\nu \mu} + M_{i \leftarrow 0}^{\nu \mu} M_{0 \leftarrow i}^{\mu \nu}] \quad (2) $$

with $\mu, \nu \in x, y, z$ and $M_{i \leftarrow 0}^{\mu \nu}$ and $M_{0 \leftarrow i}^{\mu \nu}$ representing the left and right second-order transition moments, which are equal in Hermitian theories, see Eq. S1 in the Supporting Information (SI) for their expressions.²⁸ For most organic molecules, the coupled cluster (CC) methods can be viewed as the most reliable approaches for computing $\delta^{\text{TPA}}$. In addition, the hierarchy of CC methods provides a systematic and well defined path for improving the accuracy by incorporating higher excitation orders in the wavefunction expansion. Among these methods, the approximate coupled cluster singles, doubles, and triples model (CC3),⁴⁴ in its quadratic response implementation (Q-CC3),²⁸ stands out as the most reliable for TPA calculations when a monodeterminantal description of the electronic state is appropriate. However, its computational cost, scaling as $O(N^7)$ (where $N$ is the number of orbitals), limits its application to small systems, preventing the study of large molecular assemblies. This is a significant limit since the molecules of interest for actual TPA applications are actually quite large $\pi$-conjugated molecules. Notably, just a few Q-CC3 reference TPA strengths are available in the literature. Paterson and collaborators studied a series of tiny molecules namely, water, ethylene, formaldehyde, diacetylene,²⁸ whereas Nanda and collaborators computed the TPA of twisted ethene, as well as the lowest excited $A_g$ states of butadiene, and hexatriene.⁴⁵ Moreover, Siriatayant and Andruniów determined Q-CC3 TPA strengths for the first excited state of a model chromophore for retinal.⁴³ All of these reference data have been computed using the Dalton code.⁴⁶

When Q-CC3 becomes beyond computational reach coupled cluster singles and doubles method (Q-CCSD)⁴⁷ is considered an excellent alternative, as Paterson and collaborators showed that Q-CCSD, with a scaling of $O(N^6)$ gives a maximum of 20% error on $\delta^{\text{TPA}}$ computed at Q-CC3 level.

Alternatively Q-CC2⁴⁸ was exhaustively employed for TPA calculations, though it has been reported to exhibit a general overestimation trend (by factor between 1.4 and 2) as compared to Q-CC3,²⁸,³⁹,⁴² which might however become smaller when considering medium-sized molecular systems.⁴³ Its computational scaling is $O(N^5)$, yet the resolution of identity (RI) scheme⁴⁹,⁵⁰ implemented for instance in Turbomole⁵¹,⁵² helps improving its computational efficiency, making Q-CC2 applicable to quite large systems.⁵³⁻⁵⁵

An alternative to quadratic response CC methods is provided by the expectation-value approach,⁵⁶ widely used for property computations, with the equation-of-motion CC theories (E-CC)⁵⁶ or the algebraic diagrammatic construction in the intermediate state representation (I-ADC)⁴¹,⁵⁷ schemes. An efficient E-CCSD method⁴² is implemented in Q-Chem⁵⁸ employing both the RI and Cholesky decomposition (CD)⁵⁹⁻⁶¹ schemes, which allows significant time saving and reduced memory requirements compared to standard Q-CCSD, with a decrease of the performances of about 5% only on small molecules.⁴² This method has been successfully
applied as a reference for studying medium-sized systems, e.g., stilbene, phenanthrene,⁶² as well as fluorescent protein chromophores.³⁹,⁴² Recently, Nanda and Krylov, developed the so-called Double Electron Attached CCSD (DEA-CCSD) method, with the aim of improving the accuracy of E-CCSD for states having a double excitation character, such as the A₉ states of ethene, butadiene, and hexatriene.⁴⁵ Within the ADC family the I-ADC2 and I-ADC(3/2) methods, respectively scaling *O*(*N*⁵) and *O*(*N*⁷), have been implemented in Q-Chem by Knippengerg and collaborators.⁴¹ Using Q-CC3 references for hydrogen fluoride, water, and ethylene, they demonstrate that I-ADC2 performs slightly better than Q-CC2.⁴¹ I-ADC(3/2) outperforms both second-order methods, but it results in an average percentage error four times larger than its Q-CCSD counterpart.⁴¹ It is also worth mentioning that methods employing the polarization propagator approach, deriving from a Hermitian theory, prevent the occurrence of nonphysical behavior in the calculation of properties. Although Hermitian coupled-cluster (CC) formulations have been proposed,⁶³ as far as the authors knowledge no implementation of quadratic response theory allowing for computing two-photon transition strengths is currently available. As a further alternative method, let us underline the recent implementation of the computation of nonlinear optical properties with the Bethe-Salpeter equation approach by Rauwolf and co-workers.⁶⁴

The quadratic response formalism can also be applied to compute TPA strengths using Time-Dependent Density Functional Theory (TD-DFT),⁶⁵,⁶⁶ such approach is currently employed in various software packages, e.g., Dalton,⁴⁶ GAMESS,⁶⁷ Turbomole,⁵¹,⁵² and VeloxChem.⁶⁸ TD-DFT can be routinely employed for absolute TPA calculations in medium and large dyes.²²,²⁶,⁵⁴,⁶⁹⁻⁷² As always, the accuracy of TD-DFT estimates strongly depends on the selected exchange correlation functional (XCF). It has indeed been observed how the choice of a proper XCF is particularly challenging for computing the TPA, as standard XCF considerably underperforms with respect to CC methods.⁷¹ CAM-B3LYP⁷³ is often considered suitable for δTPA, as it reproduces the Q-CC2 trends,⁵⁴,⁷¹ although it generally underestimates Q-CC2's δTPA by a factor of 3.³⁸,³⁹,⁶⁹,⁷⁰ Recent studies suggest that novel highly parameterized XCF, such as MN15,⁷⁴ improve CAM-B3LYP results for push-pull systems.³⁸,⁷²,⁷⁵,⁷⁶ Besides quadratic response has been recently implemented for large scale DFT methods,⁷⁷⁻⁷⁹ enabling a qualitative estimation of TPA strengths on extensive systems (up to thousands of atoms).

In this work, with the aim of including TPA data to the QUEST database⁸⁰⁻⁸³ of highly accurate reference values for excited state energies and properties we obtained a collection of δTPA computed at Q-CC3 level for a series of small- and medium-size benchmark molecules. With these data at hand we benchmark several wavefunction methods which could be used when Q-CC3 becomes unreachable: Q-CCSD, E-CCSD, Q-CC2 and I-ADC2, as well as TD-DFT with XCF which are frequently used in the study of medium and large systems.

## 2. COMPUTATIONAL DETAILS

This study includes 30 organic molecules containing from 1 to 6 non-hydrogen atoms, see Scheme 1. The correspond-

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

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

### 4 (non-H) atom molecules

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

### 5 (non-H) atom molecules

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

### 6 (non-H) atom molecules

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

Scheme 1: Schematic representation of the molecules considered in this study.

ing 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₁. 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. For example, the first excited state of acrolein is denoted as 1A'.

The importance of double excitations in each state was assessed using the %T₁ value, which provides the percentage of contribution of all single excitations in an excited state. States with %T₁ smaller than 80% are considered to have a significant double excitation character. We reported the %T₁ values computed at CC3/aug-cc-pVTZ level for transition belonging to small-molecules and at CC3/aug-cc-pVDZ for medium size ones.

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 a strongly mixed character, rendering their assignment challenging. Given that experimental measurements typically cannot distinguish states with tiny 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 δTPA reported should be multiplied by two. A basis set analysis is presented below for Q-CC3 calculations, comparing the δTPA obtained with the 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.

Aug-cc-pVDZ was employed for the benchmarks, coupled with the corresponding auxiliary basis sets when necessary.⁸⁴ States with δTPA Q-CC3 < 3 a.u. are considered irrelevant for the purposes of the paper and they have not been 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 such reason acetylene and streptocyanine-1 have been excluded from the statistical analysis, indeed the 1A<sub>u</sub> and 2A<sub>u</sub> states of the former are TPA forbidden, while the 1B<sub>2</sub> state of the latter (the only state studied) have δTPA value below assumed threshold. All CC calculations use the frozen-core approximation, with the length gauge for the interaction operator. Standard convergence criteria (10<sup>-7</sup> a.u.) was 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 parameters: 10<sup>-9</sup> a.u. for ground-state energies, 10<sup>-8</sup> a.u. for CC equations, and 10<sup>-7</sup> a.u. for EOM Davidson convergence. The CD of two-electron integrals uses the standard threshold, 10<sup>-3</sup> a.u., as it was observed to be sufficient for the TPA strengths.⁴⁵,⁸⁵ DEA-CCSD calculations have been performed on systems with a significant double excitation character on Q-Chem using the same thresholds as its EOM-CCSD version (indeed the RI-CD approximations are not yet implemented for this method), and employing the charge +2 for reproducing the double excitations.⁴⁵,⁸⁵ I-ADC2 calculations are achieved with Q-Chem using default convergence criteria. The coupled cluster implementations employed in Dalton and Q-chem is non-Hermitian, therefore, 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 left and right moments in Eq. 2 are negative and represent at least 20% of the total value of δTPA are excluded from the analysis. The excluded cases are documented in Table S1 in the SI. 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 are also excluded and detailed in Table S2 in the SI. 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 have been 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

respectively 20% and 50% of exact exchange and three range-separated functionals, namely CAM-B3LYP⁷³ which has an exact exchange going from 19 to 65% and a range separation parameter ω of 33 Bohr⁻¹, and two different versions of the LC-BLYP⁸⁹ for which the exact exchange ranges from 0 to 100% and the range separation parameter of 0.47 and 0.33 Bohr⁻¹.

The statistical indicators employed below are the mean average percentage error (MAPE), the mean signed error (MSE), mean absolute error (MAE), root mean squared deviation (RMSD), maximum positive and negative error (MAX+, and MAX−), which expressions are reported in Eqs.3-8:

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

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

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

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

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

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

in which x<sub>i</sub> is the value obtained with the low-cost method, t<sub>i</sub> is the reference value, and n is the number of values. The indicator of correlation, R<sup>2</sup>, 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 a significant double excitation character.

### 3.1 Reference values and basis set analysis

Basis sets effects on Q-CC3 calculations have been firstly analyzed on compounds containing 1-to-5 non hydrogen atoms. Excitation energies, and δTPA values for the excited states in which δTPA > 3 [a.u.] at the Q-CC3/aug-cc-pVTZ level are listed in Table 1.

The 50 states considered in this analysis encompass 18 valence states (π → π*, 16 transitions and n → π*, 2 transitions), among them six exhibit a significant double excitation character, each with different contribution of the double excitations in their composition: 5A<sub>g</sub> of ethylene (20.0%), 2A" (79.4%) and 5A' (75.0%) of acrolein, 1A<sub>g</sub> of butadiene (75.1%), and 1A<sub>1</sub> of cyclopentadiene (78.9%); the values in parenthesis giving the %T<sub>1</sub>. The 5A<sub>g</sub> state of ethylene, characterized by a high degree of double excitation, is an-
<table><caption>Table 1: 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 Pople's 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">δTPA [a.u.]</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="4">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>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'</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"</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'</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'</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>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 rowspan="2">Methylenecycloprop.</td><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>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 rowspan="4">Cyclopentadiene</td><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>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 rowspan="4">Furan</td><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>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 rowspan="4">Pyrrole</td><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>1</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>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 rowspan="4">Thiophene</td><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>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 rowspan="4">Streptocyanine-3</td><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+</td><td>110.9</td><td>115.8</td><td>68.3+</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+</td><td>11.5</td><td>9.8</td><td>44.7+</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></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></td><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>

<footnote>+ Mixed states (summed in the stat. analysis)</footnote>
ticipated to be largely affected not only by the basis set, as previously demonstrated for energies by one of us,⁸¹,⁹⁰ but more importantly by the level of theory. Indeed CC3 is far from an ideal method for such genuine double excitation. The $\delta^{TPA}$ values at that state are 894.1, 1636.7 and 1947.9 a.u. with aug-cc-pVTZ, aug-cc-pVDZ, and 6-31+G(d), respectively. These substantial differences lead us to conclude that basis sets smaller than aug-cc-pVTZ 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; Q-CC3 being likely well suited.

The remaining 32 states have a Rydberg character. Table 1 illustrates that for the thiophene molecule, the $1A_2$ and $2A_2$ states look strongly mixed with the 6-31+G(d) whereas the situation is clearer with larger basis sets. Indeed the values are: 110.9 a.u. ($1A_2$), and 11.5 a.u. ($2A_2$) with aug-cc-pVTZ, while 115.8 a.u. ($1A_2$), and 9.8 a.u. ($2A_2$) with aug-cc-pVDZ. The Q-CC3/6-31+G(d) method is enabled to distinguish these two states, separated by 0.06 eV only, resulting in a splitting of the $\delta^{TPA}$ values (68.3 and 44.7 a.u.). 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 2, using the aug-cc-pVTZ values as reference. According to the reported statistical indicators, the

<table><caption>**Table 2: Statistical analysis for $\delta^{TPA}$ values of small-molecules (see Table 1 and for details). Reference values are 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 [a.u.]</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 [a.u.]</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 [a.u.]</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> [a.u.]</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> [a.u.]</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 [a.u.]</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>

aug-cc-pVDZ results exhibit a MAPE of 11.5% and a MAE of 6.7 a.u. indicating a 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 both in the MAPE (22.5%) and MAE (15.1 a.u.). It is noteworthy that all error indicators of aug-cc-pVDZ results are approximately twice smaller than those of 6-31+G(d). Similar conclusions can be drawn when analyzing valence and Rydberg states separately, as also illustrated in Table 2.

For states with 6 non-hydrogens atoms, Q-CC3/aug-cc-pVTZ becomes too computational demanding. Encouraged by the performances of aug-cc-pVDZ for smaller compounds we took aug-cc-pVDZ data as references to evaluate the performances of 6-31+G(d) when larger compounds are considered as well. The $\delta^{TPA}$ values as well as energies for 34 excited states showing a $\delta^{TPA} > 3$ a.u. with aug-cc-pVDZ

are reported in Table 3, while states with smaller response values are reported in Tables S3 and S4 in the SI. This set encompasses 13 valence states ($7 \pi \rightarrow \pi^*$ and $6 n \rightarrow \pi^*$) transitions. The $1A_g$ state of hexatriene and the $1E_{2g}$ state of benzene exhibit a substantial double character with $\%T_1$ of 66.5% and 67.8%, respectively. Within this set there are also 20 Rydberg states, and one mixed state, the $3A_1$ state of pyridazine. We evaluated the performance of 6-31+G(d) with respect to aug-cc-pVDZ on 82 states: 48 states collected in Table 1 and 34 states on Table 3. The results are provided in Table 4. For the full set we found that $\delta^{TPA}$ with the 6-31+G(d) basis significantly differ from the aug-cc-pVDZ values with a MAPE of 24.7%, and a MAE of 15.1 a.u. These values are essentially unchanged with respect to the ones of Table 2, indicating that basis set differences are not strongly dependent on the considered molecules. The large $MAX_+$ in Table 4 comes from the $1A_g$ state of hexatriene, with an overestimation that might be ascribed to the worse description of double excitations with the compact basis. By excluding the 6 states where double excitations are not negligible (i.e., $\%T_1 < 80 \%$), we observe a slight increase of the MAPE but a large improvement of the others indicators. Table 4 also indicates a generally better performance of 6-31+G(d) for the valence states, although the indicators stay in the same range for Rydberg transitions. Further comment should be made on the size of $\delta^{TPA}$ and its implications for practical applications. The studied molecules exhibit rather small values of $\delta^{TPA}$ and corresponding $\sigma^{TPA}$. For instance, one of the largest $\delta^{TPA}$ values was found for the $1A_g$ transition in hexatriene (1358.3 a.u.), which translates to a $\sigma^{TPA}$ 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.⁹¹ 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 Fig. 1, we provide the error distribution patterns for all methods included in this study. The values of MAE, MSE, and SDE are given in a.u. as all values in the following.

Obviously, CCSD methods significantly outperform the other methods in all 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, a MSE of 0.1, and a 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, a MSE of 3.9, and a SDE of 17.7; these differences can be mainly ascribed to the valence states (see below). Conversely, Q-CC2 and I-ADC2 methods exhibit similar but significantly inferior performances than CCSD, with MAEs of 21.5 (22.8) for Q-CC2 (I-ADC2), over four times larger than the Q-CCSD value. Moreover, they demonstrate a quite large overestimation of $\delta^{TPA}$, resulting in a MSE of 20.2 (18.8) and a large error spread with SDE of 49.7 (34.4) for Q-CC2 (I-ADC2). Among the tested TD-DFT methods, we observe, as expected, different results based on the type of functional and the amount of exact exchange included. The most accurate methods are
<table><caption>Table 3: Q-CC3 excitation energies (ΔE) and two-photon transition strengths (δTPA) of molecules with 6 non hydrogen atoms computed with the aug-cc-pVDZ (a-DZ), and Pople's 6-31+G(d) (p-DZ) basis sets. Values of %T₁ at the Q-CC3/a-DZ level are also reported. Molecules labelled with * have a non-Abelian group symmetry, their δTPA must be multiplied by two when compared with experiments (see section 2 for more details).</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> [a.u.]</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="4">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 rowspan="5">Pyridazine</td><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>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 rowspan="4">Pyridine</td><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>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 rowspan="5">Pyrimidine</td><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="2">Tetrazine</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 rowspan="4">Triazine*</td><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>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 rowspan="2">Triazine*</td><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>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 rowspan="3">Triazine*</td><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>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>
<figure><img src="image_7.png" /><figcaption>Figure 1: Error distribution, MAE, MSE, and SDE in a.u. with respect to reference Q-CC3/aug-cc-pVDZ $\delta^{TPA}$ for states without a large contribution of double excitations.</figcaption></figure>

<table><caption>**Table 4:** Statistical analysis of $\delta^{TPA}$ computed at 6-31+G(d) for the states in Table 1 and 3. Reference values are computed at 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 [a.u.]</td><td>-9.2</td><td>-5.1</td><td>-8.6</td><td>3.8</td></tr><tr><td>MAE [a.u.]</td><td>15.1</td><td>11.3</td><td>13.4</td><td>5.6</td></tr><tr><td>RMSD [a.u.]</td><td>38.5</td><td>18.1</td><td>20.4</td><td>9.0</td></tr><tr><td>MAX<sub>+</sub> [a.u.]</td><td>306.7</td><td>62.0</td><td>62.0</td><td>5.6</td></tr><tr><td>MAX<sub>-</sub> [a.u.]</td><td>68.7</td><td>68.7</td><td>68.7</td><td>31.2</td></tr><tr><td>SDE [a.u.]</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>

BH&HLYP and CAM-B3LYP functionals, yielding 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 that show similar trends. LC-BLYP33 exhibits a MAE of 19.1, a 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 a MSE of -14.5, more than twice as large as its $\omega = 33$ Bohr$^{-1}$ counterpart. B3LYP is the worst-performing method among those tested for all indicators, with a MAE of 25.5, a MSE of 16.9, and a SDE of 54.1. Interestingly, all TD-DFT methods (excluding B3LYP) outperform or slightly underperforms the Q-CC2 and I-ADC2 wavefunction methods for all indicators for a fraction of the computational cost. This also indicates that benchmarking TD-DFT on the basis of Q-CC2 results might

not always be a fully trustworthy procedure.

In Fig. 2, the MAPE and $R^2$ indicators are reported, including a separate analysis for valence and Rydberg transitions that we discuss below.

Looking at the MAPE for all 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 MAPE 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 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 if the size of the systems affects the accuracy of the methods, we performed separate statistical analyses on systems with 1-to-4, and 5-6 non-hydrogen atoms. We collected the related indicators in Tables S16–S18 in the SI. From this analysis, we observe small variations in the trends for each method when moving from smaller to larger molecules. Interestingly, we observe how Q-CCSD, E-CCSD, Q-CC2, and I-ADC2 methods reduce their performances. For example the MAPE of Q-CCSD passes from 6.5% to 18.4%, of Q-CC2 from 40.4% to 45.1%, and of I-ADC2 from 46.1% to 60.8% when going from small to large molecules. In contrast, the performances of TD-DFT methods tend to improve, with
<figure><img src="image_8.png" /><figcaption>Figure 2: MAPE and $R^2$ of $\delta^{\text{TPA}}$ for all the wavefunction and TD-DFT methods for the sets with all singly-excited, Rydberg, and valence states. References are the $\delta^{\text{TPA}}_{\text{Q-CC}_3}$/aug-cc-pVDZ calculations with a threshold of 3 [a.u.]. The values of the indicators are reported in Tables S8–S10 in the SI.</figcaption></figure>

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 a MAPE of 36.8% for small molecules and 45.6% for the large ones.

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}}_{\text{Q-CC}_3} < 20$, $t_2 = 20 \le \delta^{\text{TPA}}_{\text{Q-CC}_3} < 80$, and $t_3 = 80 \le \delta^{\text{TPA}}_{\text{Q-CC}_3}$. In Fig. 3 we report the MAPE, and $R^2$ computed for these three groups, while the rest of the indicators are reported in Tables S11–S13 in the SI.

Looking at both MAPE and $R^2$, we can observe how the general trends of Fig. 2 are replicated for all the methods at different thresholds. However, it is noteworthy that the MAPE ($R^2$) decreases (increases) when increasing the $\delta^{\text{TPA}}_{\text{Q-CC}_3}$ threshold. Indeed, looking at 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 to 5.4% (8.1%) for the $t_3$ interval. Similar decreases are obtained for the other methods (see Tables S11–S13 in the SI). Analogous trends can be found for the correlation. 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 excell with an $R^2$ of 0.85). On the other hand in the $t_3$ interval all methods show large correlation ($R^2 > 0.91$). It is therefore a good news that as the molecules get more interesting for practical applications all methods become more trustworthy.

We also performed separate statistical analyses on valence and Rydberg states to gain deeper insights into the performances of such methods. We report the resulting indicators in Fig. 2 and in Tables S8–S10 in the SI. 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 a.u.) 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. Whereas, CC2 performs very well with an underestimation of only 4.2, within B3LYP, this state is also responsible for the large $\text{MAX}_+$ error of 349.2. For BH&HLYP, this state might be resonance enhanced (see Table S2 in the SI), and it has therefore been excluded. Focusing on the Rydberg

transitions reveals an improvement in the MAPE for all methods, 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 with MAPEs of 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 as compared to Rydberg transitions across all methods. Despite this, Q-CCSD and E-CCSD remain the most satisfying methods, 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%. It is well-known that CC2 generally performs better for valence than Rydberg transitions.⁹²,⁹³ The performance of all TD-DFT methods is 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.

As a further analysis we investigated the results for valence states with large $\delta^{\text{TPA}}$ (larger than 20 a.u. at the CC3/aug-cc-pVDZ level), since they are the more interesting for actual applications. We underline that this subset includes only 9 transitions, so these results should be interpreted with caution. The results of MAPE and $R^2$ are reported in Fig. 4, while the other indicators can be found in Tables S11–S15 in the SI. Q-CCSD is again the most adequate method with a MAPE of 24.4% and a $R^2$ of 0.93. Interestingly, for this subset, Q-CC2 and E-CCSD show similar performances: MAPE of 34.3% (44.8%), and $R^2$ of 0.88 (0.86) for CC2 (E-CCSD). We underline that the performances of Q-CCSD, E-CCSD and Q-CC2 are mainly driven by the $1B_{3g}$ state of tetrazine.
<figure><img src="image_9.png" /><figcaption>Figure 3: MAPE and R² of δ<sup>TPA</sup> for all the wavefunction and TD-DFT methods for the sets for all the singly-excited states for different δ<sup>TPA</sup> thresholds (in a.u.). References are at the Q-CC3/aug-cc-pVDZ level. The values of all the indicators are listed in Tables S11–S13 in the SI.</figcaption></figure>

In fact its Q-CC3 δ<sup>TPA</sup> of 56.0 a.u., is largely overestimated by both Q-CCSD (86.6 a.u.) and E-CCSD (118.3 a.u.), while CC2 provides a rather accurate estimate of 47.7 a.u. Indeed if we exclude this molecule the MAPE becomes: 21.1% for Q-CCSD, 37.4% for E-CCSD and 37.1% for Q-CC2. I-ADC2 still performs quite poorly with a MAPE of 96.7% and MAE of 43.1, but has a relatively acceptable R² of 0.76. TD-DFT methods perform similarly as in the whole set in terms of MAPE but they show extremely poor R² values. 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 for actual TD-DFT applications for TPA simulations.

<figure><img src="image_10.png" /><figcaption>Figure 4: MAPE and R² of δ<sup>TPA</sup> for all the wavefunction and TD-DFT methods for the set of valence states showing 20 ≤ δ<sup>TPA</sup><sub>Q-CC3</sub> [a.u.]. References are obtained at the Q-CC3/aug-cc-pVDZ level. The values of the indicators are reported in Tables S11–S15 in the SI. The legend is reported in Fig. 2</figcaption></figure>

The performances of the wavefunction methods Q-CCSD, E-CCSD, Q-CC2, I-ADC2, and DEA-CCSD (which raw values are reported in Table S6 in the SI) have been evaluated for the six electronic states with significant contributions from double excitations. TD-DFT methods were excluded due to their inability to account for double excitations, and any ac-

curate reproduction of δ<sup>TPA</sup> by TD-DFT would likely be due to significant error cancellation. Table 5 reports the absolute percentage error (APE) for these methods compared to reference Q-CC3 calculations. Values with APEs exceeding 200% are not detailed.

<table><caption>**Table 5:** Absolute Percentage Error (APE) for δ<sup>TPA</sup> of molecules with double excitations. Reference values obtained at the Q-CC3/aug-cc-pVDZ level. CycloP stands for Cyclopentadiene.</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>Acrolein</td><td>2A"</td><td>65</td><td>79</td><td>&gt;200</td><td>115</td><td>96</td></tr><tr><td></td><td>5A'</td><td>159</td><td>&gt;200</td><td>&gt;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>&gt;200</td><td>&gt;200</td></tr><tr><td>CycloP</td><td>1A<sub>1</sub></td><td>96</td><td>129</td><td>15</td><td>&gt;200</td><td>&gt;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>&gt;200</td><td>&gt;200</td><td>62</td><td>&gt;200</td><td>&gt;200</td></tr></tbody></table>

The DEA-CCSD method provides notably accurate results for the 1A₁ state of cyclopentadiene and the 1A_g 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_g state of hexatriene, DEA-CCSD is the only method that yields a reasonably accurate result (APE of 62%), whereas all other methods provide very poor values. On the other hand, for the two states of acrolein DEA-CCSD appears quite unreliable, indeed we have an APE larger 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<sub>2g</sub> state of benzene, with E-CCSD showing the lowest APE (4%), and DEA-CCSD has the highest APE (63%).

In summary, Q-CCSD generally provides acceptable results for states with not too large double excitation character, but fails for the 1A_g 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.⁸¹ As expected, methods with perturbative doubles are significantly less reliable.

## CONCLUSIONS

This study evaluated the performance of several computational methods and basis sets for calculating the vertical two-photon transition strengths ($\delta^{TPA}$) of singlet excited-state in small- and medium-sized molecules, using Q-CC3 as a benchmark. 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) reduce the performances up to 24%, thus advising 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. 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 XCF, 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 $\delta^{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 then $O(N^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.⁶⁴

## ACKNOWLEDGMENTS

This work was supported by the French Agence Nationale de la Recherche (ANR) under Contract No. ANR-21-CE07-0058-2 (CONDOR). This research used resources of the GLiCID Computing Facility (Ligerien Group for Intensive Distributed Computing, https://doi.org/10.60487/glicid, Pays de la Loire, France). R.Z. acknowledges financial support from

National Science Centre (grant no. 2018/30/E/ST4/00457).

## SUPPORTING INFORMATION AVAILABLE

Supporting informations contains: 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.

## 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₂: Eu²⁺. *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. *Journal of 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. *Analytical Chemistry* **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. Multi-photon 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. In *Handbook of Computational Chemistry*; Leszczynski, J., Ed.; Springer Netherlands: Dordrecht, 2016; pp 1–19.

(18) Meyers, F.; Marder, S.; Pierce, B.; Brédas, J. Electric Field Modulated Nonlinear Optical Properties of Donor-acceptor Polyenes: Sum-over-states Investigation of the Relationship Between Molecular Polarizabilities (*α*,*β*, and *γ*) 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. et al. 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*.

(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–95.

(28) Paterson, M. J.; Christiansen, O.; Pawłowski, F.; Jørgensen, P.; Hättig, C.; Helgaker, T.; Sąłek, 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. et al. 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. *WIREs 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) 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.

(39) 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.

(40) 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.

(41) 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.

(42) 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.

(43) 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.

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

(45) 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.

(46) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P. et al. The Dalton Quantum Chemistry Program System. *WIREs Comput. Mol. Sci.* **2014**, *4*, 269–284.

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

(48) 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.

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

(50) 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.

(51) Turbomole V7.2 2017, a Development of university of Karlsruhe and forschungszentrum Karlsruhe Gmbh, 1989-2007, turbomole Gmbh, Since 2007; Available From http://www.turbomole.com.

(52) 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. et al. Turbomole: Modular Program Suite for Ab Initio Quantum-chemical and Condensed-matter Simulations. *J. Chem. Phys.* **2020**, *152*, 184107.

(53) 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.

(54) 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.

(55) Petrusevich, E. F.; Głodek, M.; Antoniak, M. A.; Muzioł, 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 A* **2023**, *295*, 122600.

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

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

(58) Epifanovsky, E.; Gilbert, A. T. B.; Feng, X.; Lee, J.; Mao, Y.; Mardirossian, N.; Pokhilko, P.; White, A. F.; Coons, M. P.; Dempwolf, A. L. et al. software for the Frontiers of Quantum Chemistry: An Overview of Developments in the Q-chem 5 Package. *J. Chem. Phys.* **2021**, *155*, 084801.

(59) 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.

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

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

(62) 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.

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

(64) Rauwolf, N.; Klopper, W.; Holzer, C. non-linear Light-matter Interactions From the Bethe-salpeter Equation. *J. Chem. Phys.* **2024**, *160*, 061101.

(65) Sąłek, 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.

(66) Parker, S. M.; Rappoport, D.; Furche, F. Quadratic Response Properties From TDDFT: Trials and Tribulations. *J. Chem. Theory Comput.* **2018**, *14*, 807–819.
(67) 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. et al. Recent Developments in the General Atomic and Molecular Electronic Structure System. *J. Chem. Phys.* **2020**, *152*, 154102.

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

(69) 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.

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

(71) Chołuj, 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.

(72) 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.

(73) 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.

(74) 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*, 5032–5051.

(75) 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.

(76) 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.

(77) 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.

(78) 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.

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

(80) 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.

(81) 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.

(82) 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.

(83) 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.

(84) 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.

(85) 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.

(86) 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.

(87) 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.

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

(89) 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.

(90) 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.

(91) Wielgus, M.; Bartkowiak, W.; Samoc, M. Two-photon Solvatochromism. I. Solvent Effects on Two-photon Absorption Cross Section of 4-dimethylamino-4'-nitrostilbene (dans). *Chem. Phys. Lett.* **2012**, *554*, 113–116.

(92) 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.
(93) 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. *WIREs Computational Molecular Science* **2021**, 11, e1517.