A route to assess non-linear light–matter interactions from the increasingly popular GW-Bethe–Salpeter equation (GW-BSE) method is outlined. In the present work, the necessary analytic expressions within the static-screened exchange approximation of the BSE are derived. This enables a straightforward implementation of the computation of the first hyperpolarizability as well as two-photon absorption processes for molecular systems. Benchmark calculations on small molecular systems reveal that the GW-BSE method is intriguingly accurate for predicting both first hyperpolarizabilities and two-photon absorption strengths. Using state-of-the-art Kohn–Sham references as a starting point, the accuracy of the GW-BSE method rivals that of the coupled-cluster singles-and-doubles method, outperforming both second-order coupled-cluster and time-dependent density-functional theory.

Determining non-linear optical properties is an intriguing task, allowing one to access complex light–matter interactions. Common examples of this kind of interaction are two-photon absorption (2PA) processes or electric-field-induced second-harmonic generation (EHSG).1,2 Within wavefunction theory and later also within the framework of time-dependent density-functional theory (TD-DFT), the corresponding adaptions to assess non-linear optical properties have therefore been targeted by several groups. Hyperpolarizabilities and 2PA cross sections are therefore accessible for various wavefunction-based methods, including configuration-interaction-based methods as well as several coupled-cluster variants, such as the second-order method CC2, the third-order method CC3, and the equation-of-motion coupled-cluster singles-and-doubles method EOM-CCSD.3–9 Also within the framework of TD-DFT, assessing non-linear optical properties has become a common task in the past few years.10–13 Much less is known about the behavior of the increasingly popular Bethe–Salpeter equation (BSE) method in the non-linear regime. While the latter has emerged as a useful tool in molecular quantum chemistry in the past decade, it still lacks an implementation being capable of assessing non-linear optical properties. The great success of the BSE method to calculate excited states, therefore, encouraged us to outline an implementation of these properties. Compared to TD-DFT, the intrinsic advantage of being able to correctly predict local excitations and charge-transfer excitations alike is an intriguing feature of the BSE.14–19 Yet, using modern algorithmic techniques based on the resolution of the identity (RI) approximation,20–22 the BSE exhibits the same formal computational scaling as TD-DFT.23,24 Combining with a suitable O(N4) scaling GW method, as, for example, frequency-sampling contour deformation GW (fsCD-GW), allows for a broad application of GW-BSE to sizable molecular and periodic systems.25,26 The recent success of GW-BSE describing excited-state dipole moments is also encouraging to expand its possibilities.27,28 Wavefunction-based methods exhibit a significantly higher cost, starting from O(N5) for CC2 and quickly reaching O(N7) for CC3. This limits the effective range of molecular systems accessible by these methods. It was also revealed that comparisons between different wavefunction-based methods are generally less conclusive for non-linear response properties than they are for standard linear optical properties. For example, CC2, which is considered to perform sufficiently well for many linear response problems,29 deviates significantly from its EOM-CCSD parent method for the calculations of 2PA transition strengths or hyperpolarizabilities.4,5,30,31 The latter effect is suspected to be caused by a rather inaccurate description of double excitations within CC2 when compared to CC3,7,31 though it is not fully clear why CCSD performs comparably well then. For TD-DFT, also significant spreads are found depending on the chosen density-functional approximation (DFA). This spread is considerably more pronounced for non-linear properties than it is for standard linear properties.32,33 It would therefore certainly be valuable to have a method at hand that allows for an accurate assessment of sizable molecular systems on a feasible timescale as, for example, needed for recent developments targeting multiscale modeling of non-linear optical effects.34 The present work accordingly introduces a convenient way to assess non-linear light–matter interactions using the GW-BSE method. Subsequently, it will be shown that the GW-BSE method allows for an accurate assessment of non-linear light–matter interactions, rivaling even wavefunction-based methods.

