# Simulations of Two-Photon Absorption Spectra of Fluorescent Dyes: The Impact of Non-Condon Effects

Rudraditya Sarkar, Carmelo Naim, Karan Ahmadzadeh, Robert Zaleśny, Denis Jacquemin, Josep M Luis

► To cite this version:

Rudraditya Sarkar, Carmelo Naim, Karan Ahmadzadeh, Robert Zaleśny, Denis Jacquemin, et al. Simulations of Two-Photon Absorption Spectra of Fluorescent Dyes: The Impact of Non-Condon Effects. Journal of Chemical Theory and Computation, 2025, 21 (7), pp.3587-3599. 10.1021/acs.jctc.4c01545. hal-05117495

HAL Id: hal-05117495

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

Submitted on 17 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Simulations of Two-Photon Absorption Spectra of Fluorescent Dyes: The Impact of Non-Condon Effects

Rudraditya Sarkar,¶ Carmelo Naim,¶ Karan Ahmadzadeh, Robert Zaleśny,* Denis Jacquemin,* and Josep M. Luis*

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

Cite This: J. Chem. Theory Comput. 2025, 21, 3587–3599

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

Read Online

## ACCESS

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

Metrics & More

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

Article Recommendations

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

Supporting Information

**ABSTRACT:** Computer simulations play a pivotal role in interpreting experimental two-photon absorption (2PA) spectra. One of the key aspects of the simulation of these spectra is to take into account the vibrational fine structure of the bands in electronic spectra. This is typically done by employing Franck–Condon (FC) term and low-order terms in the Herzberg–Teller (HT) expansion. In this work, we present a systematic study of first-order HT effects on the vibronic structure of $\pi \rightarrow \pi^*$ electronic bands in 2PA spectra of 13 common fluorophores. We begin by evaluating the performance of several density functional approximations (DFAs) against the second-order coupled cluster singles and doubles model (CC2) for reproducing two-photon transition moments and their first- and second-order derivatives with respect to normal modes of vibration on a set of six donor–acceptor molecules. Our findings reveal that most DFAs produce inaccurate values for these derivatives, with the exception of the LC-BLYP functionals with range–separation parameters of 0.33 and 0.47. Although these functionals underestimate the HT contribution to the 2PA total intensities of the $\pi \rightarrow \pi^*$ electronic bands, they offer a reasonable qualitative reproduction of the HT vibrational fine structure of the reference spectra. We further explore HT effects on fluorescent chromophores, finding that HT contributions are secondary to FC effects, leading to small shifts of the wavelengths peaks, and minimal changes in the intensities. Additionally, the adiabatic Hessian, vertical Hessian, and vertical gradient vibronic models are assessed. The general agreement among these models confirms that the harmonic approximation is suitable for studying the selected fluorophores.

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

## 1. INTRODUCTION

Multiphoton absorption phenomena, with two-photon absorption (2PA) at their forefront, are currently in the limelight of both experimental and computational studies. The main rationale explaining this interest in 2PA processes is their potential in key technological applications, e.g., three-dimensional microfabrication,¹ optical data storage,² and bioimaging.³,⁴ Besides, 2PA can also be invaluable in spectral investigations of atoms and molecules. Indeed, as highlighted by Neusser and Schlag, the 2PA measurement setup might allow obtaining high-resolution spectra below the Doppler width.⁵ 2PA spectroscopy is also advantageous to tackle questions related to symmetry,⁶,⁷ in particular, the characterization of states that are dark with usual one-photon absorption (1PA) spectroscopies. These dark states might be strictly forbidden due to selection rules (e.g., due to the centrosymmetry of the molecule, as defined by the well-known Laporte’s rule) or simply too weak to be observed with 1PA approaches. An interesting illustration is given by the $S_0 \rightarrow S_2$ transition in Reichardt’s dye which presents a strong signal in the 2PA spectra but a very weak one in its 1PA counterpart.⁸ Because of its potential for characterizing electronic structure, 2PA spectroscopy was used during the

70s and 80s of the past century to probe the electronic structure of conjugated polyenes and other important chromophores.⁷,⁹

This large panel of applications explains why one can find a plethora of studies aiming at establishing structure–property relationships for two-photon absorbers.¹⁰⁻¹⁴ Unsurprisingly, experimental interpretations of 2PA observations have been supported by theoretical simulations since the 1970s.¹⁵⁻¹⁷ The high interpretive power of computational models remains exploited today to understand the nonlinear absorption activity.¹⁸⁻²⁰ This has been accompanied by the constant extension of the palette of methods allowing determining purely electronic two-photon transition strengths.²¹⁻²⁸ However, in order to model the band shapes in electronic 2PA spectra of molecules, one also needs to account for molecular vibrations.

**Received:** November 14, 2024

**Revised:** March 17, 2025

**Accepted:** March 18, 2025

**Published:** March 26, 2025

<figure><img src="image_9.png" /></figure>
Luo, Ågren, and collaborators have pioneered this field by performing the first ab initio calculations of the vibrationally resolved 2PA spectra.29–32 During the last two decades several groups also contributed to these efforts.33–37 These works demonstrated that computer simulations of vibrationally resolved 1PA and 2PA spectra are invaluable for interpreting the observed spectroscopic features. Indeed, absorption bands corresponding to the same electronic state might yield quite different shapes in the two types of spectra, as illustrated by the experimental data recorded for the chromophores of fluorescent proteins.38 These differences were attributed to non-Condon effects,39,40 and such effects might be potentially important in other fluorophores.41 Moreover, accurate modeling of 2PA spectra may reduce the need for complex experiments that often rely on highly sensitive detection systems and powerful laser sources.42

The primary objective of the present work is to provide a systematic evaluation of the importance of non-Condon effects in 2PA spectra across a significant range of small fluorescent probes. We propose using DFT for this analysis, as it has been observed to provide qualitatively valuable results for purely electronic 2PA with a good balance between accuracy and computational cost.20,43–50 However, the ability of DFT to capture vibronic effects has been largely overlooked. Consequently, we assess here the influence of the selected density functional approximation (DFA) on key parameters, notably, the second-order transition moments and their derivatives, comparing the data to wave function methods to appraise both their reliability and accuracy.

## 2. METHODS

**2.1. First-Principle Approaches.** The macroscopic expression of the 2PA is generally quantified as the cross section of simultaneous absorption of two isoenergetic photons, that can be expressed as51

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

where $\alpha$, $a_0$, $c$, $\omega$, $g(2\omega)$ and $\delta^{2\text{PA}}$ respectively represent the fine structure constant, the Bohr radius, the speed of light, the frequency of the incident photon, the broadening function and the rotationally averaged two-photon transition strength. In this study we assume Gaussian line-shape function (with arbitrary width) for all molecules for comparative assessment (see the Results and Discussion for details), and the given broadening values refer to half width at half-maximum (HWHM). The two-photon absorption cross section, $\sigma^{2\text{PA}}(2\omega)$, is measurable experimentally, and it is generally expressed in Göppert-Mayer units (GM) in honor of Maria Göppert-Mayer who first theorized 2PA processes,52 where 1 GM = $10^{-50}$ cm$^4 \cdot \text{s} \cdot \text{photon}^{-1}$. $\delta^{2\text{PA}}$ can be evaluated through electronic structure methods, and it can be expressed in terms of second-order transition moments ($S_{ab}$) as53

