Simulating Third-Order Nonlinear Optical ... - ACS Publications

Subscriber access provided by UNIV OF PITTSBURGH

Article

Simulating Third-order Nonlinear Optical Properties Using Damped Cubic Response Theory within Time-Dependent Density Functional Theory Zhongwei Hu, Jochen Autschbach, and Lasse Jensen J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01060 • Publication Date (Web): 03 Feb 2016 Downloaded from http://pubs.acs.org on February 14, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Simulating Third-order Nonlinear Optical Properties Using Damped Cubic Response Theory within Time-Dependent Density Functional Theory Zhongwei Hu,† Jochen Autschbach,‡ and Lasse Jensen∗,† †Department of Chemistry, The Pennsylvania State University, 104 Chemistry Building, University Park, Pennsylvania 16802-4615, USA ‡Department of Chemistry, University at Buffalo, State University of New York, Buffalo, New York 14260-3000, USA E-mail: [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract A general implementation for damped cubic response properties is presented using time-dependent Density functional theory (TDFFT) and Slater-type orbital (STO) basis sets. To directly calculate two-photon absorption (TPA) cross sections, we also present an implementation of a reduced damped cubic response approach. Validation of the implementations includes a detailed comparison between response theory and the sum-over-states (SOS) approach for calculating the nonlinear optical (NLO) properties of LiH, as well as a comparison between the simulated and experimental TPA and third-harmonic generation (THG) spectra for the dimethyl-amino-nitrostilbene (DANS) molecule. The study of LiH demonstrates the incorrect pole structure obtained in response theory due to the adiabatic approximation typically employed for the exchange-correlation (XC) kernel. For DANS, we find reasonable agreement between simulated and experimental TPA and THG spectra. Overall, this work shows that care must be taken when calculating higher-order response functions in the vicinity of one-photon poles due to the approximate kernels typically used in the simulations.

1

Introduction

Nonlinear optical (NLO) processes involving interactions with multiple photons have seen applications in optical power-limiting 1–5 optical data storage, 6–10 and optical image processing. 11–16 There is significant interest in developing new materials for all-optical switching devices by exploiting the intensity-dependence of the refractive index (IDRI), which at the microscopic level is described by the second hyperpolarizability γ(−ω; ω, ω, −ω). 17 The requirement for all-optical switching applications is NLO materials characterized by a large Re(γ) and a small Im(γ), meaning low losses caused by weak two-photon absorption (TPA) processes. 18,19 Thus, to understand the molecular response relevant for all-optical switching, it is important to describe both the real and imaginary contributions to the third-order nonlinear response. 2


Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The characterization of the third-order response is typically done using a simplified sumover-states (SOS) approach. This approach assumes that a few states dominate the response, which provides a good description near resonances but not far from resonances. 17 While this deficiency can be remedied by including all states, it is not suitable for practical use. However, the SOS expressions can be recast in a computational efficient form that includes all electronic states, which was recently used to describe TPA at the equation-of-motion coupled-cluster level of theory. 20 Alternatively, response theory provides a general framework for deriving response functions using both approximate and exact methods without the need to explicitly consider the sum over all electronic states. 21,22 For the exact wavefunction, these different formalisms give identical results but differ when approximate methods are used. 22 Standard response theory fails to predict the correct behaviors for molecular properties in the resonance region due to the divergence of response functions when optical frequencies, or the sum of them, equal an excitation energy. Although this divergence can be exploited using residual analysis to identify transition matrix elements for multi-photon transitions, 21,22 damped response theory that takes the broadening of electronic states into account can avoid the unphysical behaviors for molecular properties on resonance. Norman et al. 23,24 pioneered the extension of standard linear and quadratic response formalisms to also describe the resonance case, where an empirical damping factor was included into the standard Ehrenfest equation. The introduction of the damping factor leads to naturally broadened absorption spectra but vibronic effects are neglected which can be important for describing linear and NLO properties. 25–28 Recently, Kristensen et al. 29 reported a quasi-energy formulation of damped response theory by using complex excitation energies. This approach is intrinsically equivalent to the one reported by Norman and co-workers, 24 and has been used to describe TPA spectra based on modified damped cubic response functions. 30 In this work, we present a general implementation of damped cubic response properties using time-dependent density functional theory (TDDFT) and Slater-type orbital (STO) basis sets in the Amsterdam Density Functional (ADF) program package. 31–33 This is an

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

extension of the damped linear 34 and quadratic 35 response theory implementations in ADF. We have implemented damped cubic response theory based on the 2n + 1 rule that leads to an efficient algorithm, as well as a direct solver for the third-order density matrix by solving a set of damped third-order response equations. The latter approach facilitates the direct partitioning of the response into contributions from localized orbitals. 36 For calculating the TPA cross sections, we also implemented a reduced damped cubic response approach similar to that presented by Kristensen et al. 29 In this work, we will only consider the pure electronic response even though the vibrational response can be significant especially when considering IDRI. 37–39 As a test of our implementation, we will present a detailed comparison between response theory and SOS approach for calculating the NLO properties of LiH. In addition to this, we will also study TPA and third-harmonic generation (THG) of dimethyl-aminonitrostilbene (DANS). The NLO properties of DANS have been extensively studied from both experimental and theoretical perspectives, 40–44 and thus it serves a good benchmark system for the current implementation.

2 2.1

Theory Damped Response Theory

In damped response theory, a complex orbital energy (εi − iΓ) is involved in the timedependent Kohn-Sham (TDKS) equation, where Γ corresponds to an energy broadening term that can be related to the finite lifetime of the excited state. 45 In our previous work, we have shown that the damped response properties can be derived from the Γ-dependent TDKS equation in the presence of an external perturbation. 35 This follows the general derivations of the response equations without damping presented by Karna and Dupuis 46 for HartreeFock theory and later adapted by Van Gisbergen within a TDDFT formalism. 47 Here we will present the main results for the damped linear, quadratic and cubic response formalisms. Assuming the perturbation independence of the basis set, 48 the damped TDKS equation is 4


Page 4 of 31

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


given by 46,47 F C − iS

∂C = SC(ε − iΓ), ∂t

(1)

subject to the orthonormality constraint ∂ † C SC = 0, ∂t

(2)

where C is the coefficient matrix of the spatial orbitals (φ = χC), S is the time-independent overlap matrix of the atomic orbitals (AOs), ε is the Lagrangian multiplier matrix, and Γ is the phenomenological damping factor. Further, F is the KS matrix expressed in the AO form as 46,47 F = h + D × (2J) + νxc ,

(3)

where h is the one-electron integral matrix, J is the Coulomb matrix, D is the density matrix, and νxc is the exchange correlation (XC) potential. If we consider a closed shell molecule interacting with an external electric field E that consists of a monochromatic oscillating part and a static part, 46,47 we can expand the density matrix D in terms of the perturbation as h i 1 n ±2iωt αβ e D (±ω, ±ω) + e±iωt D = D 0 + e±iωt D α (±ω) + D α (0) E α + 2! o αβ αβ αβ αβ × D (0, ±ω) + D (±ω, 0) + D (±ω, ∓ω) + D (0, 0) E α E β + ...,

(4)