Starting from a time-dependent field that is oscillating with the general complex frequency ω′, the linear response function can be written as13,35,36
vζ;vη(ω)=tr(v̂ζγη(ω)),
(1)
where γη is the first-order reduced density matrix,
γη(x,x)=aiXiaηϕa(x)ϕi(x)+Yiaηϕi(x)ϕa(x),
(2)
which is also referred to as the transition-density matrix. Here, we adopt the generally used notation, with the subscript indices i, j, … indicating occupied orbitals, a, b, … indicating virtual orbitals, and p, q, … indicating general orbitals. The superscripts ζ, η, θ… denote the external perturbation. Within the static screened Bethe–Salpeter equation,37 the polarization vector X,Y can be obtained by solving the corresponding symplectic Bethe–Salpeter equation (BSE).21,22,37,38 v̂ζ denotes the corresponding one-electron property operator of the external perturbation ζ and x={r,σ} denotes a spin-space coordinate,
(AB)(A+B)ωζ21Rζ=Uζ,
(3)
(A+B)(AB)ωζ21Lζ=Vζ,
(4)
imposing the condition
LTR=1.
(5)
In Eq. (3), we have used the notation Riaζ=(Xiaζ+Yiaζ), Liaζ=(XiaζYiaζ), Uiaζ=(Piaζ+Qiaζ), and Viaζ=(PiaζQiaζ). Within the static-screened BSE, the matrices (A + B) and (A − B) are defined as
(A+B)ia,jb=(εaεi)δabδij+Hai,bj+,BSE,
(6)
(AB)ia,jb=Hai,bj,BSE.
(7)
Here, we have introduced the two-electron BSE-kernel contribution,
Hpq,rs+,BSE=vpq,rsWps,qr(ω=0)Wpr,qs(ω=0),
(8)
Hpq,rs,BSE=Wps,qr(ω=0)Wpr,qs(ω=0),
(9)
with vpq,rs being a Coulomb integral and W(ω) being denoted as screened exchange. The latter is obtained from v and the inverse of the dielectric function κ,
Wpq,rs(ω)=tuκpq,tu1(ω)vtu,rs.
(10)
Details of the evaluation of W and κ within the RI approximation are outlined in Refs. 22 and 39. Pζ in Eq. (3) collects the electric,
Piae=Qaie=ϕi|r|ϕa,
(11)
or magnetic,
Piam=Qaim=12cϕi|(rR)×|ϕa,
(12)
dipole integrals, where c denotes the speed of light. For higher multipole moments, the according quadrupole, octupole, etc., integrals must be used. From the solutions of Eqs. (3) and (4) for a given frequency, the corresponding property can then be obtained as
αζη=Pζ,Qζ|Xη,Yη.
(13)
The notation {X, Y} emphasizes that the corresponding supervectors are taken to evaluate the inner product. This has been well established for the BSE, and the corresponding polarizabilities from the BSE were found to be in very good agreement with experimental data.40,41
Quadratic response properties from the BSE can be derived in a similar manner as the related time-dependent DFT quantities. In second order, the quadratic response function can be written as11,
vζ;vη(ω),vθ(ω)=tr(v̂ζγηθ(ω,ω)),
(14)
where γηθ is now the second-order single-particle reduced density matrix. As noted in Ref. 13, γηθ is generally given as
γζη=KζηXζηYζηKζη
(15)
And, therefore, has contributions not only to the occupied–virtual part but also to the occupied–occupied and virtual–virtual blocks. The corresponding symmetrized subblocks of Eq. (15) can be obtained as
Kijζη,+=14aRiaζRjaηLiaζLjaη+(ζη),
(16)
Kijζη,=14aLiaζRjaηRiaζLjaη+(ζη),
(17)
Kabζη,+=14iRiaζRibηLiaζLibη+(ζη),
(18)
Kabζη,=14iLiaζRibηRiaζLibη+(ζη),
(19)
with R and L being solutions of the linear response equations (3) and (4) for the perturbing fields ζ and η.
The off-diagonal blocks Xibζη and Yajζη require a solution of the second-order Bethe–Salpeter response equations,
(AB)(A+B)(ωζ+ωη)21Rζη=Uζη,
(20)
(A+B)(AB)(ωζ+ωη)21Lζη=Vζη,
(21)
with the right-hand side being defined as
Uiaζη=12jRjaηMijζLjaηNijζ+(ζη)+12bRibηMabζLibηNabζ+(ζη)+Hia+[Kklζη,+]+Hia+[Kcdζη,+],
(22)
Viaζη=12jLjaηMijζRjaηNijζ+(ζη)+12bLibηMabζRibηNabζ+(ζη)+Hia[Kklζη,]+Hia[Kcdζη,].
(23)
In Eqs. (22) and (23), the quantities
Mpqζ=Upqζ+rsHpq,rs+,BSERrs,
(24)
Npqζ=Vpqζ+rsHpq,rs,BSELrs
(25)
have been introduced.
Closely following the related TD-DFT and time-dependent Hartree–Fock (HF) derivations, the first hyperpolarizability tensor β can subsequently be evaluated as11,13
vθ;vζ(ω),vη(ω)=tr(v̂θγζη)=tr(v̂θKζη)+Pθ,Qθ|Xζη,Yζη=tr(v̂θKζη)+Rθ,Lθ|Uζη,Vζη.
(26)
From β, the related quantities linked to EHSG42 and hyper-Rayleigh scattering (HRS)13 can be extracted. Note that there is no difference between TD-DFT, HF, and the static screened BSE for this expression. The same holds true for the two-photon absorption amplitudes, which can subsequently be evaluated as13,43–45
limωΩNvθ(ω);vζ(ω),vη(ωω)=RN,LN|Uζη,Vζη.
(27)
The bra in Eq. (27) denotes the eigenvector of the Nth excited state, which is purely real when real orbitals are used. As outlined in Ref. 9, even in the case of complex frequencies ω′ and ω″, only the real part of Eq. (27) therefore contributes to the observable. From the obtained 2PA amplitudes, the 2PA transition strength δ and the 2PA cross section σ can readily be calculated.13,31