$$ \delta^{2\text{PA}} = F\delta_F + G\delta_G + H\delta_H \quad (2) $$

where $\delta_F$, $\delta_G$, and $\delta_H$ are respectively given by

$$ \delta_F = \frac{1}{30} \sum_{a,b} S_{aa} S_{bb}^* \quad (3) $$

$$ \delta_G = \frac{1}{30} \sum_{a,b} S_{ab} S_{ab}^* \quad (4) $$

$$ \delta_H = \frac{1}{30} \sum_{a,b} S_{ab} S_{ba}^* \quad (5) $$

and $F$, $G$, $H$ are polarization variables. In what follows we assume one source of linearly polarized photons ($F = G = H = 2$). We underline that $\delta_F$, $\delta_G$, and $\delta_H$ involve products of left and right second-order transition moments in the case of non-Hermitian theories.23 Note that the components of $S_{ab}$ tensor, involving sum over states, can be evaluated using purely electronic or vibronic wave functions. In this work, we analyze the effect of non-Condon effects by using the Herzberg–Teller expansion of $S_{ab}$ tensor components with respect to vibrational normal modes (Q)

$$ S_{ab}(gQ) = S_{ab}(gQ_0) + \sum_{\nu} \frac{\partial S_{ab}(gQ_0)}{\partial Q_{\nu}} Q_{\nu} + \sum_{\nu, \mu} \frac{\partial^2 S_{ab}(gQ_0)}{\partial Q_{\nu} \partial Q_{\mu}} Q_{\nu} Q_{\mu} + \dots \quad (6) $$

where $gQ_0$ is the initial state equilibrium geometry, $S_{ab}(gQ_0)$ indicates the $S_{ab}$ tensor at the $gQ_0$ geometry, whereas $\frac{\partial^2 S_{ab}(gQ_0)}{\partial Q_\nu^2}$ and $\frac{\partial^2 S_{ab}(gQ_0)}{\partial Q_\nu^2}$ indicate the first- and second-order derivatives

with respect to the normal mode ($Q_\nu$). Detailed theoretical developments can be found in refs 33–35 and 37. The linear and higher-order terms with respect to Q account for non-Condon effects (e.g., the second-order term introduces electrical anharmonicity). In this study we simulate the vibronic spectra retaining only up to the linear term in eq 6. Moreover, the performance of various electronic-structure theories is compared for second-order normal-mode derivatives of selected $S_{ab}$ tensor elements. We obtained such parameters using three electronic-structure codes. The GAUSSIAN 16 program54 was used to perform geometry optimizations and vibrational structure calculations, while the GAMESS US26,55 and VeloxChem56 programs were used to compute the elements of the second-order transition moment tensor. Normal-mode derivatives were determined numerically using an in-house code and the numerical stability was controlled with the aid of the Rutishauser–Romberg approach.57,58 To that end, we used mesh of 10 displacements per coordinate ($\pm 2^k \Delta$, $k = 0-4$; in the case of displacements in Cartesian coordinates $\Delta$ was set to 0.005 Bohr). The analysis of Romberg triangles is performed as outlined by Medved et al.58

To perform a comparative assessment of the $S_{ab}$ derivatives, we selected as references the results obtained with the second-order coupled-cluster CC2 model59 with the resolution-of-identity approximation (RI-CC2)23 as implemented in the TURBOMOLE 7.6 program.60,61 In the TD-DFT framework, the palette of tested DFAs encompasses the semilocal generalized gradient approximation (GGA) representatives (PBE62 and BLYP63,64), the global hybrid GGA (B3LYP65 and PBE066,67), the range-separated GGA (LC-BLYP68), and the global hybrid meta-GGA (MN1S69) exchange-correlation functionals. Note that the response theory calculations of meta-GGAs were performed with the gauge-variant formulation of the kinetic energy $\tau$, however, the differences with respect to gauge-invariant results are not expected to be significant given the nature of the investigated electronic transitions.70,71 For LC-BLYP we used the software default
<caption>Scheme 1. Push-Pull Molecules Used for Analyzing the Performances of DFA (Set A), and Common Fluorophores (Set B) Considered in This Study</caption>

<figure><img src="image_10.png" /><figcaption>Set B:</figcaption></figure>

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

values, i.e., μ = 0.33 (GAMESS US) and μ = 0.47 (GAUSSIAN 16), and we additionally checked other values of the range-separation parameter μ for molecules from ref 48. We tested several different values of μ as well as its optimal tuning based on Janak's theorem.⁴⁸ Based on the statistical analysis reported in Table S1 of the Supporting Information we include μ = 0.10 in our set of DFAs. To simplify the notation, we append the decimal values of μ directly in the names of the functionals in the following, i.e., LC-BLYP10, LC-BLYP33, and LC-BLYP47.

The simulations of vibronic 2PA spectra were performed at the TD-DFT level within the mechanical harmonic approximation,³⁰,⁷² using the FCClasses program.⁷³,⁷⁴ Various vibronic models are available, each differing in the representation of the potential energy surfaces (PES) of the ground state (GS) and the target excited state (ES). Among these models, the adiabatic Hessian (AH), vertical Hessian (VH), and vertical gradient (VG) approaches are the most commonly employed. The AH method is the most computationally