where the superscript indicates the direction (x, y, z) of the perturbation and the number of superscripts indicates the order of the perturbation. 46,47 From the perturbed density matrix, we can obtain the molecular response properties as described below. 46 2.1.1

Linear Response

The molecular polarizability is given as 46

ααβ (∓ω; ±ω) = − Tr[H α D β (±ω)],

5


(5)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 31

where ω denotes the incident frequency and “Tr” stands for the trace. Here H α is the dipole moment matrix with elements given by

α Hst = hχs |ˆ µα |χt i,

(6)

and D β is the first-order perturbed density matrix given by

D β (±ω) = C β (±ω)nC 0† + C 0 nC β† (∓ω),

(7)

where n is the occupation number matrix and C β (±ω) represents the first-order perturbed MO coefficients. The first-order perturbed MO coefficients are related to the unperturbed MO coefficients as C β (±ω) = C 0 U β (±ω), where U β (±ω) is the first-order transformation matrix. Only the occupied-virtual block of the first-order transformation matrix is needed, for which the elements are given by β Uia (±ω) =

Gβia (±ω) , ε0a − ε0i ∓ (ω + iΓ)

(8)

where ε0a , ε0i are the KS one-electron energies of the virtual and occupied orbitals, respectively. Here Gβ (±ω) = C 0† F β (±ω)C 0 is the first-order KS matrix in the MO basis. 46 The first-order Lagrangian multiplier matrix needed for the higher-order response can be calculated from Gβ (±ω) and U β (±ω) as

εβ (±ω) = Gβ (±ω) + ε0 U β (±ω) − U β (±ω)ε0 ± (ω + iΓ)U β (±ω).

(9)

The first-order KS matrix in the AO basis is given as 47

β F β (±ω) = hβ + D β (±ω) × (2J) + νxc (±ω),

6


(10)

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


β where the elements of the first-order XC potential νxc (±ω) are 47

X β νxc (±ω) κλ = fxc (±ω) κλµν D β (±ω) µν ,

(11)

µν

for which the adiabatic local density approximation (ALDA) 49 is used for the XC kernel

fxc (±ω) κλµν ≈

Z

dr

Z

h i ALDA dr′ χκ (r)χλ (r) fxc (r, r′, 0) χµ (r′ )χν (r′ ).

(12)

However, as discussed previously, 47 these matrix elements are never actually constructed during the simulations since only the potentials are needed. Once a self-consistent solution of U β is obtained, the molecular polarizability can be calculated as described above. 2.1.2

Quadratic Response

The molecular first hyperpolarizability is given as 46

βαβγ (∓ωσ ; ±ω1 , ±ω2 ) = − Tr[H α D βγ (±ω1 , ±ω2 )],

(13)

where ωσ represents the sum of the two incident frequencies as ωσ = ω1 + ω2 , and D βγ is the second-order perturbed density matrix given by D βγ (±ω1 , ±ω2 ) = C βγ (±ω1 , ±ω2 )nC 0† + C 0 nC βγ† (∓ω1 , ∓ω2 )

(14)

+ C β (±ω1 )nC γ† (∓ω2 ) + C γ (±ω2 )nC β† (∓ω1 ), where C βγ (±ω1 , ±ω2 ) represents the second-order perturbed MO coefficients. In our previous study, we have shown that damped β can be calculated efficiently using only linear response terms by exploiting the 2n + 1 rule. 35 However, in order to calculate an even higher order response property, i.e., second hyperpolarizability γ, the quadratic response terms are also needed. This requires the calculation for the second-order transformation matrix U βγ (±ω1 , ±ω2 ) that connects the second-order perturbed MO coefficients with the unper7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 31

turbed MO coefficients as C βγ (±ω1 , ±ω2 ) = C 0 U βγ (±ω1 , ±ω2 ). Unlike the linear case which only requires the occupied-virtual block of U β , both the diagonal and off-diagonal blocks of U βγ are necessary and can be further divided into two parts. The constant part only contains the terms obtained from solving the linear response equations and is given by

βγ Uij,const (±ω1 , ±ω2 ) =

 n o  β γ γ β 1  U (±ω )U (±ω ) + U (±ω )U (±ω ) diagonal  1 2 2 1 ij ij ij ij  2          

Tijβγ (±ω1 , ±ω2 ) ε0j − ε0i ∓ (ω1 + ω2 + iΓ)

where Tijβγ (±ω1 , ±ω2 )

(15)

off-diagonal

all h X γ β Gβik (±ω1 )Ukj (±ω2 ) + Gγik (±ω2 )Ukj (±ω1 ) = k

(16)

i β γ γ β − Uik (±ω1 )εkj (±ω2 ) − Uik (±ω2 )εkj (±ω1 ) . The non-constant part reads

βγ Uij,non-const (±ω1 , ±ω2 ) =

   0            

diagonal (17)

Gβγ ij (±ω1 , ±ω2 ) 0 0 εj − εi ∓ (ω1 + ω2 +

iΓ)

off-diagonal

where Gβγ = C 0† F βγ (±ω1 , ±ω2 )C 0 is the second-order KS matrix in the MO basis. 46 The second-order Lagrangian multiplier matrix needed for the cubic response equations is given by εβγ (±ω1 , ±ω2 ) = Gβγ (±ω1 , ±ω2 ) + Gβ (±ω1 )U γ (±ω2 ) + Gγ (±ω2 )U β (±ω1 ) − U β (±ω1 )εγ (±ω2 ) − U γ (±ω2 )εβ (±ω1 ) + ε0 U βγ (±ω1 , ±ω2 ) − U βγ (±ω1 , ±ω2 )ε0 ± (ω1 + ω2 + iΓ)U βγ (±ω1 , ±ω2 ).

8


(18)

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The second-order KS matrix in the AO basis can be written as 47

βγ F βγ (±ω1 , ±ω2 ) = D βγ (±ω1 , ±ω2 ) × (2J) + νxc (±ω1 , ±ω2 ),

(19)

where the elements of the second-order XC potential are 47 X βγ νxc (±ω1 , ±ω2 ) κλ = fxc (±ωσ ) κλµν D βγ (±ω1 , ±ω2 ) µν µν

(20) XX β γ + gxc (±ω1 , ±ω2 ) κλµνστ × D (±ω1 ) µν D (±ω2 ) στ , µν

στ

for which both the first- and second-order ALDA kernels are employed and the latter one reads gxc (±ω1 , ±ω2 ) κλµνστ ≈

Z

Z

Z

dr′′ χκ (r)χλ (r) dr dr h i ALDA × gxc (r, r′ , r′′, 0, 0) χµ (r′ )χν (r′ )χσ (r′′ )χτ (r′′ ). ′

(21)

Once a self-consistent solution for the non-constant part of U βγ is obtained, the first hyperpolarizability can be calculated as described above by summing up both the constant and non-constant parts. 2.1.3

Cubic Response

The molecular second hyperpolarizability is given as 46

γαβγδ (∓ωσ ; ±ω1 , ±ω2 , ±ω3 ) = − Tr[H α D βγδ (±ω1 , ±ω2 , ±ω3 )],

(22)

