Benchmarking Two-photon Absorption Cross Sections: Performance of Cc2 and Cam-b3lyp

a We investigate the performance of CC2 and TDDFT/CAM-B3LYP for the calculation of two-photon absorption (TPA) strengths and cross sections and contrast our results to a recent coupled cluster equation-of-motion (EOM-EE-CCSD) benchmark study [K. In particular, we investigate whether CC2 TPA strengths are significantly overestimated compared to higher-level coupled-cluster calculations for fluorescent protein chromophores. Our conclusion is that CC2 TPA strengths are only slightly overestimated compared to the reference EOM-EE-CCSD results and that previously published overestimated cross sections are a result of inconsistencies in the conversion of the TPA strengths to macroscopic units. TDDFT/CAM-B3LYP TPA strengths, on the other hand, are found to be 1.5 to 3 times smaller than the coupled-cluster reference for the molecular systems considered. The unsatisfactory performance of TDDFT/CAM-B3LYP can be linked to an underestimation of excited-state dipole moments predicted by TDDFT/CAM-B3LYP.


Introduction
Two-photon absorption (TPA)-that is, a simultaneous absorption of two photons 1 -is a non-linear optical process with an intensity depending on the square of the incoming light. Searching for molecules with high TPA cross sections is an active field of research 2,3 and the interplay between experiment and theory is paramount for a good understanding of structure-property relations. 2,4 Quantum chemical calculation of TPA strengths has been implemented at various levels of theory in different software packages and has recently gained significant interest. The implementation in the Dalton software package 5 was described by Hättig, Christiansen and Jørgensen 6 for coupled-cluster methods and by Sałek et al. 7 for time-dependent density functional theory (TDDFT). More recently, Friese, Hättig and Ruud 8 implemented TPA based on the resolution-of-identity CC2 (RI-CC2) in a development version of Turbomole, 9 thereby enabling the calculation of TPA at the coupled-cluster level for medium-sized molecules. Further, calculation of TPA strengths was implemented in QChem 10,11 for the algebraic diagrammatic construction method by Knippenberg et al. 12 and very recently for equation-of-motion for excitation energies CCSD (EOM-EE-CCSD) by Nanda and Krylov. 13 In addition, Salem and Brown 14 have reported TDDFT TPA cross sections calculated in GAMESS. 15 All these different implementations span the spectrum from computationally efficient to very accurate, with TDDFT allowing for efficient calculations on large molecules and the coupled-cluster hierarchy of methods allowing for systematic improvement of the accuracy. CC2 and TDDFT have so far been benchmarked against higher-order coupled-cluster methods only for very small organic molecules of up to six atoms. 16 Benchmarking CC2 and TDDFT TPA strengths against higher-level methods is important in order to assess the accuracy of these methods in work on larger molecules, for which more accurate calculations are not yet feasible. The largest molecule studied using EOM-EE-CCSD by Nanda and Krylov is 4-(p-hydroxybenzylidene)-1,2-dimethylimidazolin-5-one (HBDI) with 28 atoms. 13 CC2 TPA cross sections have been reported for molecular systems of more than twice that size: a truncated molecular tweezer complex consisting of 78 atoms 8 and the YFP-Tyr203 dimer with 62 atoms. 17 TDDFT TPA cross sections have been reported for molecular systems with more than 100 atoms, such as the fullerene-buckycatcher complex in ref. 18.
Very recently, Nanda and Krylov reported TPA strengths for medium-sized (20 to 28 atoms) neutral fluorescent protein chromophores using the EOM-EE-CCSD method. 13 By comparing their calculated cross sections to published CC2 results using the same geometry of the HBI chromophore, 14 they noted a strong overestimation (a factor of 4.4) by CC2. 13 The suggested explanations for this discrepancy were a strong overestimation of TPA strengths for CC2 or a mistake in the conversion to macroscopic units. 13 In addition, a previous study by three of us on CC2 and TDDFT/CAM-B3LYP cross sections on intermolecular chargetransfer excitation in the yellow fluorescent protein also reported CC2 TPA cross sections to be larger than those calculated with TDDFT/CAM-B3LYP. 17 The aim of this study is to assess the quality of CC2 and TDDFT/CAM-B3LYP TPA strengths by comparing with the recent benchmark study by Nanda and Krylov 13 and to inspect whether CC2 calculations do indeed strongly overestimate the TPA strength. Comparison between different methods are often made on the basis of macroscopic rather than microscopic units, leading to extra challenges and complications for benchmarking work. Therefore we begin in Section 2 by addressing different factors that need to be considered in order to make a correct comparison between TPA cross sections in macroscopic units. Following this, we compare the TPA strengths of a set of fluorescent protein chromophores computed with CC2 and TDDFT/CAM-B3LYP to those obtained with EOM-EE-CCSD 13 and previously reported CC2 14 cross sections, using the same molecular geometries as in the reference works.

Conversion to macroscopic units
The macroscopic TPA cross section s TPA in cgs units can be obtained from the rotationally averaged TPA strength hd TPA i in atomic units using 19 where N is an integer value (see Section 2.3), a is the fine structure constant, a 0 the Bohr radius, o the photon energy in atomic units, c the speed of light and g(2o, o 0 , G) the lineshape function describing spectral broadening effects (see Section 2.2). The common unit for TPA cross sections is GM (originally introduced as 'maria') 20 in honour of the work of Maria Göppert-Mayer, 1 with 1 GM corresponding to 10 À50 cm 4 s photon À1 . When deriving eqn (1) without taking into account factors related to the comparison to a specific experimental setup, one obtains N = 8. 19,21 As will be elaborated in Section 2.3, this corresponds to a so-called double-beam experiment. In passing, we note that the original work of Peticolas, describing the conversion of multiphoton absorption strengths to macroscopic units, contains a typo in the equation for the transition strength for twophoton processes (eqn (10) in that work). 19 In fact, the exponent 2 of the term 2p h/V is omitted in eqn (10), which could make one believe that eqn (1) depends on 4p 2 rather than 8p 3 for a double-beam experiment. It is clear that this is a typo by comparing that eqn (10) to eqn (17) on two-photon absorption and eqn (28) on three-photon absorption in the same work. 19 This indicates that a factor (2p) 3 should appear in the final conversion for two-photon absorption and a factor (2p) 4 for three-photon absorption. This factor (2p) j+1 for j-photon absorption is also derived in ref. 21.
One should be aware that o in eqn (1) is the photon energy, which in the degenerate TPA case corresponds to half the excitation energy. Defining o as the excitation energy (2o) rather than the photon energy (as is done in ref. 14 and 22) results in TPA cross sections higher by a factor of 4.
In the following we will describe the prefactor from the rotational averaging of the TPA strength (Section 2.1), the lineshape function g(2o, o 0 , G) (Section 2.2), and the type of experiment that is compared to (Section 2.3) in more detail with particular emphasis on what one needs to be aware of when comparing TPA strengths in macroscopic units.

The prefactor of the rotational averaging
The rotationally averaged TPA strength hd TPA i in atomic units can be obtained from the TPA transition moment S as 23-25 where the sum is over the Cartesian components a and b and with the bar indicating complex conjugation. In the electric dipole approximation, the (a, b)'th component of the TPA transition moment S if between the initial state i and final state f is defined as 26 where hi|m a |ni is the transition dipole moment between the electronic states i and n along the Cartesian axis a, o ni the associated excitation energy and o 1 and o 2 the frequencies of photon 1 and 2. m a = m a À h0|m a |0i, i.e.hn| m a | f i is the difference dipole moment between the ground and excited state if n = f 27 and a transition dipole moment if n a f.
For linearly polarized light with parallel polarization (F = G = H = 2) 23 and two photons of the same energy, eqn (2) reduces to Alternatively, hd TPA i can be defined as 28 or as 14,22 in which case the prefactor of the rotational averaging is included in the conversion factor as 14,22,28,29 One should however be aware that the same eqn (7) results from the use of N = 4 in ref. 14 and 22 (where hd TPA i is defined as in eqn (6), leaving out a factor of 1/15) and from the use of N = 8 in ref. 28 (where hd TPA i is defined as in eqn (5), leaving out a factor of 1/30). Indeed, it is unclear whether N = 4 or N = 8 is used in ref. 29 because the conversion is defined as in eqn (7) without a definition of hd TPA i. This points to the importance of publishing not only the equation used for the conversion factor (e.g. eqn (1)), but also the way in which the rotational average is defined (e.g. eqn (4)).

The lineshape function
The lineshape function g(2o, o 0 , G) with broadening factor G describes broadening effects, correcting for the infinitely sharp calculated vertical excitations and allowing for comparison to experimental peaks, in which rovibrational excitations and collisional dynamics also play a role. 26 This is usually done in a phenomenological way by introducing an empirical damping parameter G to describe the spectral broadening of an excitation. In many gas-phase calculations, a Lorentzian lineshape function is used for g(2o, o 0 , G), with o the photon energy, o 0 the excitation energy and G the half width at half maximum (HWHM). The function is normalized to one. By setting o = o 0 /2, the maximum of the Lorentzian is obtained as When this maximum is inserted for g(2o, o 0 , G) in eqn (1), one obtains for the TPA cross section at the calculated excitation energy with hd TPA i defined as in eqn (2) or (4). We note that the exponent of p is 3 in eqn (1) and 2 in eqn (10), which is sometimes overseen when the equation is printed. 22 In solution, a Gaussian lineshape function is most commonly used for g(2o, o 0 , G), 30 again normalized to one and with o the photon energy, o 0 the excitation energy and G the HWHM. The maximum of the Gaussian is obtained at o = o 0 /2 as When this maximum is inserted for g(2o, o 0 , G) in eqn (1), one obtains Thus, when using the same broadening factor G, the maxima of the Lorentzian and Gaussian broadening functions have a different value. The Lorentzian function has a broader base and the Gaussian function has a higher maximum by a factor of ffiffiffiffiffiffiffiffiffiffiffi p ln 2 p % 1:48, see Fig. 1a.
One should be aware that G can also be defined as the full width at half maximum (FWHM) instead of the HWHM (FWHM = 2ÁHWHM). The use of G = 0.2 eV as the FWHM is equivalent to a HWHM of 0.1 eV when one substitutes G for G/2 31-34 in eqn (8) or (11). If one uses 0.1 eV as FWHM, 35 the resulting TPA cross section s TPA is twice as high compared to work using 0.1 eV as HWHM, see Fig. 1b and eqn (9) and (12). We also point out that the use of eqn (10) with N = 8 18,[36][37][38] can in principle result from either the use of eqn (1) with N = 8 and G defined as the HWHM or with N = 4 and G defined as the FWHM.
The broadening effects are different for each excited state, 26 which can also be taken into account in theoretical work. 30,32,36,39 What is usually done, however, is choosing a single empirical parameter for G, often chosen to be 0.1 eV, 13,14,18,21,22,28,29,37,38 but also other broadening factors have been used, usually taken from the broadening of a specific peak in an experimental spectrum. 2,31,34,35 We note that an alternative approach used for the calculation of TPA strengths is given by damped response theory, 40 in which an imaginary factor ig is added to the optical frequency o to incorporate broadening effects. In work on damped response theory with the complex polarization propagator such as ref. 40, g is used for the HWHM and G for the FWHM.
One should be cautious when comparing TPA cross sections computed with different broadening factors as they affect also the maximum of the peak, see Fig. 1c. The discussion here makes it sufficiently clear that giving the value of G and specifying the type of broadening and whether G is HWHM or FWHM is important for the reproducibility of published TPA cross sections.

Comparison with experiment
The use of different integers N in eqn (1) is among other reasons due to comparison to different experimental setups. When deriving eqn (1), one obtains N = 8. 19,21 This corresponds to an experimental setup with two laser sources (a so-called doublebeam experiment 41 ), thus in principle allowing for two photons with different energies. The use of N = 8 in eqn (1), however, is not correct when comparing calculated cross sections to the photon dissipation rate in single-beam TPA experiments. [30][31][32][33]39,42,43 In fact, the photon dissipation rate is in this case twice the calculated twophoton absorption rate because two photons together excite the molecule. Thus, one needs N = 16 in eqn (1) to compare to singlebeam experiments. 31,33,39,42,43 When two lasers are used as the photon sources-as intended in e.g. ref. 19 and 21-the photon dissipation rate of each of the lasers equals the two-photon absorption rate in the sample.
Furthermore, the TPA transition moment (see eqn (3)) can be defined in different ways. [31][32][33]39,[41][42][43] The correct definition of the two-photon transition moment S if to compare to singlebeam experiments (such as the z-scan technique 3 ) is half the one in eqn (3), [31][32][33]39,[41][42][43] where there is only one photon with frequency o. The use of eqn (14) instead of eqn (3) leads to hd TPA i-which depends on terms of S times % S-lower by a factor of 4. To correct for this factor of 4 in hd TPA i, a factor of N = 4 (instead of N = 16) in the conversion to macroscopic units has been used in combination with hd TPA i calculated using S in eqn (3) to obtain a s TPA that can be compared to single-beam experiments. 30,32,34 This is the strategy implemented in e.g. the Dalton program,5,44 where S is defined as in eqn (3). The choice of the right value of N in eqn (1) is summarized in Fig. 2.

Computational details
We present CC2 and TDDFT/CAM-B3LYP one-and two-photon absorption calculations for the four neutral fluorescent protein chromophores in Fig. 3. The HBI (4-( p-hydroxybenzylidene)imidazolidin-5-one) and HBDI (4-( p-hydroxybenzylidene)-1,2dimethylimidazolin-5-one) molecules are models for the neutral green fluorescent protein chromophore. The chromophore of the cyan fluorescent protein (CFP) is similar to HBI with an indole ring in place of the phenol ring. The chromophore of the phosphorescent yellow protein (PYP) is p-coumaric acid. The chromophores of fluorescent proteins are frequently studied using TPA both experimentally 45,46 and computationally 13,14,17,22,37,38,47,48 due to their relevance in biological imaging.
All calculations are gas-phase calculations on vertical excitations. Calculations with the approximate coupled-cluster singles and doubles model (CC2) 49 were done in Turbomole 9 using the resolution-of-identity (RI) approximation. 8,50,51 The RI approximation has been shown to give insignificant errors compared to the basis-set error for excitation energies 50 and one-photon oscillator strengths. 51 1s-orbitals of non-hydrogen atoms were kept frozen in the correlation treatment. TDDFT calculations were done using Dalton 5,7,44 with the CAM-B3LYP exchangecorrelation functional 52 with a = 0.19, b = 0.46 and m = 0.33. The range-separated functional CAM-B3LYP has an average error that is independent of the degree of charge transfer in the excitations investigated. 29,47,53 Furthermore, it yields excitation energies 47 and oscillator strengths, 54 also considered in this work, in very good agreement with coupled-cluster results for neutral chromophores. It should be noted that the parameters of the functional are optimized using a test set of only neutral chromophores, 52 making the transferability to charged systems unknown.
Basis sets from the Pople (6-31+G*) 55-57 and Dunning (aug-cc-pVDZ, d-aug-cc-pVDZ and aug-cc-pVTZ) 58 families were used to ensure reliable comparison to previous works. 13,14 The auxiliary basis sets aug-cc-pVDZ (for molecular basis sets 6-31+G*, aug-cc-pVDZ and d-aug-cc-pVDZ) and aug-cc-pVTZ (for molecular basis set aug-cc-pVTZ) were used for all RI-CC2 calculations. 59 To compare to the work of Nanda and Krylov, 13 we removed the d-functions from the second set of augmentation functions from the d-aug-cc-pVDZ basis set. We will refer to this modified basis set as mod-d-aug-cc-pVDZ in Section 4.
Molecular structures were taken from reference works to ensure that differences in calculated one-photon and two-photon absorption properties result only from different methods and not from using different molecular geometries. We will therefore only compare our results to other work where exactly the same molecular geometries were used 13,14 and not to works where different geometries of the same molecules were studied. 22,37,47,48 The geometries of HBI and CFP were taken from Salem and Brown, 14 where the geometries were optimized with DFT using the PBE0 exchangecorrelation functional and the 6-31+G** basis set. The geometries of  HBDI and PYP were taken from Nanda and Krylov, 13 where the geometries were optimized with RI-MP2 and the aug-cc-pVDZ and aug-cc-pVTZ basis sets, respectively. All chromophores have C s symmetry, which was also exploited in the CC2 calculations. We limit our study to the lowest excitation of the chromophores (lowest two excitations for PYP), which is of A 0 symmetry and pp* character. The excitation to the 2A 0 state of PYP and the considered excitations in HBI, HBDI and CFP are dominated by the HOMO-LUMO transition with both orbitals involved delocalized over the whole conjugated system of the chromophore. The 1A 0 excitation of PYP has contributions from the HOMO and LUMO as well as from an occupied p and a virtual p* orbital that are both located on the phenol ring. Rotational averaging of the TPA transition moments was done using eqn (4) with S ab defined as in eqn (3). Conversion to GM units was done using eqn (1) with N = 4 and a Lorentzian lineshape function with a HWHM of 0.1 eV was used to ensure comparability to reference work. 13 Since we only report cross sections at resonant conditions (o = o 0 /2), the conversion corresponds to eqn (10) with N = 4 and G = 0.1 eV.

Results
Calculated excitation energies and one-and two-photon absorption properties for the four molecules in Fig. 3 are tabulated in Tables 1-4 (HBI, HBDI, CFP and PYP) together with reference values from the literature. 13,14 The calculated CC2/6-31+G* TPA cross sections s TPA for HBI (Table 1) and CFP (Table 3) are a factor of 4 lower than the CC2 cross section in ref. 14 using the same molecular geometry. This shows that the conversion factor used here and by Nanda and Krylov 13 is different from the one used by Salem and Brown 14 by a factor of 4. The overestimation of the TPA cross section of the first A 0 state of the HBI chromophore by CC2/6-31+G* in comparison to EOM-EE-CCSD/ 6-31+G* is thus not a factor of 4.4, 13 but a factor of 1.1. We find that CC2 also overestimates EOM-EE-CCSD TPA strengths hd TPA i with the mod-d-aug-cc-pVDZ basis set by a factor of 1.16 for HBDI and by factors of 1.25 and 1.34 for the 2A 0 and 1A 0 states in PYP, respectively (see Tables 2 and 4). Differences in excitation energies, oscillator strengths and TPA cross sections between the aug-cc-pVDZ and mod-d-aug-cc-pVDZ basis sets calculated with CC2 are negligible (see Tables 1, 2

and 4).
TDDFT/CAM-B3LYP TPA strengths hd TPA i, on the contrary, are much lower than the coupled-cluster results. The difference depends on the molecule and is for different basis sets a factor of 3.4 to 3.6 for HBI, 2.3 to 2.6 for HBDI, 3.1 to 3.2 for CFP and 1.9 for PYP compared to CC2 calculations. The underestimation of CAM-B3LYP/aug-cc-pVDZ in comparison to EOM-EE-CCSD/d-aug-cc-pVDZ is a factor of 2.00 for HBDI and 1.50 for PYP for hd TPA i and a factor of 3.00 for HBI for s TPA (for which no hd TPA i is given in ref. 13). The TDDFT/CAM-B3LYP cross sections are more or less converged when using the aug-cc-pVDZ basis set, as was also shown in a previous study of oneand two-photon properties on the neutral and anionic structures of the GFP chromophore using the crystal structure geometry. 37 TPA cross sections are slightly larger with the 6-31+G* basis set, which was previously found to be the case for the neutral but not for the anionic GFP chromophore. 37 In contrast to the TPA strengths, the CAM-B3LYP oscillator strengths are in good agreement with the coupled-cluster results with very similar values for HBDI, a slight underestimation for HBI and CFP and a modest overestimation for PYP, none of which exceeds 20%. This is in agreement with a more thorough analysis of Caricato et al., who found that CAM-B3LYP provides the best agreement with EOM-CCSD oscillator strengths among the set of density functionals investigated for 11 small organic neutral molecules. 54

Discussion
The most important observation from the numerical results presented in the previous section is that CC2 TPA strengths agree well with the benchmark EOM-EE-CCSD results, whereas  the TDDFT/CAM-B3LYP TPA strengths are significantly underestimated.
In a thorough comparison of TPA strengths of small molecules of up to six atoms by Paterson et al., CC2 was compared to the higher-order methods CCSD and CC3. 16 For some of the investigated transitions, CC2 was observed to overestimate CCSD and CC3 TPA strengths, but not by more than a factor of 2. CCSD, on the other hand, was shown on average to overestimate TPA strengths slightly in comparison to CC3. 16 Combined with the results presented in this work, this seems to indicate that CC2 performs generally well in the calculation of TPA strengths with the possible exception of specific transitions that require an accurate description of double excitations.
It is interesting to understand whether the underestimation of TDDFT/CAM-B3LYP TPA strengths is specific to the limited number of excitations in the (neutral) molecules investigated here, or more general. While the number of studies reporting TPA cross sections calculated with TDDFT is high, few of them address the accuracy compared to more accurate methods. This is in part related to the lack of efficient publicly available codes to compute coupled-cluster TPA strengths. At least four previous studies have included a comparison of CC2 and TDDFT/ CAM-B3LYP TPA strengths. Paterson et al. 16 included different density functionals in their comparison to the coupled-cluster hierarchy of methods on small molecules of up to six atoms. They found that CAM-B3LYP generally outperforms LDA, BLYP and B3LYP but overestimates most of the TPA strengths. 16 Friese, Ruud and Hättig discuss the agreement between CC2 and TDDFT/CAM-B3LYP mainly in a qualitative way, but report cross sections of CAM-B3LYP/cc-pVDZ both higher and lower than CC2/cc-pVDZ for qualitatively similar electronic transitions. 8 In a previous work, three of us have inspected a p-stacking system of an extended anionic HBDI chromophore and a tyrosine residue found in the yellow fluorescent protein. 17 CC2/aug-cc-pVDZ TPA cross sections were found to be much higher than CAM-B3LYP/ aug-cc-pVDZ cross sections. In particular, the TPA cross section of the pp* excitation located on the chromophore was lower by a factor of 6.3 for the CAM-B3LYP calculations. Comparing this with the factor of 2.3 in this work (Table 2), we note a significant difference that is presumably due to using the neutral or anionic HBDI chromophore. Salem and Brown 14 report CC2/6-31+G* TPA cross sections for the same pp* transition of five neutral chromophores, including HBI and CFP. Their CAM-B3LYP/ 6-31+G** gas-phase TPA cross sections (in their ESI) are smaller by a factor of 1.3 for anionic HBI, 3.6 for neutral HBI, 3.3 for CFP and 2.5 and 4.5 for two different protonation states of the chromophore of the blue fluorescent protein. Similar or larger underestimations are reported for the B3LYP and PBE0 exchangecorrelation functionals.
The underestimation of CAM-B3LYP TPA strengths can be understood by its underestimation of excited-state dipole moments. [60][61][62] Eriksen et al. found that the ground-state dipole moment of para-nitroaniline is overestimated by CAM-B3LYP, while the excited-state dipole moment of the charge-transfer state is dramatically underestimated, resulting in a much too small difference dipole moment. 60 Zaleśny and co-workers found that the problems with excited-state dipole moments are not limited to a specific molecule or a particular functional, but constitute a general problem for TDDFT. 61,62 An underestimation of the excited-state dipole moment-combined with accurate or even overestimated ground-state dipole momentsleads to underestimated difference dipole moments. 62 The difference dipole moments enter in the expression for the TPA transition moment in eqn (3) through the term hn| m| f i with n = f, which is often the dominating term for noncentrosymmetric molecules. 63 Thus, underestimated difference dipole moments lead to underestimated TPA strengths. For n a f the element hn| m| f i in eqn (3) is a transition dipole moment between two excited states. These elements are likely to be inaccurate for TDDFT as well since they crucially depend on the electron density of the excited states, but they are less studied in the literature. Concluding the discussion of the TDDFT/CAM-B3LYP TPA strengths, we can say that TDDFT/ CAM-B3LYP generally seems to underestimate TPA strengths for these molecules with the magnitude varying for different functionals, molecules and protonation states. There seems to be a need for benchmarking TDDFT/CAM-B3LYP TPA strengths of a larger set of medium-sized chromophores-including different protonation states-against a more accurate method such as EOM-EE-CCSD in order to obtain a better estimation of the error inherent to using TDDFT/CAM-B3LYP to calculate the TPA strength. Table 4 Excitation energy 2o and DE, oscillator strength f, TPA strength hd TPA i and TPA cross section s TPA for the 1A 0 and 2A 0 states (1A 0 state for CAM-B3LYP) of p-coumaric acid (PYP chromophore). The states are labeled to correspond to the order of the EOM-EE-CCSD states in ref. 13 and are therefore not in order of ascending energy for the CC2 calculations. All calculations are performed using the same molecular geometry taken from ref. In this study we have limited the calculations to neutral chromophores of fluorescent proteins that allow for comparison to available benchmark data. We have not tried to reproduce experimental TPA spectra. To do this in a proper way, one should not only take into account the form of the lineshape function for a particular excitation 32,36,39 (see Section 2.2) and the type of experiment that it is compared to (see Section 2.3), but also temperature effects, 48 non-Condon transitions 64 and the environment. 37,38 In particular the latter has been shown to affect both experimental 45,46 and calculated 37,38 TPA cross sections in a profound way. Indeed, it has been shown that different fluorescent proteins with the same chromophore can have maxima in the TPA spectrum that differ by a factor of 5 in intensity. 45 The discussion of the conversion factor in Section 2 clearly shows the importance of providing the (correct) details about the conversion used for determining the TPA cross sections from TPA strengths in atomic units. The factor of 4 difference in TPA cross sections between Salem and Brown 14 and Nanda and Krylov 13 is probably related to the former using excitation energies instead of photon energies for o. For the sake of reproducibility, a good practice is to provide excitation energies and TPA cross sections in atomic units in addition to macroscopic units, as we have done in Tables 1-4. Not only should the correct formulae for rotational averaging and conversion to GM units be given explicitly and applied consistently, one should also state the type of broadening used and whether the broadening factor G is the full or half width at half maximum. When comparing to experiments, one needs to ensure that the right conversion factor is used for the type of experiment employed (see Fig. 2). [31][32][33]42,43 When comparing to experiments that use the relative fluorescence method, 4,34,46 one could consider computing the TPA cross section for the reference molecules 65 as well and compare the relative cross sections to relative fluorescence intensities from experiment. This approach corrects for systematic errors in the computational and experimental method in much the same way as for chemical shifts in NMR spectroscopy, which can be more reliably compared to experiment than absolute shielding constants. 66 In summary, care should be exercised when comparing calculated and experimental TPA cross sections, whereas for comparison between computed results it suffices to use atomic units or to ensure that the same conversion protocol is used as in the reference work.

Conclusion
We have compared CC2 and TDDFT/CAM-B3LYP TPA strengths against recent EOM-EE-CCSD benchmark data for a set of neutral fluorescent protein chromophores. This provides a much-needed comparison of two methods that allow TPA calculations on larger (biological) systems (CC2 and TDDFT/ CAM-B3LYP) against a higher-order method (EOM-EE-CCSD). We have found that CC2 TPA strengths are in good agreement with the EOM-EE-CCSD cross sections with a small overestimation up to a factor of 1.4 that depends on the molecule. TDDFT/CAM-B3LYP TPA cross sections, however, are underestimated compared to coupled-cluster results by a factor of 1.5 to 3 for this set of molecules, whereas the one-photon oscillator strengths do not deviate by more than 20%. The underestimation of TDDFT/CAM-B3LYP TPA cross sections can be understood in terms of the known underestimation of excited-state dipole moments of TDDFT methods.
We have demonstrated how an erroneous assumption on CC2 TPA cross sections was based on a conversion of the TPA strengths to macroscopic units that was different from the one used in the reference work. This stresses the importance of providing all relevant details in the conversion to macroscopic units in computational work on TPA, namely (a) the excitation energy and TPA strength in atomic units, (b) the formulae for both the conversion and the rotational averaging of the transition moments, (c) the proper choice of the conversion for comparison to different types of experiments (see Fig. 2) and (d) the type of lineshape function used, including the broadening factor and whether it corresponds to the full or half width at half maximum.