The structures of water, methanol, dimethyl ether, formaldehyde, and diacetylene were optimized using the PBE046,47 functional in combination with the def2-TZVPP basis set.48 For these molecules, quasiparticle energies are subsequently obtained using either the single-shot G0W0 method or the eigenvalue-self-consistent evGW method. Within the G0W0, the quasiparticle energies are obtained from the zeroth-order approximation to Green’s function G and the screened exchange W directly, employing the Kohn–Sham (KS) orbitals.49 Within the evGW approximation, the quasiparticle energies are updated until the quasiparticle equations converge.50,51 The hybrid PBE0 and the local-hybrid TMHF52 density-functional approximations were chosen as the Kohn–Sham (KS) reference for GW-BSE due to their general good performance in previous studies.15,26,53,54 For the chromophores CFB, PYP, HBI, and HBDI, the same structures as in Ref. 31 were used. Furthermore, for these molecules, the fsCD-GW approach was used to obtain G0W0 and evGW quasiparticle energies, optimizing the highest 20 occupied and lowest 20 unoccupied states.22,26 2PA cross sections were evaluated according to Eqs. (12) and (13) of Ref. 31. A Gaussian lineshape function with a line broadening of 0.1 eV half-width at half-maximum (HWHM) was assumed. The RI approximation is used throughout with the according auxiliary basis sets.55,56 All GW and BSE calculations were performed using the aug-cc-pVXZ57–59 and d-aug-cc-pVXZ60 basis sets with X = D, T, Q, or 5, as specified in Sec. IV. All KS references were obtained using a large integration grid of size 4 or better.61 KS energies were converged to changes of less than 10−9 hartree and to root-mean-square (RMS) changes of less than 10−7 in the density matrix. All calculations were performed using a development version of Turbomole V7.8.62 

To test the GW-BSE hyperpolarizabilities, first the corresponding results are compared to the coupled-cluster results of Beaujean and Champagne outlined in Ref. 42. In the latter work, β was reported for water (H2O), methanol (MeOH), and dimethyl ether (DME). The evaluated value in the latter work is the dipole-projected isotropic hyperpolarizability β,42,
β=15ζx,y,zμζμηx,y,zβζηη+βηζη+βηηζ.
(28)
As outlined, this quantity is directly related to the electric-field-induced second-harmonic generation.42 For these molecules, experimental results are furthermore available at 1064 nm. A convenient benchmark for various methods is therefore possible, and the resulting calculated values are presented in Table I.
TABLE I.

Dynamic first hyperpolarizability β for H2O, MeOH, and DME calculated at 1064 nm. CC2,42 CCSD,42 CC3,42 and experimental values for H2O,63 MeOH,63 and DME64 were taken from the literature. All values in atomic units.

PBE0TMHF
PBE0CAMG0W0evGWG0W0evGWCC2CCSDCC3
H2
d-aug-cc-pVDZ −14.00 −12.52 −20.09 −17.68 −16.01 −14.58 −23.61 −16.50 −15.54 
d-aug-cc-pVTZ −18.98 −17.06 −21.13 −18.29 −18.22 −16.66 −27.94 −19.45 −18.77 
d-aug-cc-pVQZ −19.46 −17.74 −20.64 −18.24 −18.56 −16.98 −28.23 −19.77 −19.28 
d-aug-cc-pV5Z −19.51 −17.89 −20.60 −18.13 −18.26 −16.71 −28.04 −19.66  
Exp. −19.2 ± 0.9 
MeOH 
d-aug-cc-pVDZ −39.77 −34.27 −48.43 −43.97 −39.91 −37.41 −42.01 −33.73 −32.53 
d-aug-cc-pVTZ −40.80 −35.15 −44.29 −40.13 −36.31 −33.97 −41.6 −33.64  
d-aug-cc-pVQZ −40.72 −35.00 −42.81 −38.72 −34.82 −32.50    
Exp. −31.2 ± 1.6 
DME 
d-aug-cc-pVDZ −121.34 −104.92 −149.54 −132.72 −126.68 −117.25 −150.67 −105.43 −102.18 
d-aug-cc-pVTZ −119.66 −103.27 −133.39 −117.45 −110.67 −102.11 −138.76 −97.83  
d-aug-cc-pVQZ −119.30 −102.85 −126.27 −111.49 −105.00 −96.41    
Exp. −94.0 ± 0.25 
PBE0TMHF
PBE0CAMG0W0evGWG0W0evGWCC2CCSDCC3
H2
d-aug-cc-pVDZ −14.00 −12.52 −20.09 −17.68 −16.01 −14.58 −23.61 −16.50 −15.54 
d-aug-cc-pVTZ −18.98 −17.06 −21.13 −18.29 −18.22 −16.66 −27.94 −19.45 −18.77 
d-aug-cc-pVQZ −19.46 −17.74 −20.64 −18.24 −18.56 −16.98 −28.23 −19.77 −19.28 
d-aug-cc-pV5Z −19.51 −17.89 −20.60 −18.13 −18.26 −16.71 −28.04 −19.66  
Exp. −19.2 ± 0.9 
MeOH 
d-aug-cc-pVDZ −39.77 −34.27 −48.43 −43.97 −39.91 −37.41 −42.01 −33.73 −32.53 
d-aug-cc-pVTZ −40.80 −35.15 −44.29 −40.13 −36.31 −33.97 −41.6 −33.64  
d-aug-cc-pVQZ −40.72 −35.00 −42.81 −38.72 −34.82 −32.50    
Exp. −31.2 ± 1.6 
DME 
d-aug-cc-pVDZ −121.34 −104.92 −149.54 −132.72 −126.68 −117.25 −150.67 −105.43 −102.18 
d-aug-cc-pVTZ −119.66 −103.27 −133.39 −117.45 −110.67 −102.11 −138.76 −97.83  
d-aug-cc-pVQZ −119.30 −102.85 −126.27 −111.49 −105.00 −96.41    
Exp. −94.0 ± 0.25 