where ωσ represents the sum of the three incident frequencies as ωσ = ω1 + ω2 + ω3 , and D βγδ is the third-order perturbed density matrix. To obtain the damped second-order hyperpolarizability γ, we could directly solve for the third-order density matrix by solving a set of damped third-order CPKS equations. However, this requires the solution of several cubic

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

response equations. Therefore, it becomes computational advantageous to use the 2n + 1 rule 50 to rewrite the cubic response in terms of quantities that can be obtained by only solving linear and quadratic response equations. 46,51 Here, we adopt the latter approach for calculating the damped γ but we note that the direct solution of the third-order density matrix have also been implemented. Using the 2n + 1 rule, the final expression for the damped γ can be written as γαβγδ (∓ωσ ; ±ω1 , ±ω2 , ±ω3 ) = h X n Tr n U α (∓ωσ )Gβ (±ω1 )U γδ (±ω2 , ±ω3 ) + U β (±ω1 )Gα (∓ωσ )U γδ (±ω2 , ±ω3 ) P

− U α (∓ωσ )U β (±ω1 )εγδ (±ω2 , ±ω3 ) − U β (±ω1 )U α (∓ωσ )εγδ† (∓ω2 , ∓ω3 ) + U α (∓ωσ )Gγδ (±ω2 , ±ω3 )U β (±ω1 ) + U β (±ω1 )Gγδ (±ω2 , ±ω3 )U α (∓ωσ ) − U α (∓ωσ )U γδ (±ω2 , ±ω3 )εβ (±ω1 ) − U β (±ω1 )U γδ (±ω2 , ±ω3 )εα (∓ωσ ) + U γδ† (∓ω2 , ∓ω3 )U α (∓ωσ )εβ (±ω1 ) + U γδ† (∓ω2 , ∓ω3 )U β (±ω1 )εα (∓ωσ ) − U γδ† (∓ω2 , ∓ω3 )Gα (∓ωσ )U β (±ω1 ) − U γδ† (∓ω2 , ∓ω3 )Gβ (±ω1 )U α (∓ωσ ) h + Tr hxc (r, r′ , r′′ , r′′′ , ±ω1 , ±ω2 , ±ω3 )D α (∓ωσ )D β (+ω1 )D γ (±ω2 )D δ (±ω3 ) oi X n + gxc (r, r′ , r′′ , ±ω2 , ±ω3 )D α (∓ωσ )D β (±ω1 )D γδ (±ω2 , ±ω3 ) ,

oi

P

where

(23)

P

represents a summation over corresponding terms obtained by permuting (±ω1 , β) P and (±ω2 , ±ω3 , γδ). For example, P U α (∓ωσ )Gβ (±ω1 )U γδ (±ω2 , ±ω3 ) is equivalent to h i α β γδ γ βδ δ βγ U (∓ωσ ) G (±ω1 )U (±ω2 , ±ω3 )+G (±ω2 )U (±ω1 , ±ω3 )+G (±ω3 )U (±ω1 , ±ω2 ) , where P

(∓ωσ , α) is not involved. It is important to note that the relationship εαβ (+ω1 , +ω2 ) =

εαβ† (−ω1 , −ω2 ) given in Ref. 46 does not hold when the damping factor is included. Instead, the following is true εαβ† (−ω1 , −ω2 ) = εαβ (+ω1 , +ω2 ) + iΓW αβ (+ω1 , +ω1 ), where W αβ (+ω1 , +ω2 ) = U α (+ω1 )U β (+ω2 ) + U β (+ω2 )U α (+ω1 ). See the Supporting Information for the derivation. For Equation (23) above, most terms can be obtained directly by solving either the linear or quadratic response functions. The only one that requires additional cal-

10


Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


culation is the third-order XC kernel, hxc (r, r′ , r′′ , r′′′ , ±ω1 , ±ω2 , ±ω3 ), which can be acquired using ALDA as h i Tr hxc (r, r′ , r′′ , r′′′ , ±ω1 , ±ω2 , ±ω3 )D α (∓ωσ )D β (±ω1 )D γ (±ω2 )D δ (±ω3 ) = Z d3 rhALDA (r)ρα (r, ∓ωσ )ρβ (r, ±ω1 )ργ (r, ±ω2 )ρδ (r, ±ω3 ), xc

(24)

where ρα , ρβ , ργ and ρδ are the first-order perturbed densities from the damped linear response theory. 47

2.2

The SOS Approach

In the SOS approach, the finite lifetime of excited states is introduced in terms of the complex excitation energy. Adopting the phenomenological damping factor, the complex excitation energy can be obtained as Ei → Ei − iΓ. We note that this Γ used here is equivalent to that adopted for the response formalism as all single-particle states are assumed to have a common energy broadening parameter in damped response theory. 34,35 Consequently, the damped first-, second- and third-order response properties can be calculated as 21,22,41,52–55 1 X h0| µα |mi hm| µβ |0i h0| µβ |mi hm| µα |0i ααβ (−ω; ω) = , + h ¯ (ωm0 − ω − iΓ) (ωm0 − ω + iΓ)

(25)

m6=0

1 βαβγ (−ωσ ; +ω1 , +ω2) = 2 P (β, γ; +ω1, +ω2 ) h ¯ X X h0| µα |mi hm| µγ |ni hn| µβ |0i h0| µγ |mi hm| µβ |ni hn| µα |0i + × (ωm0 − ωσ − iΓ)(ωn0 − ω1 − iΓ) (ωm0 + ω2 + iΓ)(ωn0 + ωσ + iΓ) m6=0 n6=0 h0| µγ |mi hm| µα |ni hn| µβ |0i , + (ωm0 + ω2 + iΓ)(ωn0 − ω1 − iΓ)

11


(26)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and

Page 12 of 31

1 γαβγδ (−ωσ ; +ω1, +ω2 , +ω3 ) = 3 P (β, γ, δ; +ω1, +ω2, +ω3 ) h ¯ X XX α h0| µ |mi hm| µδ |ni hn| µγ |pi hp| µβ |0i × (ωm0 − ωσ − iΓ)(ωn0 − ω1 − ω2 − iΓ)(ωp0 − ω1 − iΓ) m6=0 n6=0 p6=0 h0| µδ |mi hm| µα |ni hn| µγ |pi hp| µβ |0i (ωm0 + ω3 + iΓ)(ωn0 − ω1 − ω2 − iΓ)(ωp0 − ω1 − iΓ) h0| µβ |mi hm| µγ |ni hn| µα |pi hp| µδ |0i + (ωm0 + ω1 + iΓ)(ωn0 + ω1 + ω2 + iΓ)(ωp0 − ω3 − iΓ)

+

h0| µβ |mi hm| µγ |ni hn| µδ |pi hp| µα |0i + (ωm0 + ω1 + iΓ)(ωn0 + ω1 + ω2 + iΓ)(ωp0 + ωσ + iΓ) XX h0| µα |mi hm| µδ |0i h0| µγ |ni hn| µβ |0i − (ωm0 − ωσ − iΓ)(ωm0 − ω3 − iΓ)(ωn0 − ω1 − iΓ) m6=0 n6=0

(27)