demanding, as it approximates the PES by considering the minima and harmonic vibrations of both the GS and ES equilibrium geometries. In contrast, the vertical models approximate the ES minimum based on the ES Hessian (VH) or gradient (VG) at the GS equilibrium geometry.⁷² In this study, we use the VG model to assess the performance of the DFAs and evaluate the importance of non-Condon effects in simulating the 2PA spectra of the first set of molecules, Set A in Scheme 1; moreover, we conduct a qualitative comparison of the three models in the 2PA simulations of the second set of molecules, i.e., Set B in Scheme 1. While these differences do not provide direct guidance on which method to select—since this decision often relies either on convergence of the different models or on direct comparisons with experimental data—they offer valuable insights into the suitability of the harmonic approximation. Indeed would the harmonic approximation be exact for the studied PESs, there would be no differences between AH and VH methods. Hence, when anharmonic
<table><caption>Table 1. Statistical Analysis of Normal Mode Derivatives of Second-Order Transition Moments (in a.u.) for Molecules 1A–6A<sup>a</sup></caption><thead><tr><th>∂S<sub>ab</sub>/∂Q<sub>ν</sub></th><th>average</th><th>median</th><th>max (ν)</th><th>average</th><th>median</th><th>max (ν)</th><th>average</th><th>median</th><th>max (ν)</th></tr></thead><tbody><tr><td></td><td colspan="3">1A</td><td colspan="3">2A</td><td colspan="3">3A</td></tr><tr><td>xx</td><td>0.114</td><td>0.023</td><td>1.130 (58)</td><td>0.161</td><td>0.075</td><td>3.003 (65)</td><td>0.163</td><td>0.026</td><td>2.450 (53)</td></tr><tr><td>yy</td><td>0.007</td><td>0.002</td><td>0.062 (60)</td><td>0.004</td><td>0.002</td><td>0.041 (65)</td><td>0.004</td><td>0.002</td><td>0.037 (53)</td></tr><tr><td>zz</td><td>0.002</td><td>0.001</td><td>0.038 (63)</td><td>0.002</td><td>0.001</td><td>0.014 (65)</td><td>0.003</td><td>0.001</td><td>0.025 (55)</td></tr><tr><td>xy</td><td>0.024</td><td>0.010</td><td>0.325 (60)</td><td>0.015</td><td>0.007</td><td>0.090 (74)</td><td>0.013</td><td>0.003</td><td>0.124 (62)</td></tr><tr><td>xz</td><td>0.008</td><td>0.000</td><td>0.070 (65)</td><td>0.006</td><td>0.004</td><td>0.032 (65)</td><td>0.007</td><td>0.002</td><td>0.043 (58)</td></tr><tr><td>yz</td><td>0.002</td><td>0.000</td><td>0.043 (64)</td><td>0.002</td><td>0.001</td><td>0.009 (65)</td><td>0.002</td><td>0.001</td><td>0.013 (34)</td></tr><tr><td></td><td colspan="3">4A</td><td colspan="3">5A</td><td colspan="3">6A</td></tr><tr><td>xx</td><td>0.136</td><td>0.006</td><td>1.764 (62)</td><td>0.200</td><td>0.070</td><td>3.193 (61)</td><td>0.217</td><td>0.003</td><td>3.489 (59)</td></tr><tr><td>yy</td><td>0.003</td><td>0.000</td><td>0.029 (54)</td><td>0.005</td><td>0.003</td><td>0.041 (61)</td><td>0.003</td><td>0.000</td><td>0.033 (59)</td></tr><tr><td>zz</td><td>0.001</td><td>0.000</td><td>0.012 (58)</td><td>0.002</td><td>0.001</td><td>0.018 (61)</td><td>0.001</td><td>0.000</td><td>0.013 (66)</td></tr><tr><td>xy</td><td>0.012</td><td>0.000</td><td>0.133 (59)</td><td>0.018</td><td>0.008</td><td>0.096 (61)</td><td>0.013</td><td>0.000</td><td>0.117 (67)</td></tr><tr><td>xz</td><td>0.003</td><td>0.000</td><td>0.025 (24)</td><td>0.006</td><td>0.004</td><td>0.035 (61)</td><td>0.002</td><td>0.000</td><td>0.024 (24)</td></tr><tr><td>yz</td><td>0.001</td><td>0.000</td><td>0.016 (37)</td><td>0.002</td><td>0.001</td><td>0.013 (61)</td><td>0.001</td><td>0.000</td><td>0.012 (39)</td></tr></tbody></table>

effects become important, the three methods can yield substantially different spectral shapes. For calculations using the VH model, all frequencies were considered real (positive) at the excited state PES.

Following the terms present in the HT expansion of the $S_{ab}$ moments, we analyze three distinct types of vibronic 2PA spectra: (i) Franck–Condon spectra (FC, including the $S_{ab}(^6Q_0)$ term only, i.e., neglecting any dependence on $Q_ν$), (ii) HT spectra (FC + HT, both $S_{ab}$ and linear terms with respect to $Q_ν$ are included), (iii) pure HT spectra (HT, only linear term with respect to $Q_ν$ is included, FC contributions are discarded). Since for the investigated states, the purely electronic transitions are allowed ($\delta^{2PA} \neq 0$), the pure HT spectra have no physical meaning. The latter simulations were performed to pinpoint the effect of the normal-mode derivatives of the second-order transition moment on the vibronic band shapes. All vibronic calculations were performed for temperature 0 K using the time-dependent (TD) formulation and rely on Cartesian coordinates. FC + HT spectra were generated using the so-called HTi approach of FCClasses (i.e., derivatives were evaluated around the initial state geometry), while pure HT contributions have been computed by setting $S_{ab}(^6Q_0) = 0$ for all tensor components. The Pople's 6–31+G(d) atomic basis set<sup>75</sup> is employed for all the computations of geometries, frequencies, and 2PA transitions, unless otherwise stated. The vibronic spectra simulations using AH, VG, and VH were performed using the same set of 2PA properties ($S$, $dS/dQ$) evaluated at the equilibrium geometry of the initial state.

**2.2. Molecular Sets.** For assessing various DFAs, we selected a series of $π$-conjugated push–pull molecules constituting Set A in Scheme 1. These molecules exhibit moderate 2PA responses, reaching up to dozens of GM at absorption band maximum. As a reference, we note that the topologically similar donor–$π$-acceptor dye 4-dimethylamino-4′-nitrostilbene shows 2PA cross-section as large as 116 GM in CHCl<sub>3</sub> solution.<sup>76</sup> It should be highlighted that all transitions have a $π \to π^*$ character. Such transitions can typically be quite accurately described by single-determinant methods, including the RI-CC2 approach used as reference herein. Indeed, it has been shown by some of us that it reliably reproduces vertical

$\delta^{2PA}$ when dealing with such states.<sup>77</sup> The B3LYP/cc-pVTZ optimized geometries of the molecules of set A can be found in ref 48.

To evaluate non-Condon effects, we investigated a series of substituted fluorophores shown in Set B in Scheme 1, exhibiting meaningful $\sigma^{FC}$ values (greater than 1 GM) for one of their lowest ESs. This set includes eight common organic dyes, two dyes carrying a BF<sub>2</sub> group that have demonstrated large experimental 2PA,<sup>78</sup> and three chromophores found in fluorescent proteins. The corresponding geometries optimized at the LC-BLYP33/6-31+G* are reported in the Supporting Information. As mentioned above, all considered excitations correspond to $π \to π^*$ transitions to the first excited state in all molecules. The computed vibrationally resolved 1PA of these molecules, along with their experimental spectra in low-polarity solvents (when available), and corresponding references can be found in Table S8 and Figure S3 in the Supporting Information. These spectra do not show significant HT effects as applying the FC + HT method does not significantly alter the 1PA spectra obtained using the FC model. Furthermore, this approximation qualitatively captures the shape of the experimental spectra as illustrated in Figure S3 of the Supporting Information.

## 3. RESULTS AND DISCUSSION

The accuracy of the computed vibronic 2PA spectra depends on the following key factors: (i) the distribution of vibronic transitions; (ii) the electronic two-photon transition strength at the equilibrium geometry and, presumably, (iii) non-Condon effects (involving geometric derivatives of second-order transition moments). All three aspects are considered below. First, we evaluate the performance of DFAs in predicting transition moments and their derivatives with respect to normal modes on molecules of set A. Subsequently, we assess the magnitude of HT effects and the impact of the choice of the vibronic model on molecules of set B.

**3.1. Benchmark of Density Functional Approximations.** The numerical evaluation of geometrical derivatives of second-order transition moment is a computationally expensive task, especially at the RI-CC2 level, explaining why we focus on a subset of normal modes. The selection is performed

<footnote><sup>a</sup>All results were obtained at the LC-BLYP47/6-31+G(d) level of theory.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 1.** First derivative of the second-order transition moment $S_{xx}$ with respect to normal modes for molecules 1A-6A calculated at the LC-BLYP47/6-31+G(d) level of theory. Red vertical lines correspond to the subset of normal modes selected for the in-depth analyses—see text for details.</figcaption></figure>