Table I reveals that especially evGW-BSE is a rather accurate approach to the dynamic second hyperpolarizability. In particular, when combined with the local-hybrid TMHF, the agreement with CC3 is remarkable. For larger basis sets, the experimental value is approached in all three cases. In contrast to this, hyperpolarizabilities obtained from TD-DFT, as represented by the PBE0 and CAM-B3LYP65 functionals, significantly deviate from the CC3 and experimental values. While CAM-B3LYP overall performs better than PBE0, it falls short compared to evGW-BSE@TMHF. G0W0-BSE consistently overestimates β, a feature that it shares with CC2. The significant overestimation of β for the latter method has already been observed earlier.42,66 While generally CC2 is assumed to be at least on par with GW-BSE for linear response properties,29 this assumption clearly is no longer valid for non-linear response properties. The performance of evGW-BSE can therefore be labeled to be surprisingly good, coming close to that of CCSD, only being significantly outperformed by CC3. Given that GW-BSE shares its computational complexity with TD-DFT, this is a gratifying result. Unfortunately, as revealed especially for DME, the GW-BSE ansatz, however, shares the strong dependence on the chosen basis set with the coupled cluster methods. This basis-set dependence is significantly more pronounced than for TD-DFT, and generally at least triple-ζ basis sets are required to yield acceptable results. Quadruple-ζ basis sets still significantly improve upon the triple-ζ results, yet the difference is considerably less pronounced than for the first step. Furthermore, it is noted that evGW is not able to reduce the dependence on the chosen KS state as efficiently as it does for excitation energies. However, deviations of up to 20% are obtained, which is in line with deviations found for oscillator strengths.67 Linear and non-linear properties are therefore more sensitive than excitation energies with respect to the chosen GW method, with the aligning effect of evGW being reduced.

The 2PA transition strengths for a selected number of states of water, for the 1A2 state of formaldehyde, and the 1Πg state of diacetylene are listed in Table II for GW-BSE, TD-DFT, CC2, and reference CC3 results.

TABLE II.

Comparison of two-photon absorption strengths δ between PBE0, CAM-B3LYP, GW-BSE, and the CC3 coupled-cluster approach. The d-aug-cc-pVTZ basis set was used, with the exception of CC3 for C4H2, where the aug-cc-pVTZ was used. CC2, CC3, and CAMB3-LYP results were taken from Ref. 7. All values in atomic units.

PBE0TMHF
StatePBE0CAMG0W0evGWG0W0evGWCC2CC3
H2
2A1 40.41 84.13 386.75 10.07 384.71 9.38 181.25 39.28 
3A1 81.59 172.39 120.45 294.05 33.60 321.61 155.90 224.95 
2B1 16.69 42.04 70.56 52.07 72.46 59.25 48.30 41.11 
1A2 41.69 47.29 38.30 30.82 39.07 33.73 55.88 42.10 
CH2
1A2 0.43 0.24 0.26 0.03 0.67 0.01 0.41 0.20 
C4H2 
g 141.57 114.08 67.12 61.07 62.42 65.74 92.27 92.61 
PBE0TMHF
StatePBE0CAMG0W0evGWG0W0evGWCC2CC3
H2
2A1 40.41 84.13 386.75 10.07 384.71 9.38 181.25 39.28 
3A1 81.59 172.39 120.45 294.05 33.60 321.61 155.90 224.95 
2B1 16.69 42.04 70.56 52.07 72.46 59.25 48.30 41.11 
1A2 41.69 47.29 38.30 30.82 39.07 33.73 55.88 42.10 
CH2
1A2 0.43 0.24 0.26 0.03 0.67 0.01 0.41 0.20 
C4H2 
g 141.57 114.08 67.12 61.07 62.42 65.74 92.27 92.61 

The agreement for 2PA transition strengths is considerably worse when compared to the hyperpolarizability results in Sec. IV A. Deviations of a few 100% are not uncommon, even between different coupled-cluster methods. The latter is rather surprising but has been reported in the literature before.7 It was later suggested that the agreement between CC2 and CC3 may, in general, be better than for the selected small molecules and that the observed large deviations are partly caused by important contributions from double excitations.31 This would also affect TD-DFT and GW-BSE results as neither of these two methods take the influence of double excitations into account by any means. Given these restrictions, the agreement of evGW-BSE with CC3 is generally good, improving upon CC2. Similar to the hyperpolarizabilities, CC2 has a tendency to overestimate 2PA transition strengths. GW-BSE in contrast has a tendency to underestimate the 2PA transition strengths. The only clear outlier from this scheme is the 3A1 state of H2O, which is shown to be extremely sensitive to the choice of basis set.7 As previously reported in Ref. 7, TD-DFT is “hit or miss,” and depending on the chosen density-functional approximation, results can span an order of magnitude. Clearly, all GW-BSE-based Ansätze can at least provide a more compact spread of resulting 2PA transition strengths, yielding an overall easier to use toolbox simply by making the choice of the underlying functional less relevant. Comparing different GW methods again outlines that evGW is only partly able to remove the dependence on the KS starting point. This can, however, be related with evGW already failing to align the corresponding excitation energies as outlined in Tables S II and S IV in the supplementary material. Consequently, pronounced variations in the obtained 2PA strengths are obtained in cases where evGW@DFT also vary significantly. As previously observed for hyperpolarizabilities, general deviations of up to 20% have to be expected even when using evGW, and the aligning effect of evGW is strongly diminished.