h0| µα |mi hm| µδ |0i h0| µγ |ni hn| µβ |0i (ωm0 − ω3 − iΓ)(ωn0 + ω2 + iΓ)(ωn0 − ω1 − iΓ) h0| µδ |mi hm| µα |0i h0| µβ |ni hn| µγ |0i + (ωm0 + ωσ + iΓ)(ωm0 + ω3 + iΓ)(ωn0 + ω1 + iΓ) h0| µδ |mi hm| µα |0i h0| µβ |ni hn| µγ |0i , + (ωm0 + ω3 + iΓ)(ωn0 − ω2 − iΓ)(ωn0 + ω1 + iΓ)

+

respectively. Here, P (β, γ, δ; +ω1, +ω2 , +ω3) is a permutation operator that takes all the possible combinations among (β, +ω1), (γ, +ω2 ) and (δ, +ω3 ) into account, h0| µ |mi is the transition dipole moment between the ground state and the mth excited state, hm| µ |ni = hm| µ |ni − h0| µ |0i is the fluctuation of the dipole moment between mth and nth excited states from its ground state value, and ωm0 is the transition frequency that relates to the excitation energy between the ground and mth excited states as h ¯ ωm0 = Em0 = Em − E0 .

2.3

Two-photon absorption cross section

Assuming two degenerate photons that are linearly polarized, 56 the TPA cross section (σTPA ) can be calculated using the imaginary part of the third-order response property as 28,57 Nπ 3 αf2 ω 2h ¯3 X h σTPA (ω) = Im γααββ (−ω; ω, ω, −ω) 15e4 αβ i + γαββα (−ω; ω, ω, −ω) + γαβαβ (−ω; ω, ω, −ω) , 12


(28)

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where αf is the fine structure constant, and the σTPA unit is given as Göppert-Mayer (1 GM = 10−50 cm4 s photon−1 ). 58 The integer value N is related to the experimental setup and in this work N = 4 is used for all simulated TPA spectra. 57 The specific second hyperpolarizability used in Eq. (28) governs both TPA and IDRI processes, and in the vicinity of one-photon poles the third-order response becomes negative due to saturated absorption. 54,55 In the SOS approach, the expression for IDRI second hyperpolarizability can be simplified if ω is far from any one-photon resonances as

TPA γαβγδ (−ω; ω, ω, −ω)

βγ αδ ω − iΓ, ω − iΓ S0n ω + iΓ, ω + iΓ 1 X S0n = 3 , ωn0 − 2ω − iΓ h ¯ n

(29)

where we have introduced the TPA transition moments

αδ S0n

ω ± iΓ, ω ± iΓ =

X h0| µα |mi hm| µδ |ni m

ωm0 − (ω ± iΓ)

h0| µδ |mi hm| µα |ni . + ωm0 − (ω ± iΓ)

(30)

The full derivation is presented in the Supporting Information, but it is important to note that this simplification requires the elimination of many terms from Eq. (27) that becomes important in the vicinity of one-photon poles. Far from any one-photon resonances, the damping terms in Eq. (30) can be neglected and the TPA transition moments can be calculated using standard response theory.

21,22

However, standard response functions are

divergent in the case of a double resonance where there is a one-photon resonance at half the energy of the TPA transition. The damped response formalism presented by Kristensen et al. 29 allows for the calculation of TPA intensities also in the case of double resonance effect by using a modified damped cubic response formalism. Using the same strategy, the reduced IDRI second hyperpolarizability can be obtained using the damped response theory

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 31

framework presented in this work as TPA γαβγδ (−ω; ω, ω, −ω) = h n Tr n U α (−ω)Gδ (−ω)U βγ (ω, ω) + U δ (−ω)Gα (−ω)U βγ (ω, ω)

− U α (−ω)U δ (−ω)εβγ (ω, ω) − U δ (−ω)U α (−ω)εβγ† (−ω, −ω) + U α (−ω)Gβγ (ω, ω)U δ (−ω) + U δ (−ω)Gβγ (ω, ω)U α (−ω) (31) α

βγ

δ

δ

βγ

α

− U (−ω)U (ω, ω)ε (−ω) − U (−ω)U (ω, ω)ε (−ω) + U βγ† (−ω, −ω)U α (−ω)εδ (−ω) + U βγ† (−ω, −ω)U δ (−ω)εα (−ω) − U βγ† (−ω, −ω)Gα (−ω)U δ (−ω) − U βγ† (−ω, −ω)Gδ (−ω)U α (−ω) h i + Tr gxc (r, r′, r′′ , ω, ω)D α(−ω)D δ (−ω)D βγ (ω, ω) ,

oi

where only the two-photon terms are included. This reduced form is consistent with the SOS expression given in Eq. (29), and will be used for all simulated TPA spectra in this work unless specified.

3

Computational Details

We have implemented the damped cubic response theory into a locally modified version of the ADF program package. 31–33 Geometry optimization was performed using the Becke-Perdew (BP86) 59,60 XC functional with a triple-ζ polarized Slater type (TZP) basis set from the ADF library. Response properties of LiH were calculated using the local density approximation (LDA) through ADF with a single-ζ (SZ) Slater type basis set, and Dalton 61 with a STO-3G basis set. For LiH, we also performed full configuration interaction (FCI) calculations using Dalton with the STO-3G basis set. Response properties of DANS were calculated using ADF with the statistical average of orbital model exchange-correlation potentials (SAOP) 62,63 and a triple-ζ polarized Slater type (TZP) basis set. To account for the incomplete fit basis sets in ADF, the keyword “FitAOderiv” was invoked for all response calculations. Unless

14


Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


specified, the finite lifetime of the electronic excited states is included phenomenologically using a damping parameter of 0.0034 a.u. (∼800 cm−1 ), which was previously found to be acceptable. 34,35 The isotropic average of the response properties was calculated for the polarizability as

α(−ω; ¯ ω) =

1 αxx (−ω; ω) + αyy (−ω; ω) + αzz (−ω; ω) , 3

(32)

for the first hyperpolarizability as 64,65 1 X ¯ β(−ω ; +ω , +ω ) = βzαα (−ωσ ; +ω1 , +ω2 ) σ 1 2 5 α

(33)

+ βαzα (−ωσ ; +ω1 , +ω2 ) + βααz (−ωσ ; +ω1 , +ω2 ) , and for the second hyperpolarizability as 64,66,67 1 X γααββ (−ωσ ; +ω1 , +ω2 , +ω3) γ¯ (−ωσ ; +ω1 , +ω2 , +ω3) = 15 αβ + γαββα (−ωσ ; +ω1 , +ω2, +ω3 ) + γαβαβ (−ωσ ; +ω1 , +ω2, +ω3 ) . (34) Unless otherwise stated, in the following all SOS and response calculations have been done using Dalton and ADF, respectively.

4 4.1

Results and Discussion Comparing response theory and the SOS approach for LiH

As an initial test of our implementation, we will present a detailed comparison between response theory and the SOS approach for calculating the NLO properties of LiH. The excitation energies, transition dipole moments, and dipole fluctuations operators needed in the SOS calculations were obtained using a STO-3G basis set and the LDA XC-functional.

15