both with respect to the number of modes and the number of components of the $S_{ab}$ tensor. Our goal is to select no more than 5 vibrational normal modes and one tensor component for each molecule from set A. The selection procedure of the subset of normal modes was performed as follows. First, the complete set of first derivatives with respect to normal modes $\frac{\partial S_{ab}}{\partial Q_\nu}$ was obtained for molecules 1A-6A with the LC-BLYP47 functional. The statistical analysis in terms of average, median, and maximum value (with corresponding mode number $\nu$) for all tensor components $\frac{\partial S_{ab}}{\partial Q_\nu}$ was performed. The summary presented in Table 1 indicates a predominant contribution from the longitudinal $S_{xx}$ tensor element. Additionally, the significant difference between the median and the average values for all studied molecules suggests a nonstandard distribution. The median-mean discrepancy and the large difference between the maximum and average values indeed indicate that most of the normal modes possess low $\frac{\partial S_{ab}}{\partial Q_\nu}$ values.

This is also seen in Figure 1 where the $\frac{\partial S_{xx}}{\partial Q_\nu}$ values are plotted for each normal mode $\nu$.

The selection of the dominant longitudinal component $S_{xx}$ and its derivatives with respect to three (molecules 1-5) or four (molecule 6) vibrational modes for assessing DFAs was performed examining the data presented in Table 1 and Figure 1. In more detail, we selected the first derivatives of $S_{xx}$ with respect to the following vibrational modes: $\nu_{42}$, $\nu_{57}$, and $\nu_{58}$ (1); $\nu_{65}$, $\nu_{67}$, and $\nu_{71}$ (2); $\nu_{53}$, $\nu_{60}$, and $\nu_{61}$ (3); $\nu_{62}$, $\nu_{64}$, and $\nu_{65}$ (4); $\nu_{61}$, $\nu_{63}$, and $\nu_{67}$ (5); $\nu_{47}$, $\nu_{59}$, $\nu_{65}$, and $\nu_{77}$ (6). In the latter case, we included the values of $\frac{\partial S_{xx}}{\partial Q_\nu}$ with respect to four vibrational normal modes because the first derivatives with respect to the two in-plane scissoring bending vibrational modes ($\nu_{47}$ and $\nu_{65}$) present high and similar values. The selected vibrational modes are highlighted in red in Figure 1, while the atomic displacements corresponding to the considered vibrational modes can be found in Figure S1 and Table S2 in the Supporting Information.

The $S_{xx}$ values calculated at the equilibrium geometry of the selected molecules employing several DFAs as well as RI-CC2

with the cc-pVDZ atomic basis set can be found in Table S3 in the Supporting Information. The corresponding statistical analysis is given in Table S4 in the Supporting Information and the ordering of the DFAs follows the ranking provided by the mean absolute percentage error (MAPE). It turns out that the *meta*-GGA global hybrid MN15 is superior to other DFAs with a MAPE of 15%, while both B3LYP and PBE0 provide less satisfactory results with respective MAPEs of 31% and 25%. Interestingly, conventional RSH functionals, LC-BLYP33 and LC-BLYP47, perform worse with respective MAPEs of 41% and 48%. In contrast, LC-BLYP10 is more satisfying with a MAPE of 26%, consequently outperforming B3LYP. Finally, as expected, a poor performance is found for the two evaluated semilocal GGA functionals, PBE and BLYP, with MAPE of 52%.

The values of the $\frac{\partial S_{xx}}{\partial Q_\nu}$ obtained with different DFAs and RI-CC2 are provided in Table S5 in the Supporting Information, while the corresponding statistical analysis is presented in Table 2. Here, we used normal modes determined at LC-BLYP47/6-31+G(d,p) level of theory for computing derivatives at all other levels of theory, i.e., numerical derivatives were evaluated based on displacements along normal modes obtained from the above-mentioned method.

<table><caption>**Table 2.** Statistical Analysis of $\frac{\partial S_{xx}}{\partial Q}$ (a.u.) Performed for 19 Vibrational Modes of Molecules 1A-6A<sup>a</sup></caption><thead><tr><th>DFA</th><th>MAPE [%]</th><th>MAE [a.u.]</th><th>SDE [a.u.]</th><th>RMSE [a.u.]</th><th>MAX AE [a.u.]</th></tr></thead><tbody><tr><td>LC-BLYP33</td><td>38</td><td>0.91</td><td>1.22</td><td>1.2</td><td>3.82</td></tr><tr><td>LC-BLYP47</td><td>40</td><td>0.94</td><td>1.26</td><td>1.23</td><td>3.97</td></tr><tr><td>MN15</td><td>53</td><td>0.86</td><td>1.31</td><td>1.3</td><td>3.86</td></tr><tr><td>PBE0</td><td>74</td><td>1.36</td><td>2.01</td><td>1.97</td><td>5.71</td></tr><tr><td>B3LYP</td><td>83</td><td>1.52</td><td>2.22</td><td>2.17</td><td>6.2</td></tr><tr><td>LC-BLYP10</td><td>90</td><td>1.76</td><td>2.39</td><td>2.34</td><td>6.18</td></tr><tr><td>PBE</td><td>129</td><td>2.44</td><td>3.28</td><td>3.2</td><td>7.49</td></tr><tr><td>BLYP</td><td>131</td><td>2.47</td><td>3.3</td><td>3.23</td><td>7.51</td></tr></tbody></table>

<footnote><sup>a</sup>The errors are calculated with respect to the RI-CC2 results. The data are ordered according to increasing MAPE. The cc-pVDZ atomic basis set was used in all calculations.</footnote>
As expected the relative errors of all DFAs with respect to RI-CC2 are much larger when considering derivatives rather than $S_{xx}$. In contrast to what was observed for $S_{xx}$, MN15 underperforms with respect to LC-BLYP33 and LC-BLYP47 in all indicators but the mean absolute error (MAE). In particular, the MAPE indicates that the most adequate functionals are LC-BLYP33 (38%) and LC-BLYP47 (40%), followed by MN15 (53%). PBE0 and B3LYP are the next functionals with large MAPEs of 74% and 83%, respectively, while LC-BLYP10 delivers a MAPE of 90%. Again, GGA functionals are unsatisfactory with very large MAPEs of ca. 130%. In short, LC-BLYP33 is the “best” functional (lowest MAPE), and therefore it will be used to account for the HT effects in the simulations of the vibronic 2PA spectra of molecules from set B.

We also computed the second derivatives of $S_{xx}$ with respect to vibrational normal modes ($\frac{\partial^2 S_{xx}}{\partial Q_\nu^2}$). These derivatives are necessary to account for the electrical anharmonicity in the simulation of vibronic 2PA spectra. This analysis was done for a subset of modes where numerical stability was satisfactory, i.e., numerical errors in $\frac{\partial^2 S_{xx}}{\partial Q_\nu^2}$ were less than 4% (we underline that this rather large value corresponds to one outlier, the numerical errors are less than 0.1% for the majority of cases). In more detail, we report the data related to 6 modes of molecules **1A** and **4A** in Table S6 in the Supporting Information, and the relative statistical analysis in Table 3.