The choice of basis set remains rather critical for non-linear properties also within the GW-BSE method, as shown in Fig. 1, again reproducing the trends also observed for coupled-cluster-based methods.

FIG. 1.

Basis-set dependence of the two-photon absorption cross sections of the 2A1, 3A1, 2B1, and 1A2 states of water and 1Πg of diacetylene C4H2. Solid lines refer to G0W0-BSE, while dashed lines refer to evGW-BSE. The TMHF functional was used as the KS reference. The 1A2 of CH2O state was omitted due to its small 2PA transition strength.

FIG. 1.

Basis-set dependence of the two-photon absorption cross sections of the 2A1, 3A1, 2B1, and 1A2 states of water and 1Πg of diacetylene C4H2. Solid lines refer to G0W0-BSE, while dashed lines refer to evGW-BSE. The TMHF functional was used as the KS reference. The 1A2 of CH2O state was omitted due to its small 2PA transition strength.

Close modal

As shown in Fig. 1, especially the 2A1 and 3A1 states of H2O are extremely sensitive to the choice of basis set. In these cases, meaningful values are only obtained from at least doubly-augmented triple-ζ basis sets. Without a set of flexible diffuse functions, the 2PA transition strengths of such states are vastly overestimated. This trend was also observed for the coupled-cluster-based CC2, CCSD, and CC3 methods in Ref. 7. While the increased flexibility is important, additionally the 3A1 states are positioned at nearly exactly twice the energy of the 1B1 state for the aug-cc-pVTZ basis set. This leads to unphysical divergences, yielding the observed pronounced overestimation of the 2PA transition strength.68,69 For this particular state, an imaginary damping parameter should therefore be used to evaluate the dynamic polarizability. The effects of damping on the 2PA transition strength of this state are shortly outlined in the supplementary material. For the 2A1 states, evGW is more robust than G0W0, additionally outlining the importance of reasonably good quasiparticle energies. The sensitivity of the remaining states is less pronounced, but still a general tendency to overestimate 2PA transition strengths can be observed for smaller basis sets. Paterson et al. argued that for small molecular systems, the choice of basis set might, however, be more important than for more sizable systems.7 This at least yields some relief in the sense that the basis-set dependence shown in Fig. 1 is likely a worst-case situation.

Given the good performance of especially the evGW-BSE approach, further tests are conducted on the 2PA of four chromophores. These chromophores are the molecules HBI [4-(p-hydroxybenzylidene)-imidazolidin-5-one], HBDI [4-(p-hydroxybenzylidene)-1,2-dimethylimidazolin-5-one], the chromophore of the cyan fluorescent protein (CFP), and p-coumaric acid (PYP). Table III lists the computed 2PA transition strengths and the 2PA cross sections.

TABLE III.

Excitation energies, oscillator strengths fosc, 2PA transition strengths δ, and 2PA cross sections σ obtained for the chromophores CFB, PYP, HBI, and HBDI. Values obtained from the evGW-BSE@DFT/aug-cc-pVTZ method. CAM-B3LYP and CC2 values taken from Ref. 31 and EOM-CCSD values taken from Ref. 70.

CFBPYPHBIHBDI
Exc. State1A′2A′1A′1A′
evGW-BSE@PBE0 
Exc. (eV) 3.28 4.16 3.44 3.47 
fosc (a.u.) 0.58 0.64 0.64 0.67 
δ(TPA) (a.u.) 1042 1779 511 238 
σ(TPA) (GM) 6.06 16.66 3.27 1.55 
evGW-BSE@TMHF 
Exc. (eV) 3.29 4.24 3.48 3.49 
fosc (a.u.) 0.55 0.65 0.62 0.63 
δ(TPA) (a.u.) 922 2016 498 225 
σ(TPA) (GM) 5.41 19.57 3.25 1.48 
Exp. σ(TPA) (GM) 2.85  3.11  
CAM-B3LYP 
δ(TPA) (a.u.) 1659 1706 983 490 
σ(TPA) (GM) 7.46 11.57 4.85 2.44 
CC2 
δ(TPA) (a.u.) 5197 3239 3527 1281 
σ(TPA) (GM) 23.24 22.51 17.60 6.37 
EOM-CCSD 
δ(TPA) (a.u.)  2631  978 
CFBPYPHBIHBDI
Exc. State1A′2A′1A′1A′
evGW-BSE@PBE0 
Exc. (eV) 3.28 4.16 3.44 3.47 
fosc (a.u.) 0.58 0.64 0.64 0.67 
δ(TPA) (a.u.) 1042 1779 511 238 
σ(TPA) (GM) 6.06 16.66 3.27 1.55 
evGW-BSE@TMHF 
Exc. (eV) 3.29 4.24 3.48 3.49 
fosc (a.u.) 0.55 0.65 0.62 0.63 
δ(TPA) (a.u.) 922 2016 498 225 
σ(TPA) (GM) 5.41 19.57 3.25 1.48 
Exp. σ(TPA) (GM) 2.85  3.11  
CAM-B3LYP 
δ(TPA) (a.u.) 1659 1706 983 490 
σ(TPA) (GM) 7.46 11.57 4.85 2.44 
CC2 
δ(TPA) (a.u.) 5197 3239 3527 1281 
σ(TPA) (GM) 23.24 22.51 17.60 6.37 
EOM-CCSD 
δ(TPA) (a.u.)  2631  978 