While this minimum basis set will not provide accurate results for the NLO properties, it allows for a straightforward inclusion of all excited states in the SOS approach. In the Supporting Information, we show that the α, β, and γ obtained using response theory with the SZ basis set are in good agreement with those obtained from the SOS approach using the STO-3G basis set. It is important to note that the agreement between response theory and the SOS approach for NLO properties is not exact at the level of TDDFT since the manifold of excited states is not explicitly resolved. Therefore, we also include SOS results obtained at the FCI level using the STO-3G basis set. For the FCI results, it is shown in the Supporting Information that there is near exact agreement between response theory and the SOS approach as expected for an exact theory.

100 80 60 40 20 0 −20 −40 (a)

σTPA ( × 10−50 cm4 s/photon)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

0.0 100 80 60 40 20 0 −20 −40 (b) 0.0 100 80 60 40 20 0 −20 −40 (c) 0.0

TDDFT Response one photon pole two photon pole

0.5

1.0

×10−3

1.5

2.0

2.5

3.0

TDDFT SOS one photon pole two photon pole

0.5

4.0

4.5

4.0

4.5

4.0

4.5

×10−1

1.0

1.5

2.0

2.5

3.0

FCI SOS one photon pole two photon pole

0.5

3.5

3.5

×10−1

1.0

1.5

2.0 2.5 Energy (eV)

3.0

3.5

Figure 1: Simulated TPA spectra for LiH using the full expression of γIDRI through (a) TDDFT response theory, (b) TDDFT SOS approach, and (c) FCI SOS approach. The vertical “- -” and “- .” lines indicate the one- and two-photon resonances due to the excited states, respectively. Values in the shaded region are scaled by the labeled factor. In Figure 1 we compare the σTPA calculated using response theory at the TDDFT level, the SOS approach at the TDDFT level, and the SOS approach at the FCI level. The 16


Page 17 of 31

calculation of σTPA has been done using the full expression of γ IDRI and thus also contains the one-photon processes. The σTPA for excitation energies up to 4.5 eV are shown, which covers the two lowest two-photon transitions and the lowest one-photon transition. This shows that good agreement is found for all three approaches in the two-photon region below 3.0 eV. However, in the vicinity of the one-photon resonance, we find that the result obtained using response theory differs significantly from the SOS results. The magnitude calculated using response theory is two orders of magnitude larger than the SOS results and the band shape is also very different. Consider the good agreement between the two SOS spectra with respect to both the magnitude and the band shape, this finding indicates a potential issue in TDDFT response theory that will be addressed below.

120 105 90 75 60 45 30 15 (a)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.0 120 105 90 75 60 45 30 15 (b) 0.0 120 105 90 75 60 45 30 15 (c) 0.0

TDDFT Response one photon pole two photon pole

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

1.5

2.0

2.5

3.0

3.5

4.0

4.5

1.5

2.0 2.5 Energy (eV)

3.0

3.5

4.0

4.5

TDDFT SOS one photon pole two photon pole

0.5

1.0

FCI SOS one photon pole two photon pole

0.5

1.0

Figure 2: Simulated TPA spectra for LiH using the reduced expression of γIDRI (γ TPA ) through (a) TDDFT response theory, (b) TDDFT SOS approach, and (c) FCI SOS approach. The vertical “- -” and “- .” lines indicate the one- and two-photon resonances due to the excited states, respectively. As the full expression of γ IDRI includes one-photon processes that are large and negative near one-photon resonances, it is interesting to perform a similar comparison using the 17



reduced response functions γ TPA that only contains the TPA process. The comparison among the TPA dominant terms calculated using response theory and the SOS approach at the TDDFT level and the SOS approach at the FCI level is shown in Figure 2. As expected, the two-photon regions are nearly identical to those obtained using the full expression and show good agreement between response theory and the SOS approach. The one-photon region above 3.0 eV is now all positive as per construction of the reduced response function, however, the σTPA obtained using response theory is still about an order of magnitude larger than that obtained using the SOS approach. 4

4 Real Imag one photon pole

2 1 0 −1 −2 −3 0.0

Real Imag one photon pole

3 γ¯OKE(−ω; ω, 0, 0) (×107 a.u.)