<table><caption>**Table 3. Statistical Analysis of $\frac{\partial^2 S_{xx}}{\partial Q^2}$ (a.u.) Performed for 6 Vibrational Modes of Molecule 1A and 4A<sup>a</sup>**</caption><thead><tr><th>DFA</th><th>MAPE [%]</th><th>MAE [a.u.]</th><th>SDE [a.u.]</th><th>RMSE [a.u.]</th><th>MAX AE [a.u.]</th></tr></thead><tbody><tr><td>LC-BLYP47</td><td>38</td><td>0.013</td><td>0.017</td><td>0.017</td><td>0.034</td></tr><tr><td>LC-BLYP33</td><td>48</td><td>0.015</td><td>0.020</td><td>0.019</td><td>0.043</td></tr><tr><td>MN15</td><td>102</td><td>0.029</td><td>0.040</td><td>0.037</td><td>0.062</td></tr><tr><td>PBE0</td><td>153</td><td>0.047</td><td>0.067</td><td>0.063</td><td>0.112</td></tr><tr><td>LC-BLYP10</td><td>155</td><td>0.046</td><td>0.064</td><td>0.059</td><td>0.091</td></tr><tr><td>B3LYP</td><td>171</td><td>0.053</td><td>0.077</td><td>0.073</td><td>0.138</td></tr><tr><td>PBE</td><td>215</td><td>0.060</td><td>0.081</td><td>0.076</td><td>0.118</td></tr><tr><td>BLYP</td><td>219</td><td>0.062</td><td>0.083</td><td>0.078</td><td>0.128</td></tr></tbody></table>

<sup>a</sup>The errors are calculated with respect to RI-CC2 results. The data are ordered according to increasing MAPE. The cc-pVDZ atomic basis set was employed in all calculations.

The data given in Table 3 indicate that the performance of the two standard LC-BLYP variants is very good in the case of second derivatives. Indeed the errors (MAE, SDE, RMSE, and MAPE) obtained with LC-BLYP47 and LC-BLYP33 are almost half of their MN15 counterparts, which ranks third in the list. In more detail, the first two functionals have a MAPE value of 38% and 48% respectively, while the MN15 MAPE is 102%. The other DFAs deliver unsatisfactory estimations of the second derivatives with MAPE larger than 100%. LC-BLYP10, PBE0, and B3LYP deliver similar inaccuracies (MAPE between 155% and 171%) whereas PBE (215%) and BLYP (219%) are particularly poor. A summary of MAPE for all properties examined in this Section can be found in Figure 2.

Let us now assess the impact of $\frac{\partial S_{xx}}{\partial Q_\nu}$ derivatives, as given by different DFAs, on the vibronic 2PA spectra of molecules **1A-6A**. A comparison of the band shapes simulated using $\frac{\partial S_{xx}}{\partial Q_\nu}$ for the selected normal modes obtained from the most accurate LC-BLYP33, least accurate BLYP, and reference RI-CC2 levels of theory is presented in Figure 3 for molecule **6A** (see also Figure S2 in the Supporting Information). The rationale behind the inclusion of the BLYP functional in this part of the study is to demonstrate how significant statistical errors determined for individual components affect the simulated 2PA band shapes. We recall that we determined all the $\frac{\partial S_{ab}}{\partial Q_\nu}$ values using LC-BLYP47/6-31+G(d) only, while these derivatives for selected subset of vibrational normal modes were also computed at the RI-CC2, LC-BLYP33, and BLYP levels. Hence, to illustrate how the choice of DFA affects the 2PA band topologies, we use the $S_{ab}$ and $\frac{\partial S_{ab}}{\partial Q_\nu}$ values determined at the LC-BLYP47/6-31+G(d) level, replacing the derivatives for the selected modes with those computed at the RI-CC2, LC-BLYP33, and BLYP levels. In other words, for a given DFA (or RI-CC2) the Franck-Condon contribution and derivatives for 3N-9 (3N-10 for **6A**) modes are determined at LC-BLYP47/6-31+G(d) level, while the remaining derivatives for 3 (4 for **6A**) key modes are determined at DFA(RI-CC2)/cc-pVDZ level.

Hence, it comes as no surprise that Figure 3 and Figure S2 in the Supporting Information show similar peak positions for all three methods, i.e., all schemes restore the FC contribution as predicted by LC-BLYP47/6-31+G(d) which are predominant for these systems. However, analyzing pure HT vibronic

<figure><img src="image_13.png" /><figcaption>**Figure 2.** MAPE of $S_{xx}$, $\frac{\partial S_{xx}}{\partial Q_\nu}$, and $\frac{\partial^2 S_{xx}}{\partial Q_\nu^2}$ for molecules **1A-6A** using the RI-CC2 results as benchmarks. The first-order derivatives are performed on the selected 24 vibrational modes of molecules **1A-6A**, while the second-order derivatives on the 6 modes of molecules **1A** and **4A** (see the text for the list of modes included). The cc-pVDZ atomic basis set is employed.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 3. Vibronic FC + HT (left) and pure HT (right) 2PA spectra for molecule 6A.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 4. Percentage contributions of HT effects (% HT) on the total intensity of the studied $\pi \rightarrow \pi^*$ 2PA electronic transitions for 1B–13B. See eqs 7 and 8.</figcaption></figure>

spectra (with no FC contributions) one notices that the Herzberg–Teller spectra for 6A reveal that the relative peak intensities might strongly differ when using RI-CC2 (or LC-BLYP33) or BLYP. These substantial changes are due to the differences in the magnitudes of $\partial S_{ab}/\partial Q_{\nu}$ derivatives for four key normal modes only. It is clear that HT contributions, which affect the relative intensities of the peaks, crucially depend on the level of theory. We stress that LC-BLYP33 predicts a very similar intensity pattern as RI-CC2 but the intensities of individual peaks are much smaller. In contrast, the shape of the pure HT spectra predicted by BLYP strongly differs from its RI-CC2 counterpart. This outcome can be easily explained using the data of Table S5 of the Supporting Information, i.e., one finds that derivative values computed using LC-BLYP33 are roughly half their RI-CC2 counterparts. Finally, we underline that some DFAs deliver wrong derivative signs compared to RI-CC2 (see Table S5 of the Supporting Information).

### 3.2. Importance of HT Contributions in Common

Dyes. Based on the above results, the LC-BLYP33 functional was used to compute the geometries and vibronic 2PA spectra of the 1B–13B molecules. To estimate the impact of HT effects on the total intensity of a 2PA electronic band, $\delta^{2PA}$, we can separate the FC and HT effects by using the expression for total intensity at 0 K$^{30}$

$$ \begin{aligned} \delta^{2\text{PA}} &= 2 \sum_{a,b} (S_{aa}S_{bb}^* + 2S_{ab}S_{ba}^*) + 2 \sum_{\nu,a,b} \frac{\hbar}{2\omega_{\nu}} \\ & \qquad \left( \frac{\partial S_{aa}}{\partial Q_{\nu}} \frac{\partial S_{bb}}{\partial Q_{\nu}} + 2 \frac{\partial S_{ab}}{\partial Q_{\nu}} \frac{\partial S_{ba}}{\partial Q_{\nu}} \right) \\ &= \delta^{\text{FC}} + \delta^{\text{HT}} \end{aligned} \qquad (7) $$

where $\delta^{\text{FC}}$ includes the FC contributions to the total 2PA intensity, while $\delta^{\text{HT}}$ provides the HT contribution. In Figure 4 we report the percentage contribution of HT effects, namely

$$ \% \text{HT} = \frac{\delta^{\text{HT}}}{\delta^{\text{FC}} + \delta^{\text{HT}}} \cdot 100 \qquad (8) $$