Only a few clear trends can be observed when comparing the different methods. The overall spread of obtained 2PA transition strength and cross sections between different methods is significant. Comparisons between EOM-CCSD, experiment, and GW-BSE hint at the latter yielding rather promising results, yielding the least pronounced overestimation of 2PA cross sections. This is remarkable considering the moderate computational complexity of GW-BSE, being significantly lower than that of CC2 or especially EOM-CCSD. The tendency of CC2 to strongly overestimate 2PA remains in place also for these more sizable chromophores. It is therefore questionable to benchmark 2PA cross sections using CC2, as previously done in the literature.71,72

G0W0-BSE yields 2PA cross sections in close agreement with evGW-BSE, generally deviating by only a few %, as shown in the supplementary material. For this kind of calculation, the importance of the choice of the exact GW treatment is therefore largely reduced. This is in contrast to hyperpolarizabilities or even the excitation energies themselves, which are moderately dependent on the GW method, with evGW often outperforming the G0W0 method significantly. While true for the investigated systems, this behavior may, however, change if there are excited states close to the single-photon energy of the 2PA process, yielding peaks in the polarizability as outlined in Sec. IV B. The polarizability will then be strongly dependent on the underlying GW method, inducing significant changes in the calculated 2PA cross section.

A closed set of formulas needed to evaluate non-linear light–matter interactions within the static-screened Bethe–Salpeter-equation method has been outlined. The resulting expressions have subsequently been used to evaluate the first hyperpolarizabilities and two-photon absorption strengths for a set of small molecules. Comparison with high-level CC3 data as well as with experimental results revealed that especially the evGW-BSE method is an intriguingly accurate yet efficient way of assessing these non-linear properties. While possessing a similar basis-set dependence as coupled-cluster methods, employing larger basis sets is more feasible for GW-BSE than for CC2 or CCSD. Given the strong tendency of CC2 to overestimate non-linear light–matter interactions, which is further amplified by a small basis sets, GW-BSE therefore constitutes a more reliable methodical choice for the investigation of non-linear properties of moderately sized molecular systems.

Tables providing the GW-BSE excitation energies and two-photon absorption strengths for water, formaldehyde, and diacetylene; G0W0-BSE results for the four chromophores CFP, PYP, HBI, and HBDI; as well as an investigation of the impact of an imaginary damping factor on the 2PA transition strength for cases close to a resonance are contained within the supplementary material.

N.R and W.K. gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Collaborative Research Centre “4f for Future” (CRC 1573, Project No. 471424360), project Q. C.H. gratefully acknowledges the Volkswagen Stiftung for financial support.

The authors have no conflicts to disclose.

Nina Rauwolf: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Wim Klopper: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Christof Holzer: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