3 γ¯OKE(−ω; ω, 0, 0) (×107 a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

2 1 0 −1 −2 −3

0.5

1.0

1.5

2.0 2.5 3.0 Energy (eV)

3.5

4.0

4.5

0.0

0.5

1.0

(a) TDDFT Response

1.5

2.0 2.5 3.0 Energy (eV)

3.5

4.0

4.5

(b) TDDFT SOS

Figure 3: Simulated OKE spectra for LiH using (a) response theory and (b) the SOS approach at the TDDFT level of theory. The vertical “- -” line indicates the one photon resonance due to the excited state. We found similar trends when comparing response theory with the SOS approach for other NLO processes involving multiple photons, such as second harmonic generation (SHG), electric field induced SHG (EFISHG), and THG. In all cases, good agreement is found except in the vicinity of one-photon poles. In opposition to this, NLO processes involving only one photon, such as electrooptic Pockels effect (EOPE), optical rectification (OR), optical Kerr effect (OKE), and electric field induced OR (EFIOR), all show good agreement between response theory and the SOS approach at all wavelengths. Part of this is illustrated in Figure 3, where we compare averaged γOKE (−ω; ω, 0, 0) spectra obtained by response theory 18


Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and the SOS approach at the TDDFT level of theory. The comparison between response theory and the SOS approach for the other NLO processes can be found in the Supporting Information. To further understand in detail the origin of the differences in the NLO properties obtained using response theory and the SOS approach near one-photon resonances, we will for simplicity consider the SHG process. Near a one-photon resonance, the SHG response is dominated by the following terms in the SOS approach SOS βSHG

X X h0| µα |mi hm| µβ |ni hn| µγ |0i h0| µα |mi hm| µγ |ni hn| µβ |0i ≈ + (ωm0 − 2ω − iΓ)(ωn0 − ω − iΓ) (ωm0 − 2ω − iΓ)(ωn0 − ω − iΓ) m6=0 n6=0

(35) h0| µγ |mi hm| µα |ni hn| µβ |0i h0| µβ |mi hm| µα |ni hn| µγ |0i . + + (ωm0 + ω + iΓ)(ωn0 − ω − iΓ) (ωm0 + ω + iΓ)(ωn0 − ω − iΓ)

If there is no double resonance effect (i.e., 2ω → ωm0 and ω → ωn0 ), as is the case for LiH, then we would expect the SHG to be dominated by a single resonance term determined by (ωn0 − ω − iΓ)−1 when ω → ωn0 . For comparison, we also derived the spectral representation of the SHG response function as obtained from damped response theory. The full details are given in the Supporting Information. The one-photon dominant term is found to be

Response βSHG

occ X virt X all exc. ∗ HXC Xno. Ym,ik XX µα0m,ik · Kkl,jb · Xn,jbµβn0,jb · Xp,li µγp0,li ≈n (ωm − 2ω − iΓ)(ωn − ω − iΓ)(ωp − ω − iΓ) p ij kl mnp6=0 b

−

∗ HXC Ym,ik µα0m,ik · Xp,kl µγp0,kl · Kli,jb · Xn,jbµβn0,jb (ωm − 2ω − iΓ)(ωp − ω − iΓ)(ωn − ω − iΓ)

∗ HXC Ym,ik µα0m,ik · Kkl,bj · Yn,jbµβn0,jb · Xp,li µγp0,li + (ωm − 2ω − iΓ)(ωn − ω − iΓ)(ωp − ω − iΓ)

(36)

∗ HXC Ym,ik µα0m,ik · Xp,klµγp0,kl · Kli,bj · Yn,jbµβn0,jb , − (ωm − 2ω − iΓ)(ωp − ω − iΓ)(ωn − ω − iΓ) where Xn and Yn are the spectral representation of the response vectors obtained from the linear response equations. It is important to note that the coupling matrix K HXC , as introduced in section 2.1.1, includes the adiabatic XC kernel and thus is frequency independent. The spectral representation reveals that response theory shows a different pole structure 19



than that found in the SOS approach, where the one-photon dominant term contains only one pole for the latter one. This is further illustrated in Figure 4, where we compare the SHG response obtained using damped response theory and the SOS approach. Similar to the cubic response properties discussed above, we find good agreement between the two approaches in the two-photon region below 3.0 eV. 1.6

0.8

1.6 Real Imag one photon pole two photon pole

1.2 β¯SHG(−2ω; ω, ω) (×104 a.u.)

1.2 β¯SHG(−2ω; ω, ω) (×104 a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

0.4 0.0 −0.4 −0.8 −1.2 0.0

0.8

Real Imag one photon pole two photon pole

0.4 0.0 −0.4 −0.8 −1.2

0.5

1.0

1.5

2.0 2.5 3.0 Energy (eV)

3.5

4.0

4.5

0.0

0.5

1.0

(a) TDDFT Response

1.5

2.0 2.5 3.0 Energy (eV)

3.5

4.0

4.5

(b) TDDFT SOS

Figure 4: Simulated SHG spectra for LiH using (a) response theory and (b) the SOS approach at the TDDFT level of theory. The vertical “- -” and “- .” lines indicate the one- and twophoton resonances due to the excited states, respectively. In the one-photon region, we find that the band structure differs significantly between the two approaches. For the SOS approach, the band shape is Lorentzian as expected based on the one-photon pole structure given in Eq. (35). The band structure obtained using damped response theory is more complicated due to the additional pole found in the spectral presentation of the response function. Thus, the differences between damped response theory and the SOS approach can be traced back to the different pole structures of the response functions. This incorrect pole structure obtained in response theory has previously been recognized by Dalgaard 68 in the context of time-dependent Hartree Fock (TDHF) response theory. Within TDDFT, the spurious pole arises from the adiabatic approximation that renders the K HXC frequency independent. Such deficiencies of the TDDFT response function have also recently been pointed out when calculating the derivative couplings between excited 20


Page 21 of 31

states. 69–71

4.2

TPA and THG of DANS

To further test our implementation of damped cubic response theory, we consider TPA and THG of DANS. The TPA of DANS has been measured experimentally using the femtosecond Z-scan technique, 42 and THG has been determined experimentally for the DANS chromophore embedded in a polymer matrix. 41 Therefore, DANS serves as a good benchmark system here. To better approximate the observed band width found in the experimental spectra, 42 we used a larger damping factor of Γ = 0.0068 a.u. for the simulations. In 4.0

4.0 Theo. Absorption

3.5 3.0 2.5 2.0 1.5

3.0 2.5 2.0 1.5

1.0

1.0

0.5

0.5

0.0

0.5

1.0

1.5 2.0 2.5 Energy (eV)

Exp. Absorption

3.5 ε (×104 M−1 cm−1)

ε (×104 M−1 cm−1 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.0

3.5

4.0

1.0

1.5

2.0

(a) Theory

2.5 3.0 3.5 Energy (eV)

4.0

4.5

5.0

(b) Experiment

Figure 5: (a) Simulated absorption spectrum for DANS in the gas phase. (b) Experimental absorption spectrum for DANS in DMSO taken from Ref. 42. Figure 5 we compare the simulated and experimental absorption spectra for DANS. The simulated spectrum is characterized by two main bands, which are found at 2.15 eV and 3.29 eV with the lowest band being relatively weaker. In contrast to this, the experimental spectrum shows opposite relative intensities for the two bands located at 2.75 eV and 4.10 eV, respectively. Such a discrepancy is likely caused by the SAOP XC potential used in the simulation. In addition to this, a weak shoulder found around 3.70 eV in the simulated spectrum is also not seen in the experimental spectrum. The underestimation of the excitation energies in the simulation is expected as the lowest excitation of DANS was found to be of 21



charge-transfer character. 42 Even though the solvent effects that would likely have caused a red shift 72 were not included in the simulation, the lowest band is still found red shifted by about 0.6 eV related to the lowest band in the experimental spectrum. In Figure 6(a) we plot the simulated TPA spectrum of DANS in the gas phase and compare with the experimental spectrum in Figure 6(b) obtained by Antonov et al. 42 in DMSO. The simulated spectrum shows a dominant TPA band at 1.09 eV compared to the 400

400 Theo. γTPA

350

S1 two photon pole


350 σTPA ( × 10−50 cm4 s/photon)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

300 250 200 150 100

300 250 200 150 100 50

50 0.7

Exp. γTPA S1 two photon pole

0.8

0.9

1.0 1.1 Energy (eV)

1.2

1.3

1.4

1.0

1.1

(a) Theory

1.2

1.3 1.4 Energy (eV)

1.5

1.6

1.7

(b) Experiment

Figure 6: (a) Simulated TPA spectrum for DANS in the gas phase. (b) Experimental TPA spectrum for DANS in DMSO taken from Ref. 42. The vertical “- -” line in (a) and (b) indicates the two-photon resonance due to the calculated and measured S1 , respectively. experimental TPA band at 1.37 eV. The maximum cross section of the simulated band is found around 229 GM, which is in good agreement with its experimental counterpart measured as 190 ± 10 GM. 42 However, this is potentially fortuitous due to the XC potential used in the simulation as well as the neglect of solvent effects. Furthermore, it is important to note that the intensities depend strongly on the damping parameter used in the simulations, and thus it is important to choose appropriate values. Here we choose to fix the Γ value based on the experimental absorption spectrum. The THG measurement of DANS (as a side-chain polymer containing the chromophore) was carried out in a vacuum cell by means of the Maker fringes technique, and the reported THG spectrum was normalized with respect to the maximum γ(3ω) tensor. 41 Thus for 22


Page 23 of 31

consistency, in Figure 7(a) we normalize the simulated THG spectrum of gas-phase DANS and compare with the experimental spectrum obtained by Beljonne et al. 41 as shown in Figure 7(b). Both simulated and experimental spectra are dominated by two main bands, 1.0

1.0 Exp. Norm. γTHG

Theo. Norm. γTHG

S1 two photon pole

S1 two photon pole

0.8

S1 three photon pole

0.8 THG/THGmax

0.6

0.4

0.2

0.0

S1 three photon pole S2 three photon pole

S2 three photon pole

THG/THGmax

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.6

0.4

0.2

0.2

0.4

0.6 0.8 1.0 Energy (eV)

1.2

1.4

1.6

0.0

0.2

0.4

(a) Theory

0.6 0.8 1.0 Energy (eV)

1.2

1.4

1.6

(b) Experiment

Figure 7: (a) Simulated THG spectrum for DANS in the gas phase. (b) Experimental THG spectrum for DANS taken from Ref. 42. The vertical “- -”, “- .”, and “. .” lines in (a)/(b) indicate the two-photon resonance due to the calculated/measured S1 , threephoton resonance due to the calculated/measured S1 , and three-photon resonance due to the calculated/measured S2 , respectively. The dashed curve in (b) is manually fitted for experimental data and shown for eye guide. Note that the experimental two- and threephoton resonance energies are all obtained according to Ref. 40. locating at 0.72 eV and 1.10 eV for the former one while 0.94 eV and 1.29 eV for the latter one. The differences in band position are again attributed to the XC potential in conjunction with the missing environmental effects in the simulation. In the simulated spectrum, the second band is found to be the most intense whereas in the experimental spectrum the lowest band is found to be the strongest. Part of this difference can be attributed to the stronger oscillator strength found in the simulation for the second absorption band. However, we find it is more interesting that there is a potential stronger double resonance effect in the simulated spectrum as indicated in Figure 7(b). This double resonance effect is caused by a three-photon resonance with S2 and a two-photon resonance with S1 . In the experimental spectrum, the overlap between these two resonances is weak due to the energy separation between the two lowest excited states and thus will cause a less pronounced double resonance 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

effect. Therefore, for accurately describing these higher-order response properties, it becomes important to not only capture the absolute band positions but also their relative positions.

5

Conclusions

In this work, we have presented a general implementation of the damped cubic response properties within a TDDFT framework. To calculate the TPA cross sections, we also have presented a reduced damped cubic response approach. By providing a detailed comparison between response theory and SOS approach for calculating the NLO properties of LiH, we illustrate the inconsistent behaviors for multi-photon NLO processes near one-photon resonances and demonstrate the incorrect pole structure obtained in response theory. The root cause of this is the adiabatic approximation that renders the XC kernel frequency independent. The TPA and THG spectra of DANS were simulated and compared with experimental spectra, where reasonable agreement was found in general. More importantly, this work shows that care must be taken when calculating higher-order response functions in the vicinity of one-photon poles when the adiabatic approximation is used.

Acknowledgement L.J. acknowledges support from the NSF award CHE-1362825 and support received from Research Computing and Cyberinfrastructure, a unit of Information Technology Services at Penn State. J.A. acknowledges support from the National Science Foundation, Grant No. CHE-1265833.

Supporting Information Available All detailed derivations as well as comparisons between response theory and the SOS approach. This material is available free of charge via the Internet at http://pubs.acs.org/. 24


Page 24 of 31

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Tutt, L. W.; Kost, A. Nature 1992, 356, 225–226. (2) Ehrlich, J. E.; Wu, X. L.; Lee, I.-Y. S.; Hu, Z.-Y.; Röckel, H.; Marder, S. R.; Perry, J. W. Opt. Lett. 1997, 22, 1843–1845. (3) W. Spangler, C. J. Mater. Chem. 1999, 9, 2013–2020. (4) Neto, N. M. B.; Oliveira, S. L.; Misoguti, L.; Mendonca, C. R.; Goncalves, P. J.; Borissevitch, I. E.; Dinelli, L. R.; Romualdo, L. L.; Batista, A. A.; Zilio, S. C. J. Appl. Phys. 2006, 99, 123103. (5) Hales, J. M.; Cozzuol, M.; Screen, T. E. O.; Anderson, H. L.; Perry, J. W. Opt. Express 2009, 17, 18478–18488. (6) Parthenopoulos, D. A.; Rentzepis, P. M. Science 1989, 245, 843–845. (7) Strickler, J. H.; Webb, W. W. Opt. Lett. 1991, 16, 1780–1782. (8) Shen, Y.; Swiatkiewicz, J.; Jakubczyk, D.; Xu, F.; Prasad, P. N.; Vaia, R. A.; Reinhardt, B. A. Appl. Opt. 2001, 40, 938–940. (9) Corredor, C. C.; Huang, Z.-L.; Belfield, K. D.; Morales, A. R.; Bondar, M. V. Chem. Mater. 2007, 19, 5165–5173. (10) Dvornikov, A. S.; Walker, E. P.; Rentzepis, P. M. J. Phys. Chem. A 2009, 113, 13633– 13644. (11) Köhler, R. H.; Cao, J.; Zipfel, W. R.; Webb, W. W.; Hanson, M. R. Science 1997, 276, 2039–2042. (12) Miller, M. J.; Wei, S. H.; Parker, I.; Cahalan, M. D. Science 2002, 296, 1869–1873.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 31

(13) Zipfel, W. R.; Williams, R. M.; Christie, R.; Nikitin, A. Y.; Hyman, B. T.; Webb, W. W. Proc. Natl. Acad. Sci. USA 2003, 100, 7075–7080. (14) Kneipp, J.; Kneipp, H.; Kneipp, K. Proc. Natl. Acad. Sci. USA 2006, 103, 17149. (15) Fuentes-Hernandez, C.; Ramos-Ortiz, G.; Tseng, S.-Y.; Gaj, M. P.; Kippelen, B. J. Mater. Chem. 2009, 19, 7394–7401. (16) Clark, T. B.; Orr, M. E.; Flynn, D. C.; Goodson, T. J. Phys. Chem. C 2011, 115, 7331–7338. (17) Kanis, D. R.; Ratner, M. A.; Marks, T. J. Chem. Rev. 1994, 94, 195. (18) Mizrahi, V.; Saifi, M. A.; Andrejco, M. J.; DeLong, K. W.; Stegeman, G. I. Opt. Lett. 1989, 14, 1140–1142. (19) Hales, J. M.; Matichak, J.; Barlow, S.; Ohira, S.; Yesudas, K.; Brédas, J.-L.; Perry, J. W.; Marder, S. R. Science 2010, 327, 1485–1488. (20) Nanda, K. D.; Krylov, A. I. J. Chem. Phys. 2015, 142, 064118. (21) Olsen, J.; Jørgensen, P. J. Chem. Phys. 1985, 82, 3235–3264. (22) Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Chem. Rev. 2012, 112, 543–631. (23) Norman, P.; Bishop, D. M.; Jensen, H. J. A.; Oddershede, J. J. Chem. Phys. 2001, 115, 10323–10334. (24) Norman, P.; Bishop, D. M.; Jensen, H. J. A.; Oddershede, J. J. Chem. Phys. 2005, 123, 194103. (25) Macak, P.; Luo, Y.; ˚ Agren, H. Chem. Phys. Lett. 2000, 330, 447. (26) Kamarchik, E.; Krylov, A. I. J. Phys. Chem. Lett. 2011, 2, 488–492. 26


Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(27) Silverstein, D. W.; Jensen, L. J. Chem. Phys. 2012, 136, 064110. (28) Silverstein, D. W.; Jensen, L. J. Chem. Phys. 2012, 136, 064111. (29) Kristensen, K.; Kauczor, J.; Kjærgaard, T.; Jørgensen, P. J. Chem. Phys. 2009, 131, 044112. (30) Kristensen, K.; Kauczor, J.; Thorvaldsen, A. J.; Jørgensen, P.; Kjærgaard, T.; Rizzo, A. J. Chem. Phys. 2011, 134, 214104. (31) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. J. Comput. Chem. 2001, 22, 931–967. (32) Fonseca Guerra, C.; Snijders, J. G.; te Velde, G.; Baerends, E. J. Theor. Chem. Acc. 1998, 99, 391–403. (33) Baerends, E.; Ziegler, T.; Autschbach, J.; Bashford, D.; Bérces, A.; Bickelhaupt, F.; Bo, C.; Boerrigter, P.; Cavallo, L.; Chong, D.; Deng, L.; Dickson, R.; Ellis, D.; van Faassen, M.; Fan, L.; Fischer, T.; Guerra, C. F.; Franchini, M.; Ghysels, A.; Giammona, A.; van Gisbergen, S.; Götz, A.; Groeneveld, J.; Gritsenko, O.; Gr¨ uning, M.; Gusarov, S.; Harris, F.; van den Hoek, P.; Jacob, C.; Jacobsen, H.; Jensen, L.; Kaminski, J.; van Kessel, G.; Kootstra, F.; Kovalenko, A.; Krykunov, M.; van Lenthe, E.; McCormack, D.; Michalak, A.; Mitoraj, M.; Morton, S.; Neugebauer, J.; Nicu, V.; Noodleman, L.; Osinga, V.; Patchkovskii, S.; Pavanello, M.; Philipsen, P.; Post, D.; Pye, C.; Ravenek, W.; Rodriguez, J.; Ros, P.; Schipper, P.; Schreckenbach, G.; Seldenthuis, J.; Seth, M.; Snijders, J.; Solà, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T.; van Wezenbeek, E.; Wiesenekker, G.; Wolff, S.; Woo, T.; Yakovlev, A. Amsterdam Density Functional. 2013; http://www.scm.com. (34) Jensen, L.; Autschbach, J.; Schatz, G. C. J. Chem. Phys. 2005, 122, 224115.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 31

(35) Hu, Z.; Autschbach, J.; Jensen, L. J. Chem. Phys. 2014, 141, 124305. (36) Ye, A.; Autschbach, J. J. Chem. Phys. 2006, 125, 234101. (37) Bishop, D. M.; Kirtman, B.; Kurtz, H. A.; Rice, J. E. J. Chem. Phys. 1993, 98, 8024– 8030. (38) Bishop, D. M.; Pipin, J.; Kirtman, B. J. Chem. Phys. 1995, 102, 6778–6786. (39) Ingamells, V. E.; Papadopoulos, M. G.; Raptis, S. G. Chem. Phys. Lett. 1999, 307, 484 – 492. (40) Torruellas, W. E.; Zanoni, R.; Stegeman, G. I.; Möhlmann, G. R.; Erdhuisen, E. W. P.; Horsthuis, W. H. G. J. Chem. Phys. 1991, 94, 6851–6856. (41) Beljonne, D.; Brédas, J. L.; Cha, M.; Torruellas, W. E.; Stegeman, G. I.; Hofstraat, J. W.; Horsthuis, W. H. G.; Mohlmann, G. R. J. Chem. Phys. 1995, 103, 7834–7843. (42) Antonov, L.; Kamada, K.; Ohta, K.; Kamounah, F. S. Phys. Chem. Chem. Phys. 2003, 5, 1193–1197. (43) Day, P. N.; Nguyen, K. A.; Pachter, R. J. Phys. Chem. B 2005, 109, 1803–1814. (44) Day, P. N.; Nguyen, K. A.; Pachter, R. J. Chem. Phys. 2006, 125, 094103. (45) Boyd, R. W. Nonlinear Optics; Academic Press: San Diego, 1992. (46) Karna, S. P.; Dupuis, M. J. Comput. Chem. 1991, 12, 487–504. (47) van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. J. Chem. Phys. 1998, 109, 10644–10656, journal article. (48) Autschbach, J. ChemPhysChem 2011, 12, 3224–3235. (49) Banerjee, A.; Harbola, M. Eur. Phys. J. D 1999, 5, 201–206. 28


Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


´ (50) Angy´ an, J. G. J. Math. Chem. 2009, 46, 1–14. (51) Jansik, B.; Salek, P.; Jonsson, D.; Vahtras, O.; ˚ Agren, H. J. Chem. Phys. 2005, 122, 054107. (52) Orr, B.; Ward, J. Mol. Phys. 1971, 20, 513–526. (53) Willetts, A.; Rice, J. E.; Burland, D. M.; Shelton, D. P. J. Chem. Phys. 1992, 97, 7590–7599. (54) Pérez-Moreno, J.; Kuzyk, M. G. J. Chem. Phys. 2005, 123, 194101. (55) Pérez-Moreno, J.; Clays, K.; Kuzyk, M. G. J. Chem. Phys. 2008, 128, 084109. (56) McClain, W. M. J. Chem. Phys. 1971, 55, 2789. (57) Beerepoot, M. T. P.; Friese, D. H.; List, N. H.; Kongsted, J.; Ruud, K. Phys. Chem. Chem. Phys. 2015, 17, 19306–19314. (58) Göppert-Mayer, M. Ann. Phys. 1931, 401, 273. (59) Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100. (60) Perdew, J. P. Phys. Rev. B 1986, 33, 8822–8824. (61) ”DALTON, a molecular electronic structure program, Release Dalton2011 (2011), see http://daltonprogram.org/”. (62) Schipper, P. R. T.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. J. Chem. Phys. 2000, 112, 1344–1352. (63) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. J. Chem. Phys. 2003, 119, 12998–13006. (64) Shelton, D. P.; Rice, J. E. Chem. Rev. 1994, 94, 3–29. (65) Jensen, L.; van Duijnen, P. T. J. Chem. Phys. 2005, 123, 074307. 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(66) Bishop, D. M. Rev. Mod. Phys. 1990, 62, 343–374. (67) Hättig, C.; Christiansen, O.; Jørgensen, P. Chem. Phys. Lett. 1998, 282, 139–146. (68) Dalgaard, E. Phys. Rev. A 1982, 26, 42–52. (69) Li, Z.; Suo, B.; Liu, W. J. Chem. Phys. 2014, 141, 244105. (70) Ou, Q.; Bellchambers, G. D.; Furche, F.; Subotnik, J. E. J. Chem. Phys. 2015, 142, 064114. (71) Subotnik, J. E.; Alguire, E. C.; Ou, Q.; Landry, B. R.; Fatehi, S. Acc. Chem. Res. 2015, 48, 1340–1350. (72) Eriksen, J. J.; Sauer, S. P.; Mikkelsen, K. V.; Christiansen, O.; Jensen, H. J. A.; Kongsted, J. Mol. Phys. 2013, 111, 1235–1248.

30


Page 30 of 31

Page 31 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

31


Simulating Third-Order Nonlinear Optical ... - ACS Publications

Recommend Documents