In Figure 4, we observe that HT contributions to the 2PA total intensity have a minimal impact for molecules 1B–13B. The HT term contributes minimally to the total $\delta$ value, with contributions ranging from 1.4% to 6.4%. Interestingly, molecules with the largest percentage of HT contributions tend to have low $\delta^{\text{FC}}$ values, stressing the fact that HT transitions are not predominant for compounds of practical interest. To further explore the effect of HT contributions, Figure 5 presents the vibronic 2PA spectra simulated using both the FC and FC + HT approximations, with the VG model and a homogeneous Gaussian broadening of 0.005 eV.

Figure 5 shows that HT contributions not only have a minor influence on the total intensity of the 2PA $\pi \rightarrow \pi^*$ electronic bands, but also that they do not play a significant role in determining the vibronic shape of the bands. The inclusion of HT couplings results in slight variations in the intensity of a few peaks compared to the FC approximation. For instance, in
<figure><img src="image_16.png" /><figcaption>**Figure 5.** Vibrational LC-BLYP33/6-31+G(d) $\sigma^{2PA}$ spectra obtained with the TD method, the PES VG model, and a broadening of 0.005 eV, for the different fluorescent dyes computed including FC couplings only (in blue) or both FC and HT terms (in red).</figcaption></figure>

molecule **3B**, we observe a significant reduction in 2PA intensity for high-wavelength modes and an increase for lower-

wavelength modes. Conversely, molecule **13B** shows a slight intensity increase in high-wavelength modes. To quantify the
<figure><img src="image_17.png" /><figcaption>**Figure 6.** Vibrational LC-BLYP33/6-31+G(d) $\sigma^{2PA}$ spectra employing the TD method, and the VG, VH, and AH PES models. A Gaussian broadening with a width of 0.100 eV is applied and both FC and HT couplings are included.</figcaption></figure>

impact of HT contributions, we computed vibrationally resolved spectra using a broader line width of 0.100 eV, which is better suited for direct comparisons with experimental data (see Figure S6 in the Supporting Information). We interpolated the peaks from these spectra and compared the

results obtained using the FC and FC + HT models. The corresponding data are reported in Table S9 in the Supporting Information. Focusing on the peak positions, we see that the slight increase in low-wavelength peaks identified in Figure 5 leads to small blueshifts in the spectra, an effect being
negligible in most cases except for molecule **3B**, for which HT couplings cause a blueshift of 8.7 nm (0.02 eV). Conversely, a few molecules (**1B**, **9B**, **12B**) exhibit slight (less than 2 nm) redshifts when HT effects are included. Including HT effects minimally impacts the intensity of the maximum peaks with variations of less than 1 GM, except for molecules **9B** and **10B**, which show slight enhancement of the 2PA signal by 3.2 GM, 1.3 GM, and 1.1 GM, respectively.

**3.2.1. Impact of the Vibronic Model.** For molecules **1B–13B** we compared the results obtained with the VG, VH, and AH PES models. To estimate the suitability of the harmonic approximation in the AH framework, we computed the RMSD values between the ground- and excited-state geometries (see Table S11 in the Supporting Information). For the studied molecules the RMSD values remain below 0.17 Å, indicating that the AH approximation is likely suitable for those cases. Figure 6 provides the corresponding vibronic spectra with a 0.100 eV broadening. As above, we interpolated the values of main peaks and identified the corresponding wavelengths, as well as the full-width half maxima ($\Delta\nu$) and the integral surfaces (I). These values are listed in Table S10 in the Supporting Information.

By and large, the three methods yield qualitatively similar spectra with relatively minor topological differences. Transition wavelength shifts are observed when moving from the simpler VG model to the more sophisticated VH and AH methods. Specifically, for the VH model a redshift is noted for all the molecules with respect to VG values, ranging from 5.8 to 30.2 nm (0.01 to 0.07 eV) for VH model. Similarly, for AH model we observe a redshift for 11 molecules ranging from 3.8 (0.02 eV) to 26.7 nm (0.05 eV). Globally such a trend is expected.⁷⁹⁻⁸¹ Conversely, for molecules **5B** and **12B**, AH exhibits a modest blueshift of 1.8 nm (0.01 eV) for **12B** and 7.1 nm (0.02 eV) for **5B**. Dye **10B** shows a notable variation in the maximum $\sigma^{2\text{PA}}$ for AH yielding intensities smaller by 10.3 GM, while for VH it is smaller of 3.7 GM. However, these values are less significant in percentage given that the $\sigma^{2\text{PA}}$ for this molecule is ca. 100 GM. Moreover, in molecule **5B** the AH spectra deviates significantly from the other two models, showing a smaller $\sigma^{2\text{PA}}$ of 3.8 GM, compared to 7.0 GM (VG) and 6.7 GM (VH). Additionally, AH produces a broader spectral profile (see below). For the remaining molecules, variations among the three methods are negligible, with VG consistently yielding the largest responses.

Focusing now on the full-width at half-maximum (fwhm) values, we notice that VH and AH tend to produce larger values than VG, a trend previously observed in 1PA spectra.⁸² Although these differences are rather negligible for most molecules, they become significant for molecule **5B**, with fwhm differences of 76.1 nm for AH. This aligns with the intensity changes discussed above and correlates with the larger RMSD observed between ground and excited-state structures, indicating a more pronounced geometry relaxation in the excited state. A few exceptions exist, where VG yields slightly larger fwhm values, specifically for molecules **9B** (relative to VH only) and **11B** (relative to AH), yet these differences remain minimal (<1.4 nm).

## 4. CONCLUSIONS

In this work, we systematically assessed the importance of non-Condon effects in simulating 2PA spectra. Using a series of representative molecules, we evaluated the performance of various DFAs against RI-CC2 in reproducing the tensor

elements of two-photon transition moments and their first- and second-order derivatives with respect to normal modes. The present findings align with previous studies regarding electronic contributions to two-photon transition strengths,⁴⁸,⁴⁹ i.e., the global hybrid *meta*-GGA MN15 functional outperforms all the other studied range-separated hybrid, hybrid, and GGA functionals for the investigated molecules. The analysis of the first- and second-order derivatives of second-order transition moment with respect to the normal modes indicate an unsatisfactory accuracy across all DFAs, except for LC-BLYP33 and LC-BLYP47, which deliver an error margin in the range of 40–50% compared to the RI-CC2 reference. In particular, LC-BLYP33 can likely be recommended for its ability to accurately predict the peak positions due to the HT contributions to the 2PA spectrum, though it tends to underestimate the HT contribution to the total intensity of the 2PA $\pi \rightarrow \pi^*$ electronic bands. We emphasize that while DFAs for 2PA still fall short of the accuracy achieved by wave function-based methods, they can provide qualitatively reliable results.

Building on these findings, we explored the influence of non-Condon effects on a series of commonly used dyes and fluorescent chromophores within proteins. Across different models of dipole expansion, HT contributions were found to be secondary to their FC counterparts: the largest observed peak shifts are around 9 nm while the intensity changes are trifling. Given the overall insignificant effect of linear Herzberg–Teller term on the shape of the spectra of the investigated fluorophores, one may assume the electrical anharmonicity effects to be negligible. Additionally, we evaluated the AH, VH, and VG vibronic models for simulating the 2PA vibronic spectra, finding that these models produced consistent results for all cases, indicating that the harmonic approximation is generally suitable for such molecular systems. The next logical steps would be to extend this analysis to larger systems and benchmark normal-mode derivatives of second-order transition moments at higher levels of theory, such as CCSD or CC3, which is unfortunately beyond reach at present. To achieve this, further research into the development of analytical gradients for two-photon moments is essential, as it would greatly enhance the efficiency of 2PA spectra calculations, allowing for efficient computations on more complex molecular systems.