1.
D. P.
Shelton
and
J. E.
Rice
,
Chem. Rev.
94
,
3
(
1994
).
2.
M.
Pawlicki
,
H.
Collins
,
R.
Denning
, and
H.
Anderson
,
Angew. Chem., Int. Ed.
48
,
3244
(
2009
).
3.
H.
Sekino
and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
3022
(
1993
).
4.
H.
Larsen
,
J.
Olsen
,
C.
Hättig
,
P.
Jørgensen
,
O.
Christiansen
, and
J.
Gauss
,
J. Chem. Phys.
111
,
1917
(
1999
).
5.
F.
Pawłowski
,
P.
Jørgensen
, and
C.
Hättig
,
Chem. Phys. Lett.
391
,
27
(
2004
).
6.
S.
Høst
,
P.
Jørgensen
,
A.
Köhn
,
F.
Pawłowski
,
W.
Klopper
, and
C.
Hättig
,
J. Chem. Phys.
123
,
094303
(
2005
).
7.
M. J.
Paterson
,
O.
Christiansen
,
F.
Pawłowski
,
P.
Jørgensen
,
C.
Hättig
,
T.
Helgaker
, and
P.
Sałek
,
J. Chem. Phys.
124
,
054322
(
2006
).
8.
J. H.
Andersen
,
K. D.
Nanda
,
A. I.
Krylov
, and
S.
Coriani
,
J. Chem. Theory Comput.
18
,
6189
(
2022
).
9.
K.
Kristensen
,
J.
Kauczor
,
A. J.
Thorvaldsen
,
P.
Jørgensen
,
T.
Kjærgaard
, and
A.
Rizzo
,
J. Chem. Phys.
134
,
214104
(
2011
).
10.
S. J. A.
van Gisbergen
,
J. G.
Snijders
, and
E. J.
Baerends
,
J. Chem. Phys.
109
,
10644
(
1998
).
11.
F.
Furche
,
J. Chem. Phys.
114
,
5982
(
2001
).
12.
P.
Sałek
,
O.
Vahtras
,
J.
Guo
,
Y.
Luo
,
T.
Helgaker
, and
H.
Ågren
,
Chem. Phys. Lett.
374
,
446
(
2003
).
13.
S. M.
Parker
,
D.
Rappoport
, and
F.
Furche
,
J. Chem. Theory Comput.
14
,
807
(
2018
).
14.
X.
Blase
and
C.
Attaccalite
,
Appl. Phys. Lett.
99
,
171909
(
2011
).
15.
X.
Gui
,
C.
Holzer
, and
W.
Klopper
,
J. Chem. Theory Comput.
14
,
2127
(
2018
).
16.
X.
Blase
,
I.
Duchemin
,
D.
Jacquemin
, and
P.-F.
Loos
,
J. Phys. Chem. Lett.
11
,
7371
(
2020
).
17.
Y.
Cho
,
S. J.
Bintrim
, and
T. C.
Berkelbach
,
J. Chem. Theory Comput.
18
,
3438
(
2022
).
18.
J.
Li
,
D.
Golze
, and
W.
Yang
,
J. Chem. Theory Comput.
18
,
6637
(
2022
).
19.
Y.
Yao
,
D.
Golze
,
P.
Rinke
,
V.
Blum
, and
Y.
Kanai
,
J. Chem. Theory Comput.
18
,
1569
(
2022
).
20.
X.
Ren
,
P.
Rinke
,
C.
Joas
, and
M.
Scheffler
,
J. Mater. Sci.
47
,
7447
(
2012
).
21.
K.
Krause
and
W.
Klopper
,
J. Comput. Chem.
38
,
383
(
2017
).
22.
C.
Holzer
,
A. M.
Teale
et al,
J. Chem. Phys.
150
,
204116
(
2019
).
23.
M. P.
Ljungberg
,
P.
Koval
,
F.
Ferrari
,
D.
Foerster
, and
D.
Sánchez-Portal
,
Phys. Rev. B
92
,
075422
(
2015
).
24.
N. C.
Bradbury
,
M.
Nguyen
,
J. R.
Caram
, and
D.
Neuhauser
,
J. Chem. Phys.
157
,
031104
(
2022
).
25.
I.
Duchemin
and
X.
Blase
,
J. Chem. Theory Comput.
16
,
1742
(
2020
).
26.
M.
Kehry
,
W.
Klopper
, and
C.
Holzer
,
J. Chem. Phys.
159
,
044116
(
2023
).
27.
I.
Knysh
,
J. D. J.
Villalobos-Castro
,
I.
Duchemin
,
X.
Blase
, and
D.
Jacquemin
,
J. Phys. Chem. Lett.
14
,
3727
(
2023
).
28.
J.
Villalobos-Castro
,
I.
Knysh
,
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
,
J. Chem. Phys.
159
,
024116
(
2023
).
29.
C.
Suellen
,
R. G.
Freitas
,
P.-F.
Loos
, and
D.
Jacquemin
,
J. Chem. Theory Comput.
15
,
4581
(
2019
).
30.
J.
Kongsted
,
A.
Osted
,
K. V.
Mikkelsen
, and
O.
Christiansen
,
J. Chem. Phys.
120
,
3787
(
2004
).
31.
M. T. P.
Beerepoot
,
D. H.
Friese
,
N. H.
List
,
J.
Kongsted
, and
K.
Ruud
,
Phys. Chem. Chem. Phys.
17
,
19306
(
2015
).
32.
M. A.
Salem
and
A.
Brown
,
J. Chem. Theory Comput.
10
,
3260
(
2014
).
33.
M.
Chołuj
,
M. M.
Alam
,
M. T. P.
Beerepoot
,
S. P.
Sitkiewicz
,
E.
Matito
,
K.
Ruud
, and
R.
Zaleśny
,
J. Chem. Theory Comput.
18
,
1046
(
2022
).
34.
B.
Zerulla
,
D.
Beutel
,
C.
Holzer
,
I.
Fernandez-Corbaton
,
C.
Rockstuhl
, and
M.
Krstić
,
Adv. Mater.
(published online
2023
).
35.
S.
Ghosh
and
B.
Deb
,
Chem. Phys.
71
,
295
(
1982
).
36.
C.
Jamorski
,
M. E.
Casida
, and
D. R.
Salahub
,
J. Chem. Phys.
104
,
5134
(
1996
).
37.
X.
Blase
,
I.
Duchemin
, and
D.
Jacquemin
,
Chem. Soc. Rev.
47
,
1022
(
2018
).
38.
C.
Holzer
and
W.
Klopper
,
J. Chem. Phys.
150
,
204116
(
2019
).
39.
C.
Holzer
,
J. Chem. Theory Comput.
19
,
3131
(
2023
).
40.
D.
Rocca
,
D.
Lu
, and
G.
Galli
,
J. Chem. Phys.
133
,
164109
(
2010
).
41.
M.
Kehry
,
Y. J.
Franzke
,
C.
Holzer
, and
W.
Klopper
,
Mol. Phys.
118
,
e1755064
(
2020
).
42.
P.
Beaujean
and
B.
Champagne
,
J. Chem. Phys.
145
,
044311
(
2016
).
43.
J.
Olsen
and
P.
Jørgensen
,
J. Chem. Phys.
82
,
3235
(
1985
).
44.
O.
Christiansen
,
P.
Jørgensen
, and
C.
Hättig
,
Int. J. Quantum Chem.
68
,
1
(
1998
).
45.
C.
Hättig
,
O.
Christiansen
, and
P.
Jørgensen
,
J. Chem. Phys.
108
,
8331
(
1998
).
46.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
47.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
48.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
49.
L.
Hedin
,
J. Phys.: Condens. Matter
11
,
R489
(
1999
).
50.
M.
van Schilfgaarde
,
T.
Kotani
, and
S.
Faleev
,
Phys. Rev. Lett.
96
,
226402
(
2006
).
51.
T.
Kotani
,
M.
van Schilfgaarde
, and
S. V.
Faleev
,
Phys. Rev. B
76
,
165106
(
2007
).
52.
Y. J.
Franzke
and
C.
Holzer
,
J. Chem. Phys.
157
,
034108
(
2022
).
53.
C.
Holzer
,
Y. J.
Franzke
, and
M.
Kehry
,
J. Chem. Theory Comput.
17
,
2928
(
2021
).
54.
F.
Caruso
,
M.
Dauth
,
M. J.
van Setten
, and
P.
Rinke
,
J. Chem. Theory Comput.
12
,
5076
(
2016
).
55.
C.
Hättig
,
Phys. Chem. Chem. Phys.
7
,
59
(
2005
).
56.
F.
Weigend
,
Phys. Chem. Chem. Phys.
8
,
1057
(
2006
).
57.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
58.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
59.
D. E.
Woon
and
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
98
,
1358
(
1993
).
60.
D. E.
Woon
and
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
100
,
2975
(
1994
).
61.
O.
Treutler
and
R.
Ahlrichs
,
J. Chem. Phys.
102
,
346
(
1995
).
62.
Y. J.
Franzke
,
C.
Holzer
,
J. H.
Andersen
,
T.
Begušić
,
F.
Bruder
,
S.
Coriani
,
F.
Della Sala
,
E.
Fabiano
,
D. A.
Fedotov
,
S.
Fürst
,
S.
Gillhuber
,
R.
Grotjahn
,
M.
Kaupp
,
M.
Kehry
,
M.
Krstić
,
F.
Mack
,
S.
Majumdar
,
B. D.
Nguyen
,
S. M.
Parker
,
F.
Pauly
,
A.
Pausch
,
E.
Perlt
,
G. S.
Phun
,
A.
Rajabi
,
D.
Rappoport
,
B.
Samal
,
T.
Schrader
,
M.
Sharma
,
E.
Tapavicza
,
R. S.
Treß
,
V.
Voora
,
A.
Wodyński
,
J. M.
Yu
,
B.
Zerulla
,
F.
Furche
,
C.
Hättig
,
M.
Sierka
,
D. P.
Tew
, and
F.
Weigend
,
J. Chem. Theory Comput.
19
,
6859
(
2023
).
63.
P.
Kaatz
,
E. A.
Donley
, and
D. P.
Shelton
,
J. Chem. Phys.
108
,
849
(
1998
).
64.
V. W.
Couling
and
D. P.
Shelton
,
J. Chem. Phys.
143
,
224307
(
2015
).
65.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
66.
S.
Sirimatayanant
and
T.
Andruniów
,
J. Chem. Phys.
158
,
094106
(
2023
).
67.
D.
Jacquemin
,
I.
Duchemin
,
A.
Blondel
, and
X.
Blase
,
J. Chem. Theory Comput.
12
,
3969
(
2016
).
68.
S. M.
Parker
,
S.
Roy
, and
F.
Furche
,
J. Chem. Phys.
145
,
134105
(
2016
).
69.
D.
Dar
,
S.
Roy
, and
N. T.
Maitra
,
J. Phys. Chem. Lett.
14
,
3186
(
2023
).
70.
K. D.
Nanda
and
A. I.
Krylov
,
J. Chem. Phys.
142
,
064118
(
2015
).
71.
D. H.
Friese
,
C.
Hättig
, and
K.
Ruud
,
Phys. Chem. Chem. Phys.
14
,
1175
(
2012
).
72.
M. T. P.
Beerepoot
,
M. M.
Alam
,
J.
Bednarska
,
W.
Bartkowiak
,
K.
Ruud
, and
R.
Zaleśny
,
J. Chem. Theory Comput.
14
,
3677
(
2018
).

Supplementary Material