### ■ ASSOCIATED CONTENT

**Supporting Information**

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

Additional info: benchmark of LC-BLYP variants for two-photon absorption strengths; selection of vibrational modes; performances of density functional approximations in the simulations of two-photon transition strengths/spectra; theoretical and experimental one-photon spectra; comparison of two-photon absorption spectra for different vibronic models; Cartesian coordinates of molecules (PDF)

### ■ AUTHOR INFORMATION

**Corresponding Authors**

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

Email: Robert.Zalesny@pwr.edu.pl

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

Josep M. Luis – Institute of Computational Chemistry and Catalysis and Department of Chemistry, University of Girona, 17003 Girona, Catalonia, Spain; orcid.org/0000-0002-2880-8680; Email: Josepm.Luis@udg.edu

## Authors

Rudraditya Sarkar – Institute of Computational Chemistry and Catalysis and Department of Chemistry, University of Girona, 17003 Girona, Catalonia, Spain; Department of Chemistry, School of Science, Gandhi Institute of Technology and Management (GITAM), Hyderabad 502329, India; orcid.org/0000-0001-7915-8625

Carmelo Naim – Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes, France

Karan Ahmadzadeh – Hefei National Research Center for Physical Sciences at the Microscale, University of Science and Technology of China, Hefei, Anhui 230026, China

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

## Author Contributions

†R.S. and C.N. contributed equally to this work.

## Notes

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

The financial support from the National Science Centre (Poland) is acknowledged (grant no. 2018/30/E/ST4/00457). C.N. and D.J. acknowledge the French Agence Nationale de la Recherche (ANR) for financial support in the framework of contract no. ANR-21-CE07-0058-2 (CONDOR). J.M.L. thank the Spanish Ministry of Science (PID2022-140666NB-C22) and Generalitat de Catalunya (2021SGR00623) for financial support. The authors also thanks the Wroclaw Center for Networking and Supercomputing (Poland), as well as the GLiCID Computing Facility (Ligerien Group for Intensive Distributed Computing, 10.60487/glicid, Pays de la Loire, France) for the computational resources. R.S. thanks the European Council and the Spanish Ministry of Universities for a “María Zambrano” grant (REQ2021_C 29).

## REFERENCES

(1) Cumpston, B. H.; Ananthavel, S. P.; Barlow, S.; Dyer, D. L.; Ehrlich, J. E.; Erskine, L. L.; Heikal, A. A.; Kuebler, S. M.; Lee, I.-Y. S.; McCord-Maughon, D.; et al. Two-Photon Polymerization Initiators for Three-Dimensional Optical Data Storage and Microfabrication. *Nature* **1999**, *398*, 51–54.

(2) Parthenopoulos, D. A.; Rentzepis, P. M. Three-Dimensional Optical Storage Memory. *Science* **1989**, *245*, 843–845.

(3) Hoover, E. E.; Squier, J. A. Advances in Multiphoton Microscopy Technology. *Nat. Photonics* **2013**, *7*, 93–101.

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

(5) Neusser, H. J.; Schlag, E. W. High Resolution Spectroscopy Below the Doppler Width. *Angew. Chem., Int. Ed.* **1992**, *31*, 263–273.

(6) McClain, W.; Harris, R. *Excited States*; Academic Press: New York, 1977; Vol. 3.

(7) McClain, W. M. Two-Photon Molecular Spectroscopy. *Acc. Chem. Res.* **1974**, *7*, 129–135.

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

(9) Birge, R. R.; Bennett, J. A.; Pierce, B. M.; Thomas, T. M. Two-Photon Spectroscopy of the Visual Chromophores. Evidence for a Lowest Excited $^1\text{ag}$-like $\pi\pi^*$ State in All-trans-retinol (vitamin A). *J. Am. Chem. Soc.* **1978**, *100*, 1533–1539.

(10) Albota, M.; Beljonne, D.; Brédas, J.-L.; Ehrlich, J.; Fu, J.-Y.; Heikal, A.; Hess, S.; Kogej, T.; Levin, M.; Marder, S.; et al. Design of Organic Molecules with Large Two-Photon Absorption Cross Sections. *Science* **1998**, *281*, 1653–1656.

(11) Lee, W.-H.; Lee, H.; Kim, J.-A.; Choi, J.-H.; Cho, M.; Jeon, S.-J.; Cho, B. Two-Photon Absorption and Nonlinear Optical Properties of Octupolar Molecules. *J. Am. Chem. Soc.* **2001**, *123*, 10658–10667.

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

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

(14) Pawlicki, M.; Collins, H. A.; Denning, R. G.; Anderson, H. L. Two-Photon Absorption and the Design of Two-Photon Dyes. *Angew. Chem., Int. Ed.* **2009**, *48*, 3244–3266.

(15) Birge, R.; Pierce, B. A Theoretical Analysis of the Two-Photon Properties of Linear Polyenes and the Visual Chromophores. *J. Chem. Phys.* **1979**, *70*, 165–178.

(16) Marchese, F.; Seliskar, C.; Jaffé, H. H. The Use of CNDO in Spectroscopy. XV. Two-Photon Absorption. *J. Chem. Phys.* **1980**, *72*, 4194–4203.

(17) Dick, B.; Hohlneicher, G. Importance of Initial and Final States as Intermediate States in Two-Photon Spectroscopy of Polar Molecules. *J. Chem. Phys.* **1982**, *76*, 5755–5760.

(18) Alam, M. M.; Chattopadhyaya, M.; Chakrabarti, S.; Ruud, K. Chemical Control of Channel Interference in Two-Photon Absorption Processes. *Acc. Chem. Res.* **2014**, *47*, 1604–1612.

(19) Alam, M. M.; Beerepoot, M. T. P.; Ruud, K. Channel Interference in Multiphoton Absorption. *J. Chem. Phys.* **2017**, *146*, 244116.

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

(21) Olsen, J.; Jørgensen, P. Linear and Nonlinear Response Functions for an Exact State and for and MCSCF State. *J. Chem. Phys.* **1985**, *82*, 3235–3264.

(22) Hättig, C.; Christiansen, O.; Jo/rgensen, P. Multiphoton Transition Moments and Absorption Cross Sections in Coupled Cluster Response Theory Employing Variational Transition Moment Functionals. *J. Chem. Phys.* **1998**, *108*, 8331–8354.

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

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

(25) Salék, 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.

(26) Zahariev, F.; Gordon, M. S. Nonlinear Response Time-Dependent Density Functional Theory Combined with the Effective Fragment Potential Method. *J. Chem. Phys.* **2014**, *140*, 18A523.
(27) Parker, S. M.; Rappoport, D.; Furche, F. Quadratic Response Properties From TDDFT: Trials and Tribulations. *J. Chem. Theory Comput.* **2018**, *14*, 807–819.

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

(29) Luo, Y.; Ågren, H.; Knuts, S.; Minaev, B.; Jørgensen, P. Response Theory Calculations of the Vibrationally Induced $ ^{1}A_{ig} - ^{1}B_{Iu} $ Two-Photon Spectrum of Benzene. *Chem. Phys. Lett.* **1993**, *209*, 513–518.

(30) Macak, P.; Luo, Y.; Ågren, H. Simulations of Vibronic Profiles in Two-Photon Absorption. *Chem. Phys. Lett.* **2000**, *330*, 447–456.

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

(32) Minaev, B.; Wang, Y.-H.; Wang, C.-K.; Luo, Y.; Ågren, H. Density Functional Theory Study of Vibronic Structure of the First Absorption Qx Band in Free-Base Porphin. *Spectrochim. Acta Part A* **2006**, *65*, 308–323.

(33) Bishop, D. M.; Luis, J. M.; Kirtman, B. Vibration and Two-Photon Absorption. *J. Chem. Phys.* **2002**, *116*, 9729–9739.

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

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

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

(37) Lin, N.; Santoro, F.; Rizzo, A.; Luo, Y.; Zhao, X.; Barone, V. Theory for Vibrationally Resolved Two-Photon Circular Dichroism Spectra. Application to (R)-(+)-3-methylcyclopentanone. *J. Phys. Chem. A* **2009**, *113*, 4198–4207.

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

(39) Drobizhev, M.; Makarov, N. S.; Tillo, S. E.; Hughes, T. E.; Rebane, A. Describing Two-Photon Absorptivity of Fluorescent Proteins with a New Vibronic Coupling Mechanism. *J. Phys. Chem. B* **2012**, *116*, 1736–1744.

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

(41) Zaleśny, R.; Murugan, N. A.; Tian, G.; Medved', M.; Ågren, H. First-Principles Simulations of One- and Two-Photon Absorption Band Shapes of the Bis(BF₂) Core Complex. *J. Phys. Chem. B* **2016**, *120*, 2323–2332.

(42) de Reguardati, S.; Pahapill, J.; Mikhailov, A.; Stepanenko, Y.; Rebane, A. High-Accuracy Reference Standards for Two-Photon Absorption in the 680–1050 nm Wavelength Range. *Opt. Express* **2016**, *24*, 9053–9066.

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

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

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

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

(47) Knysz, 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*, No. e202200137.

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

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

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

(51) Cronstrand, P.; Luo, Y.; Ågren, H. In *Response Theory and Molecular Properties (A Tribute to Jan Lindberg and Poul Jørgensen)*; Jensen, H., Ed.; *Adv. Quantum Chem.*; Academic Press, 2005; Vol. 50, pp 1–21.

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

(53) Monson, P. R.; McClain, W. M. Polarization Dependence of the Two-Photon Absorption of Tumbling Molecules with Application to Liquid 1-Chloronaphthalene and Benzene. *J. Chem. Phys.* **1970**, *53*, 29–37.

(54) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; et al. *Gaussian 16*. Revision C.01; Gaussian Inc: Wallingford CT, 2016.

(55) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. *J. Comput. Chem.* **1993**, *14*, 1347–1363.

(56) 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* **2020**, *10*, No. e1457.

(57) Rutishauser, H. Ausdehnung Des Rombergschen Prinzips. *Num. Math* **1963**, *5*, 48–54.

(58) Medved, M.; Stachová, M.; Jacquemin, D.; André, J.-M.; Perpète, E. A. A Generalized Romberg Differentiation Procedure for Calculation of Hyperpolarizabilities. *J. Mol. Struct.* **2007**, *847*, 39–46.

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

(60) TURBOMOLE V7.6 2020; a Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, Since, 2007, Available From https://www.turbomole.org.

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

(62) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for Mixing Exact Exchange with Density Functional Approximations. *J. Chem. Phys.* **1996**, *105*, 9982–9985.

(63) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. *Phys. Rev. A* **1988**, *38*, 3098–3100.

(64) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula Into a Functional of the Electron Density. *Phys. Rev. B* **1988**, *37*, 785–789.

(65) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. *J. Chem. Phys.* **1993**, *98*, 5648–5652.
(66) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0Model. *J. Chem. Phys.* **1999**, *110*, 6158–6170.

(67) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew-Burke-Ernzerhof Exchange-Correlation Functional. *J. Chem. Phys.* **1999**, *110*, 5029–5036.

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

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

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

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

(72) Avila Ferrer, F. J.; Santoro, F. Comparison of Vertical and Adiabatic Harmonic Approaches for the Calculation of the Vibrational Structure of Electronic Spectra. *Phys. Chem. Chem. Phys.* **2012**, *14*, 13549–13563.

(73) Cerezo, J.; Santoro, F. Fcclasses3, Vibrationally-Resolved Spectra Simulated at the Edge of the Harmonic Approximation. *J. Comput. Chem.* **2023**, *44*, 626–643.

(74) Santoro, F.; Cerezo, J. Fcclasses3, A Code for Vibronic Calculations. V3.0.3–0-gb28e48a, 2023.

(75) Rassolov, V. A.; Ratner, M. A.; Pople, J. A.; Redfern, P. C.; Curtiss, L. A. 6–31G* Basis Set for Third-Row Atoms. *J. Comput. Chem.* **2001**, *22*, 976–984.

(76) Wielgus, M.; Bartkowiak, W.; Samoc, M.; Two-Photon Solvatochromism, I. 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.

(77) Naim, C.; Zaleśny, R.; Jacquemin, D. Two-Photon Absorption Strengths of Small Molecules: Reference CC3 Values and Benchmarks. *J. Chem. Theory Comput.* **2024**, *20*, 9093–9106.

(78) Bednarska, J.; Zaleśny, R.; Wielgus, M.; Jędrzejewska, B.; Puttreddy, R.; Rissanen, K.; Bartkowiak, W.; Ágren, H.; Ośmiałowski, B. Two-Photon Absorption of BF2-Carrying Compounds: Insights From Theory and Experiment. *Phys. Chem. Chem. Phys.* **2017**, *19*, 5705–5708.

(79) Humeniuk, A.; Bužančić, M.; Hoche, J.; Cerezo, J.; Mitrić, R.; Santoro, F.; Bonačić-Koutecký, V. Predicting Fluorescence Quantum Yields for Molecules in Solution: A Critical Assessment of the Harmonic Approximation and the Choice of the Lineshape Function. *J. Chem. Phys.* **2020**, *152*, 054107.

(80) Bousquet, M. H. E.; Papineau, T. V.; Veys, K.; Escudero, D.; Jacquemin, D. Extensive Analysis of the Parameters Influencing Radiative Rates Obtained through Vibronic Calculations. *J. Chem. Theory Comput.* **2023**, *19*, 5525–5547.

(81) do Casal, M. T.; Veys, K.; Bousquet, M. H. E.; Escudero, D.; Jacquemin, D. First-Principles Calculations of Excited-State Decay Rate Constants in Organic Fluorophores. *J. Phys. Chem. A* **2023**, *127*, 10033–10053.

(82) Petrusevich, E. F.; Bousquet, M. H. E.; Ośmiałowski, B.; Jacquemin, D.; Luis, J. M.; Zaleśny, R. Cost-Effective Simulations of Vibrationally-Resolved Absorption Spectra of Fluorophores with Machine-Learning-Based Inhomogeneous Broadening. *J. Chem. Theory Comput.* **2023**, *19*, 2304–2315.