Open-Ended Recursive Approach for the Calculation of Multiphoton Absorption Matrix Elements.

We present an implementation of single residues for response functions to arbitrary order using a recursive approach. Explicit expressions in terms of density-matrix-based response theory for the single residues of the linear, quadratic, cubic, and quartic response functions are also presented. These residues correspond to one-, two-, three- and four-photon transition matrix elements. The newly developed code is used to calculate the one-, two-, three- and four-photon absorption cross sections of para-nitroaniline and para-nitroaminostilbene, making this the first treatment of four-photon absorption in the framework of response theory. We find that the calculated multiphoton absorption cross sections are not very sensitive to the size of the basis set as long as a reasonably large basis set with diffuse functions is used. The choice of exchange-correlation functional, however, significantly affects the calculated cross sections of both charge-transfer transitions and other transitions, in particular, for the larger para-nitroaminostilbene molecule. We therefore recommend the use of a range-separated exchange-correlation functional in combination with the augmented correlation-consistent double-ζ basis set aug-cc-pVDZ for the calculation of multiphoton absorption properties.

Introduction

The calculation of transition properties between ground and excited states has long been a challenge for quantum chemistry. In 1985, Olsen and Jørgensen showed that the transition moments between the ground and excited states of a molecular system can be obtained from the single residues of its ground-state response functions.[1] They provided expressions for the calculation of one-, two-, and three-photon transition matrix elements, corresponding to the description of one-, two-, and three-photon absorption, respectively.

Two-photon absorption (TPA), in particular, has received a lot of attention. Despite the early prediction by Göppert-Mayer in 1931,[2] a measurement of this nonlinear optical effect was not reported until 1961.[3] The long time between the theoretical prediction and the first measurement of this phenomenon can be explained by the need for a high laser intensity to match the corresponding small absorption cross section, which depends quadratically on the intensity of the incident light. Several interesting applications for TPA, in particular, and multiphoton absorption (MPA), in general, have been proposed since their first experimental realization, including 3D data storage, multiphoton microscopy, photodynamic cancer therapy, and drug delivery.[4−7] These applications require the design of materials with a high MPA cross section.[7,8]

Higher-order MPA properties suffer from even smaller absorption cross sections than TPA. In general, the j-photon absorption strength depends on the radiation intensity to j-th order if j is the number of photons to be absorbed.[9] Nevertheless, modern high-intensity laser pulse techniques can achieve the intensity necessary for the observation of these higher-order absorption properties. Frequency up-conversion of infrared semiconductor lasers is an important field in which multiphoton absorption is applied.[10,11] This technique results in a source of coherent light at short wavelengths that is pumped by laser diodes emitting in the infrared region. In recent years, several chromophores have been synthesized that show stimulated emission after multiphoton absorption processes. In these studies, three-,[11,12] four-,[13] and even five-photon absorption have been measured.[14,15]

Rational design of molecules with favorable absorption properties is often aided by computational studies. Therefore, a lot of effort has been put in enabling the calculation of such properties in quantum-chemical program packages. TPA has been implemented for self-consistent field (SCF)[16,17] and multiconfigurational SCF (MCSCF) methods[16] as well as for several coupled cluster models.[18−21] Implementations of three-photon absorption (3PA), however, are few.[22−24] Implementations of higher-order absorption properties have not yet been realized. However, general expressions for four-photon absorption cross sections have been published by Andrews and Ghoul[9] and for coupled cluster theory, a general scheme for the derivation of j-photon absorption strengths has been presented by Hättig, Christiansen, and Jørgensen.[25]

A reason for the lack of implementations of four-photon absorption (or higher-order multiphoton absorption) is the complexity of the corresponding expressions. Olsen and Jørgensen have shown in their 1985 study that j-photon absorption matrix elements can be obtained from the residue of the response function of at least order j + 1.[1] Thus, the matrix elements for one-photon absorption are obtained from the residue of the linear response function and those for two-photon absorption from the quadratic response function, which correspond to the second and third derivatives of the energy with respect to external perturbations, respectively. Because linear response functions and their residues are of moderate complexity, implementations of one-photon absorption are common. Implementations of two-photon absorption are fewer, especially in the framework of response theory. The calculation of third-, fourth-, and fifth-order absorption matrix elements requires residues of the cubic, quartic, and quintic response functions, corresponding to the fourth-, fifth-, and sixth-order derivatives of the energy, respectively. The complexity of calculating such high-order properties increases rapidly with the order of the property, making a tailored implementation of properties at each higher order increasingly complicated.

To tackle the increasing complexity of high-order response functions, Ringholm, Jonsson, and Ruud recently presented an open-ended code for the calculation of response properties at the Hartree–Fock (HF) and density-functional theory (DFT) level. By its recursive nature, the code is capable of identifying the contributions to response functions to any order and assemble them.[26] The main aim of the present work is to present an implementation of single residues within this recursive scheme that will enable the calculation of absorption properties to arbitrary order and to present the first computational study of four-photon absorption.

The rest of this paper is organized as follows: In Section 2, we present the theory of single residues of response functions at the HF and DFT level; in Section 3, we present the integration of our findings into the recursive scheme of Ringholm et al.[26] In Section 4, the formation of (isotropic) transition cross sections from multiphoton transition matrix elements is discussed and universal schemes for the handling of multiphoton absorption to arbitrary order will be presented. In Section 5, we give the computational details, and in Section 6, we present the first comprehensive computational study of four-photon absorption (4PA), including a detailed study of the dependence of the 4PA cross section on the choice of basis set and exchange–correlation functional. We have calculated the 4PA cross sections for the para-nitroaniline (PNA) and para-nitroaminostilbene (PNAS) molecules. These are two rather small push–pull substituted aromatic chromophores that are frequently used in theoretical studies of multiphoton absorption as model compounds for systems with excited states displaying charge-transfer character.[21,24,27] Additionally, we have calculated OPA, TPA, and 3PA to compare the MPA behavior of the calculated states. In Section 7, we give some concluding remarks.

Theory

In this section, we first give a brief overview of response theory in a density-matrix-based formulation (Section 2.1), before we discuss how the single residues of these response functions are formed (Section 2.2).

Response Functions in Density-Based Response Theory

A general theory for response functions in a density-based formulation was presented by Thorvaldsen et al. in 2008,[28] and we therefore limit ourselves here to a brief overview. We also limit our discussion to perturbations that do not influence the basis set (e.g., an applied electric dipole field), without showing in detail the more comprehensive form these expressions take in the general case. A detailed discussion of single residues of response functions involving perturbation-dependent basis sets is planned for future work.

Time-Dependent Quasienergy

In response theory, the time dependence of a perturbed system is expanded in a Fourier sum of periodic perturbations and described by a quasienergy expression. The quasienergy plays a role that is analogous to the energy in the static case and is in the time-independent case reduced to the regular electronic energy. The explicit time dependence of the quasienergy is commonly removed by time averaging, yielding the time-averaged quasienergy, which is a key step in the determination of the response functions. A comprehensive treatment of the quasienergy approach has been given by Christiansen, Jørgensen, and Hättig.[29] The application of this approach to density-matrix-based Kohn–Sham theory was presented by Thorvaldsen et al.[28]

It is not possible to express the time-dependent Schrödinger equation (TDSE) in terms of the density matrix only, since the TDSE is not symmetric with respect to time differentiation. For this reason, the starting point for density-matrix-based response theory is the quasienergy derivative Lagrangian defined with respect to an applied perturbation a. The first derivative of the time-averaged quasi-energy iswith Ẽ0, being the energy differentiated to first order with respect to the perturbation a and to zeroth order with respect to the density D̃. The tilde denotes both time dependence and evaluation at a general field strength.

Quasienergy Derivative Lagrangian

Response functions in density-based response theory are obtained from the derivatives of the general quasienergy derivative Lagrangianwhere the Lagrangian multipliers λ̃ are defined aswhere S̃ is the overlap matrix and we have introduced the short-hand notation A⊖

The perturbed densities D are discussed in Section 2.1.3. The matrix Ỹ represents a constraint defined by the time-dependent SCF (TDSCF) equationwhere F̃ is the Fock matrixwith h̃ being the one-electron kinetic energy and nuclear attraction integrals and Ṽ being the integrals of a one-electron perturbation operator. G̃γ(D) represents the contribution of the two-electron repulsion integrals contracted with the unperturbed density D and with the exchange contribution scaled by the factor γ, whereas F̃xc is the exchange–correlation contribution to the Fock matrix, which only enters at the DFT level of theory. We note that Hartree–Fock theory is recovered by setting γ = 1 and removing F̃xc, whereas pure DFT is recovered by choosing γ = 0.

The second set of Lagrangian multipliers ζ̃ is defined aswhere the notation A⊕ in a similar manner as eq 4 is defined as

The perturbed Fock matrix F will be defined in eq 22. The matrix Z̃ is given byand represents the so-called idempotency condition for the density and is the second constraint that must be satisfied.

Response functions can now be formulated as perturbation-strength derivatives of eq 2 evaluated at zero perturbation strengths. These expressions depend on perturbed Fock and density matrices F and D. The Wigner rules of perturbation theory can be applied to the perturbation parameters[30] to evaluate the response functions in a flexible and efficient way.

Because the perturbation designated as perturbation a is already included in the quasienergy derivative Lagrangian, the expressions for the quasienergy Lagrangian derivatives will not be symmetric with respect to perturbation a and the other perturbations. Letting a response property be characterized by the perturbation tuple abc..., where a,b,c... are the individual perturbations, when specifying the highest order needed for the different perturbation parameters, two numbers must be provided: one for perturbation tuples involving perturbation a (this maximum order will be called k in the following) and one for perturbations not involving a (called n in the following). The integers k and n can be chosen freely as long as they match the conditions k + n = N – 1, with N being the total number of perturbations, and k ∈ [0,(N + 1)/2], where the fraction for the maximum of the interval is rounded down for even N.

Perturbed Densities

The perturbed parameters in density-matrix-based response theory are expressed in terms of perturbed densities D, where the superscript denotes differentiation with respect to the perturbation strengths of an arbitrary collection K of perturbations. The collection of perturbations K can either consist of one single perturbation (in our case always a spatial component of the electric field) or of several perturbations (several components of the electric field). In the formulation of Thorvaldsen et al., the perturbed density is defined as[28]

The perturbed parameters X that are used to set up the perturbed densities are solutions of linear sets of response equations of the formwhere E[2] and S[2] are the generalized Hessian and metric matrices, respectively. For details about these quantities, we refer to refs (28) and (31). The frequency ω is defined as the sum of the frequencies of the individual perturbations forming the set of perturbations K. Throughout this work, we refer to this set of perturbations as a “perturbation tuple” in accordance with the nomenclature used in ref (26). The matrix MRHS is the so-called right-hand side vector, which for the one-electron perturbations considered here is obtained as the derivative of the TDSCF equation eq 5. Explicit expressions for the right-hand side vectors will be given in the next section.

Formulation of Residues

The response function has a pole whenever the frequency ω approaches an excitation energy, at which point the corresponding response property diverges. The calculation of excitation energies, which will not be discussed in detail here, is thus based on the determination of the poles of the linear response function. The residues of the response function at the pole correspond to the transition properties between the ground state and the corresponding excited state.[1]

In general, the residue of a response function is obtained by multiplying the response function with (ω – ω), where ω is the excitation energy from the ground state to the excited state p, and forming the limit of the resulting product with respect to the frequency of perturbation K approaching ω. For limω, the difference (ω – ω) is zero, which means that all terms depending on (ω – ω) in the product vanish. Only the terms in the response function containing (ω – ω)−1 remain as (ω – ω) then cancels to unity.

The theory for constructing residues from the response functions in the density-matrix-based framework has been presented by Thorvaldsen et al.,[28] and in the following, we show this technique explicitly. The residue of eq 11 can be written

The denominator on the right-hand-side of eq 11 can, in the spectral representation, be written aswhere X represents the excitation eigenvector to an excited state q. Combining eqs 11, 12, and 13, we obtain for the residue of the perturbed parameters

Therefore, for the residue of X, the solution of the corresponding response equations (eq 11) reduces to the multiplication of X with a scalar that can be obtained from X† and MRHS.

Applying this technique to the whole response function, only those terms remain that depend on X or D, and the residues can therefore be obtained by applying the following procedure to the response function:[28]

Remove all terms that do not depend on D either explicitly or implicitly.

Replace D by D, where D is defined as

If the residue of a density of order k is formed, then this residue has to replace the corresponding perturbed densities in all expressions for perturbed densities of order > k.

In the following, this technique will be applied to the linear, quadratic, cubic, and quartic response functions.

Residues of the Linear Response Function

The linear response function corresponds to the second derivative of the energy with respect to an external perturbation. It describes molecular properties such as the polarizability, for which the perturbations are two electric dipole operators. The residue of the linear response function in the density-matrix-based framework has been reported by Thorvaldsen et al.[28] and is here repeated for completenessAs shown by Thorvaldsen et al., the residue in eq 16 is the product of the left and right one-photon transition matrix element, which are related through complex conjugation.

Residues of the Quadratic Response Function

The quadratic response function describes properties such as the first hyperpolarizability. From the residues of the quadratic response function, two-photon absorption matrix elements can be obtained.[1,16]

Making the rule choice k = n = 1, the response function is given bywith the intermediate quantitieswhere the index 1′ means that perturbed parameters are taken into account only up to first order in the corresponding quantity. All perturbed densities depend on the frequency of the corresponding perturbation.

Kjærgaard et al. have presented an expression for two-photon absorption matrix elements in AO-based response theory.[31] In the approach we use here, however, these matrix elements are obtained using the procedure we have described above, so that

Note that ω = ω due to the residue formation and that the expression in eq 23 is valid for both HF and DFT but is limited to one-electron perturbations. The expressions will therefore not contain any integral terms that are higher than first order in the integral derivatives. Terms depending on the density to higher than quadratic order can be nonzero only at the DFT level since only the exchange–correlation contribution has a dependence on the density beyond quadratic order. For HF, all terms that are higher than second order in the density vanish. For instance, the second contribution to the F1 intermediate (vide infra) is nonzero at the DFT level but vanishes at the HF level.

The quantities used in eq 23 are defined asThe residue of the perturbed density D is defined aswith the residue of the perturbed parameters defined aswhere the notation F̆ denotes all contributions to the Fock matrix that are independent of the perturbation parameters. The notation we use here differs in some aspects from the one used in refs (26) and (28) because we restrict ourselves to perturbations that do not affect the basis set. Formally, the expression in eq 23 is similar to the one by Kjærgaard et al.[31] as it only contains perturbation parameters up to first order. Nevertheless, an important difference between the two formulations is that the expression by Kjærgaard et al. is symmetric with respect to the perturbations involved, whereas in our case, the perturbation labeled as perturbation a plays a privileged role, and it is therefore not possible to symmetrize the expressions completely.

Another expression for the quadratic response function can be obtained by using the rule choice (k,n) = (0,2)[28]for which the residue can be formed with a frequency sum approaching the excitation energy as

Equation 33 contains one term less than eq 23, but it involves the residue of a doubly perturbed density not present in eq 23. Equation 33 is linear in D, as is the case for the linear response function and the first-order density in eq 16.

In the density-based formulation, it is rather straightforward to determine the residue of a doubly perturbed density. We will show this by applying the technique shown in eqs 28 and 29 to the residue of a doubly perturbed density, which, according to eq 15, is defined aswhere the residue of the perturbed parameters is

The contraction Tr(X†MRHS) gives a scalar and represents the right two-photon absorption matrix element. The left and the right absorption matrix elements are in general adjoint to each other in variational theory,[29] and therefore, only one of the matrix elements has to be determined.

The vector MRHS is the right-hand-side vector for the solution of the second-order response equation and can be obtained from the second derivative of the TDSCF equations.[28] For perturbations not affecting the basis set, MRHS is defined as

From eq 36 we also note that, as in the formulation by Kjærgaard et al., the second-order transition matrix elements can be determined using perturbed parameters to first order only, although the underlying expression was formulated without using the 2m + 1 rule. This illustrates that residue formation can reduce the level up to which response equations have to be solved. We note from eq 36 that this expression only depends on perturbed densities to first order, and hence, it can be calculated once the perturbation parameters to first order are determined. This is also the case for the residue formulation in eq 23.

Both eqs 23 and 33 can be assumed to be a product of a left and a right transition matrix element whose characteristics depend on the formulation. Equation 23 consists of a left second- and a right first-order transition matrix element, whereas for eq 33, it is the converse (left first- and right second-order transition matrix element).

The decomposition in transition matrix elements is due to the structure of the residues of the perturbed parameters and densities (eqs 14 and 15): the residue of the densities can be decomposed into a product of an excitation density and a scalar that is obtained as a contraction of the excitation eigenvector and the corresponding right-hand-side vector

Examining eqs 23 and 33, we note that all terms in these equations are linear in D or D, respectively. Since the perturbation dependence in these density residues depends only on multiplication with a scalar, it can be extracted from the equations, yieldingfor eq 23 andfor eq 33. The scalars that are extracted from the expressions are the right transition matrix elements, whereas the remainder of the expression is the left transition matrix element. The second-order transition matrix elements are written as

These matrix elements are related by complex conjugation in variational theory,[29] and therefore, only one of the matrix elements has to be determined. Since both approaches presented here are equivalent in their requirements of perturbed parameters and solution of linear equation systems, either expression can be chosen.

Residues of the Cubic Response Function

The cubic response function describes properties such as the second hyperpolarizability. Additionally, it is the lowest-order response function that allows the calculation of third-order transition properties from the corresponding residues. Furthermore, products of left and right second-order transition matrix elements can be obtained from the cubic response function, but because these matrix elements can also be obtained from the quadratic response function, there is no need for this. Three-photon absorption matrix elements can be obtained from the residues of the cubic response function.

In the following, we will discuss the extraction of the third-order transition matrix elements from the cubic response function in two different formulations.

Setting k = 1 and n = 2, the cubic response function can be written as

The intermediate quantities used here are defined in Section 1 of the Supporting Information.

With ω approaching an excitation energy, we obtain for the residuewith the intermediate quantities being defined in Section 2 of the Supporting Information.

The perturbed densities of the type D depend on two perturbations with the frequency of one of them approaching the excitation energy. At this level of residue formation, the influence of the residue formation on the right-hand side vector MRHS must be taken into account, in contrast to the expressions for the right-hand side vectors for the quadratic response function, eqs 30 and 36.

The response parameters X arewith the right-hand side vector

The evaluation of the residue of the density D is analogous to the procedure used for the quadratic response function (see eqs 28, 29, and 30), yielding a right transition matrix element of first order and a left transition matrix element of third order , defined as

The quantities F2, Y2′, and Z2′ are formed from F2, Y2′, and Z2′, respectively, by exchanging D and D with D and D, respectively.

The notation D means that the perturbed density depends on a perturbation whose frequency approaches the excitation energy at a lower perturbation level, in this case perturbation d, which contributed to the original perturbed density D from which D has been formed. The residue D in eq 44 can be considered a product of D and the right first-order transition matrix elementEquation 51 is a linear equation system. Its right-hand-side vector MRHS can be derived from eq 47 by factorizing out the scalar X†MRHS because the right-hand-side vector for the corresponding perturbed parameters eq 46 is linear in X and thus according to eq 37 corresponds to D.

As D and D require the solution of a linear equation system, they must be treated as doubly perturbed densities. From a computational point of view, thus depends on first- and second-order perturbed densities.

Using k = 0 and n = 3, corresponding to the n + 1 rule, the cubic response function is

Forming the residue with ω + ω + ω → ω approaching the excitation energy, we obtain

We note that eq 44 is linear in D, which, in analogy to the approach used for the second-order transition matrix elements, can be expressed as

For the residue of the third-order response parameters X, we find thatwhich enables us to determine the right third-order transition matrix element aswith MRHS being defined asWe note that MRHS depends only on first- and second-order perturbed matrices. Therefore, also in this case, we have to solve the first- and second-order response equations to get the third-order transition matrix elements. There is therefore no formal difference in the computational requirements for and , but there is nevertheless one difference in the computational requirements for the two approaches, because only the excitation eigenvectors and MRHS are needed for , making this easier to implement than .

Residues of the Quartic Response Function

The quartic response function describes properties such as the third hyperpolarizability. The residue of the quartic response function enables us to treat absorption matrix elements of fourth order. Since the expressions for the quartic response function are prohibitively long independently of which rule is used, they are not given explicitly here, and we restrict ourselves to defining the residues we want to discuss. For the formulation using k = 2 and n = 2, the residue can be written as a product of a left four- and a right one-photon-absorption matrix element aswhere 2′ means that the corresponding quantity depends only on derivatives up to second order. The intermediates from eq 59 are defined in Section 3 of the Supporting Information.

The complexity of eq 59 draws attention to the importance of a recursive implementation scheme for such high-order properties, and we return to this point in the next section.

Examining eq 59, we note that the four-photon transition matrix elementcan be obtained using perturbed densities to first and second order, or from intermediates to the same order as the third-order transition matrix elements.

Using the m + 1 rule, the residue of the quartic response function for electric field perturbations can be written

Equation 61 is again linear in D, which is defined asX can be writtenwith the right-hand-side vector MRHS defined asUsing the m+1 rule, the four-photon transition matrix elements require the calculation of perturbed parameters up to third order (for MRHS). This means that a formulation using the 2m+1 rule is more efficient since it only requires second-order parameters. Nevertheless, the m+1 formulation has the benefit of an easier implementation and can therefore serve for testing and debugging the more efficient 2m+1 formulation.

Discussion of the Residue Expressions

For the formulation of the residues, we can formulate the following general rules, which are valid to any order:

The residues can be decomposed into a left and a right transition matrix element. In the examples discussed here, one of the two matrix elements forming each residue is a first-order transition matrix element. However, this is not a requirement.

Left or right jth-order transition matrix elements are obtained from the residue of the jth-order response function.

Up to the cubic response function, the computational requirements do not depend on the choice of 2m+1 or m+1 rule because the same order of the perturbed density matrices is needed. For higher-order response functions the choice of the rule has an impact on the computational requirements of the calculations.

The requirements for perturbed parameters (i.e. orders to which response equations have to be solved) are listed in Table 1 using the most computationally efficient 2m + 1 rule for the cases discussed above as well as for the fifth-order transition matrix elements not shown here. As can be seen from Table 1, the computational requirements for jth-order transition matrix element calculations follow the 2m + 1 rule of regular response functions.

Table 1

Response Function and Parameter Requirements for Transition Matrix Elements of Different Order Using the 2m + 1 Rule

order of transition matrix element	order of response function	order of params.
1	2	0
2	3	1
3	4	2
4	5	2
5	6	3

Implementation

Section 2 has illustrated the rapid increase in complexity of the response functions and their residues with increasing order of the perturbations. To tackle the otherwise similarly rapidly increasing programming effort, Ringholm et al. presented an open-ended recursive scheme for calculating response functions to arbitrary order.[26,28,32]

As seen in the previous section, the residues can be obtained from the response functions by removing all terms that do not depend on the frequency that matches the excitation energy and by replacing the corresponding perturbed density with its residue in the remaining terms. As we will show, this makes the modification of the recursive scheme by Ringholm et al. straightforward. For a detailed description of the original open-ended scheme, we refer to ref (26).

Calculation and Handling of Excitation Eigenvectors

The first requirement for a residue calculation is the calculation of the excitation energies and the excited-state eigenvectors X, which can be determined as a solution of the generalized eigenvalue problem (see, e.g., ref (31).)obtainable from, for example, the response solver[33] in the DALTON program.[34,35] The eigenvectors X determined in this procedure can be handled the same way as other perturbation parameters.

The calculation of the excitation eigenvectors is performed in a separate step before the residues are calculated. After the determination of X, the excitation densities D are formed using eq 39. Both X and D are stored in linked lists,[36] and these linked lists are used by the code for storing and handling perturbed intermediates in the response function calculations.

Modifications to the Response Function Calculation

The response function calculation is based on a set of five algorithms used to identify the quantities that need to be calculated and for their proper assembly into the final result. Only two of these algorithms (Algorithms 2 and 3) need to be modified to enable open-ended calculations of residues of the response function. In addition, a new algorithm (Algorithm 6) is needed and will be described. In order to provide a proper reference frame for the modifications that we will describe, we recapitulate the main purpose and functionality of Algorithms 2 and 3. For a complete description of these and the other algorithms, we refer to ref (26). Algorithm 2 manages the calculation of the perturbed F and D matrices that have been identified as necessary in order to calculate the contributions to the response function. In Algorithm 3, relevant contributions to the perturbed F or energy-type contributions to the response tensor are identified according to which property is to be calculated and the choice of (k,n) rule made. The necessary modifications to these two algorithms will be described in the following.

Modifications to Algorithm 2

Algorithm 2 manages the calculation of the perturbed intermediates and therefore needs some way of recognizing whether response equations have to be solved (eq 11) for some given set of perturbations, or whether instead the residue of the perturbation parameter is to be formed (eq 14). This is accomplished by comparing the sum of the frequencies in the perturbation tuple considered with the excitation energy. If these two numbers are sufficiently close to each other, only the right-hand-side vector MRHS is created and contracted with X, forming the corresponding right transition matrix element—that is, the second term on the right-hand side of, for example, eqs 37 and 38.

The perturbed parameters X are only used to calculate the perturbed density matrix. In residue calculations, this is also done to create the perturbed densities from X according to eq 15. The details of the modified algorithm are summarized in Figure 1.

Figure 1

Pseudocode of Algorithm 2 for the calculation of F and D in residue mode. Based on Algorithm 2 from ref (26).

Modifications to Algorithm 3

Algorithm 3 is used to calculate contributions to perturbed Fock matrices, the perturbed energy-type contributions and, with minor modifications, to the other contributions to the response tensor. It identifies the relevant contributions to each of these quantities, consisting of perturbed or unperturbed Fock matrices or energy expressions contracted with various tuples of perturbed or unperturbed density matrices for the former two quantities, and various other terms for the latter.

In residue calculations, Algorithm 3 has to be supplemented by a structure that recognizes whether the different contributions are zero or not due to the formation of a residue. Consider for instance the perturbed Fock matixwhich is needed for the calculation of three-photon transition moments (see Supporting Information, Section 1). All its terms contribute to the cubic response function in eq 43. Forming the residue, we have to determine which of these terms contribute to the residue—that is, whether the perturbed density depends on a frequency that approaches the excitation energy ω. From the contributions in eq 66, for example, the first ((D)) and the fourth ((D,D)) enter the residue as (D) and (D,D(), whereas the third term ((D)) does not contain an ω-dependent density and therefore vanishes in the residue formation.

As seen from this example, either the surviving term can depend exactly on the frequency approaching the excitation energy as for D or the corresponding frequency can be contained in the frequencies that the corresponding intermediate depends on, as was the case for D. In order to identify both these cases, it is not sufficient to check only whether the frequencies add up to the excitation energy, as done in Algorithm 2. Instead, an additional recursive algorithm denoted Algorithm 6 is needed.

Algorithm 6

Since the residue code presented here is intended to be as universal as possible, the identification of the correct contributions to the residues is a computational challenge. For a perturbation tuple with respect to which a given perturbed density matrix for a residue is to be formed, the sum of some or all of the associated frequencies of the perturbation tuple can coincide with the excitation energy. Determining the corresponding contribution is therefore made more difficult by the possibility that it is not necessarily only one frequency approaching the excitation energy but also combinations of t frequencies, with t being an integer between 1 and N, with N being the number of perturbations considered.

Identifying these different appearances of the excitation energy can be achieved by a recursive function that examines all frequencies and all sums of frequencies associated with any subset of the perturbation tuple considered. A pseudocode of the algorithm is shown in Figure 2.

Figure 2

Pseudocode of Algorithm 6 to recognize the excitation energy from the frequencies of a perturbed intermediate (b,t,w,j), the result is returned in the logical variable recognized..

At the first invocation of this procedure, the number of frequencies from which the excitation energy in the frequency tuple is formed (t) has to be known. Another variable to be provided as an argument to this algorithm is the perturbation tuple b, which contains all information about the perturbation considered.

The algorithm forms all possible t-element sums from the perturbation frequencies and compares them with the excitation energy. This is achieved by a recursive structure that calls itself t times for every frequency between the first and the (N – t)th position in the list, and it can therefore form a t-fold sum of frequencies on the earlier positions in the list.

Postprocessing of the Transition Matrix Elements

The response theory described in Section 2 yields multiphoton transition matrix elements. Several steps are needed to extract observable quantities from them. The matrix elements need to be multiplied to form transition strength tensors and these tensors need to be rotationally averaged if we consider isotropic samples. The rotationally averaged tensors can then be converted into multiphoton absorption cross sections. In the following, we discuss this postprocessing of the transition matrix elements step by step.

Formation of the Transition Strength Tensors and Rotational Averaging

It is important to note that there is a difference between transition matrix elements as described in Section 2 and the transition strength tensors or absorption cross sections that can be compared to an experimental observation. In contrast to the transition matrix elements that can be obtained from the single residues of the corresponding response functions as listed in Table 1, the j-photon transition strength tensor depends on a single residue of the 2jth-order response function. For example, the two-photon transition strength tensor depends on the residue of the cubic response function although the transition matrix elements from which it can be constructed are obtained from the quadratic response function.[1] Hence, the tensors that have to be handled in the postprocessing of the transition matrix elements are of higher rank. The 4PA transition strength tensor is of rank 8. In general, these transition strength tensors are obtained as tensor products of outer form of two transition matrix elements. This step is usually not conducted explicitly for the treatment of isotropic samples but combined with the rotational averaging of the transition strength tensor.

The rotational averaging of tensors representing molecular response functions or their residues has been investigated intensively, especially in the 1970s and early 1980s. The first expression for rotational averaging of the TPA transition strength was formulated by Monson and McClain in 1970.[37] Andrews and co-workers have presented formulas for the rotational averaging of tensors up to eighth rank.[38,39] More work on this has been presented by Wagnière.[40] In order to be able to treat multiphoton absorption to arbitrary order, we recently presented universal equations for rotational averaging of tensors of even rank for linearly polarized light.[41]

The rotational averaging is based on the contraction of the transition matrix elements to one isotropic value, which can be correlated to an observable. For one-, two-, three-, and four-photon absorption, these expressions are[38,41]where and are the left- and right transition matrix elements, respectively.

Formation of Absorption Cross Sections

The rotationally averaged transition strengths can now be used to determine multiphoton absorption cross sections. For TPA, this is a well-known procedure that has been described first by Peticolas in 1967.[42] Two-photon absorption cross sections are usually given in Göppert-Mayer units (GM) with

Units for higher-order multiphoton absorption cross sections can be defined in analogy to this expression. In Table 2, we give an overview of the units commonly used in experimental work.[12−15] In order to derive the universal expression for the j-photon absorption cross section, we start with an expression for the multiphoton absorption probability from the ground state to a final state p, which is formulated analogously to eqs 10 and 28 in ref (42), correcting at the same time a couple of misprints in the original paper:In this expression, ω is the circular frequency of photon i, e is the elementary charge, m is the electron mass, and ρ(E) is the density of states of photons with energy E. n is the number of photons of energy i, and the fraction nρ(E)/V is the number of photons per volume V in the energy interval between E and E + dE where d determines the width of the interval.

Table 2

Units for Multiphoton Absorption Cross Sections, Given in the cgs Unit System

	unit
TPA	(cm4·s)/photon
3PA	(cm6·s2)/photon2
4PA	(cm8·s3)/photon3
jPA	(cm2j·sj–1)/photonj–1

The formulation in ref (42) is based on a transition strength tensor ⟨δ⟩, which is written in terms of the momentum operator p. The brackets ⟨⟩ represent the rotational average as described in the previous subsection. In the following, we derive an expression using the position operator from this approach following the lines of Peticolas.[42] To do this, we first simplify the one-photon transition moment using the dipole approximation[42]

Using the relation between the length and velocity gaugeswhere r is the position operator, we can transform the transition strength tensor to a quantity that is expressed in terms of the position operatorwhere δ denotes the square of the j-photon transition matrix elements discussed in the previous section.

Combining eqs 72 and 75, we get

Equation 76 can be further simplified by introducing the photon fluxand the energy fluxcorresponding to photons with frequency ω with c0 being the speed of light.[42] The product of these two quantities can be used to substitute the nρ(E)/V-terms in eq 72 according to

Using eq 79 in eq 76 and expanding with respect to ℏ[∏ω], we obtainwhich is a generalized and slightly modified version of eq 17 in ref (42). It now enables us to determine the transition probabilities for one-, two-, three-, and four-photon absorption according towhere eq 82 resembles eq 17 in ref (42). All expressions show transition probabilities of onephotonin the presence of the j – 1 other photons. This becomes clear if we consequently interpret the term ℏω not only as energy but as energy per photon. All these expressions can then be interpreted as having the unit photon/s.

In our expressions, as well as in the expressions in ref (42), the j photons are treated differently, as the transition probability is formulated for one photon (which is formally the photon with index j) in the presence of the other photons. Equations 82–84 can be transformed to expressions for the absorption cross sections following ref (42) by dividing by the photon flux F for all involved “types” of photons. To substitute the energy flux corresponding to photon j, which is always present in the expression, we use the equality

Furthermore, we have to take into account the broadening of the absorption band.[42] Experimentally, the absorption cross sections do not occur at discrete energies but rather in absorption bands. We do not discuss the origin of this broadening here; referring instead to the literature,[42,43] we only note that it is customary to represent this broadening by a line shape function g(∑ω) that satisfies the conditionWe thus note that this conditions leads to g(∑ω) having the dimension of a reciprocal frequency.

Expanding eq 86 with ω, we note that the transition probabilities in eqs 80–84 can be interpreted as integrals over a whole spectrum. The transition probability at a given frequency is therefore ωg(∑ω)dω. This is the starting point for the determination of the absorption cross section, which is also formulated for every frequency.

Using eq 85, we can now formulate a general expression for the j-photon absorption cross section asand the absorption cross sections for two-, three-, and four-photon absorption thus become

Equations 89–91 contain a large number of quantities that are constant. Moreover, we note that ⟨δ⟩ in eq 75 has a systematic dimensionality ofwhich enables us to gather most of the components of the cross section in one prefactor for the conversion from the rotationally averaged transition strength tensor to the cross section. These prefactors consist of powers of 2π and the conversion factors of the speed of light (reciprocal value of the fine structure constant, 0.007297353), length (Bohr radius, 5.2918 × 10–9 cm/au) and time (2.42 × 10–17 s/au) from atomic units to the centimeter–gram–second (cgs) system. Therefore, we get the jPA cross sections in the units in Table 2 as

Note that the prefactors differ in the literature as many authors combine them with the prefactor of the rotational averaging (see eqs 67–70). When scaled with 1050, eq 93 yields the TPA cross section in the well-established Göppert-Mayer units. In analogy to this, we here scale the units for 3PA and 4PA in eqs 94 and 95 with 1080 and 10110, respectively.

Computational Details

Molecular Structures

Geometry optimizations have been performed in vacuo for the para-nitroaniline (PNA) and para-nitroaminostilbene (PNAS) molecules (see Figure 3 for chemical structures) with the B3LYP[44] functional using the cc-pVQZ basis set.[45] For PNAS, this leads to a structure with all atoms in one plane. For PNA, the aromatic ring and the nitro group are in a single plane, whereas the hydrogens of the amino group are on the same side slightly out-of-plane. All geometry optimizations were performed using Gaussian[46] with default convergence criteria.

Figure 3

Chemical structures of the molecules under consideration.

Multiphoton Calculations

All calculations of multiphoton transition moments have been performed using our implementation of single residues in an open-ended response code,[26] which is used as an independent module in the DALTON program.[34,35] We have used the correlation-consistent polarized basis sets cc-pVDZ and cc-pVTZ,[45] as well as the more diffuse aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets.[47] In our calculations of multiphoton absorption cross sections, the HF method as well as the BLYP,[48] B3LYP,[44] and CAM-B3LYP[49] density functionals have been used. The three density functionals represent the generalized gradient approximation (GGA), hybrid, and range-separated hybrid density functional classes, respectively. These density functional classes differ by their treatment of the exchange contribution. GGA functionals only contain approximate exchange, whereas hybrid functionals have a mixture of approximate and Hartree–Fock exchange. In range-separated hybrid functionals, the contribution of exact Hartree–Fock exchange is variable depending on the distance, which gives these functionals a larger amount of flexibility in the description of charge-transfer excitations. It has been shown that the CAM-B3LYP functional gives a description of two-photon absorption comparable to the approximate coupled cluster model CC2.[21]

For the line shape function g(∑ω), a Lorentzian function with a width constant of 0.1 eV centered at the excitation energy and normalized to unity in the frequency domain was used, as frequently done in the literature.[50] As the Lorentzian is evaluated only at its center (the line shape function discussed in Section 4 is evaluated at the sum of all photon energies which is the same as the excitation energy), this reduces to a multiplication by a constant factor.

Applications

In this section, we discuss the results of our calculations. The PNA and PNAS molecules have been chosen for a first study to explore the dependence of the multiphoton absorption cross sections on the choice of basis set and exchange-correlation functional. We will discuss all one- and multiphoton absorption cross sections, with the 4PA results being the first calculations of this property in the literature.

para-Nitroaniline

PNA is a popular molecule for testing new computational models for evaluating optical properties because of its small size and its strong push–pull character, especially when it comes to solvatochromism and nonlinear optical properties.[51−53]

The characterization of the different excited states is crucial to understand the corresponding multiphoton processes and to evaluate the computational method properly. We will therefore first consider the excitation energies and the character of the different excited states. For all systems and computational levels, the five lowest excited states have been investigated.

We present the excitation energies and the character of the excited states based on an analysis of the dominant orbital transitions in Table 3. We see that there are differences between the results obtained with different exchange-correlation functionals. We note that the first five states in CAM-B3LYP and HF are identical in character, albeit with some differences in the ordering. We label these states by capital letters A, B, C, D, and E, where the ordering is defined by the CAM-B3LYP calculation (see Table 3). Four of these states can also be identified in the B3LYP calculations, whereas in the BLYP calculations only two of the states A–E are among the first five states calculated. Three of the states obtained with BLYP differ significantly in character from the ones obtained with the other exchange-correlation functionals.

Table 3

Excitation Energies (ΔE in eV) and Character of the Five Lowest States of PNA (A–E) Calculated with the aug-cc-pVDZ Basis Set at Different Levels of Theorya

HF			BLYP			B3LYP			CAM-B3LYP
ΔE			ΔE			ΔE			ΔE
3.209	A	CT	1.636	A	CT	1.982	A	CT	2.176	A	CT
4.598	B		2.990		CT,Ry	3.918	B		4.138	B
5.029	D	CT	3.613			4.003	C	CT	4.223	C	CT
5.246	C	CT	3.624		CT	4.269	D	CT	4.611	D	CT
5.322	E		3.843	C	CT	4.397			4.724	E

States of charge-transfer character are marked “CT”; Rydberg states are marked “Ry”.

Apart from one state in the BLYP calculation, none of the states have a significant Rydberg character. States A–E are always dominated by one orbital transition. The final orbital of all these transitions is the one marked (f) in Figure 4. The initial orbitals in these transitions are the orbitals (a), (b), (c), (d), and (e), respectively. The order of these orbitals varies slightly between the methods, but their shape is always the same. States A, C, and D show partial charge transfer from the nitro group into the rest of the aromatic system.

Figure 4

Orbitals dominating the excitation of PNA. (a)–(e): occupied orbitals dominating the states of characteristics A to E; (f): virtual orbital to which all excitations take place. All orbitals have been plotted at a contour value of 0.05. The orbitals were calculated using the CAM-B3LYP functional and the aug-cc-pVDZ basis set. Orbital plots have been rendered using the MOLDEN program.[54]

Comparing our state characteristics with others recently published in the literature, we find that especially the characteristics of the states from the CAM-B3LYP functional are in good agreement with results from coupled cluster calculations. Kosenkov and Slipchenko recently studied the PNA molecule in the gas phase using EOM-CCSD and the basis set 6-31+G(d).[55] For the first three excitations, they find the same order of states as in our CAM-B3LYP calculations, while for higher excitations the order is somewhat different. In particular, they also find low-lying Rydberg states, which we do not observe here among the first five states using CAM-B3LYP.

In Figure 5, we show the results for the multiphoton absorption cross section calculations for the four different methods on the PNA molecule.

Figure 5

Multiphoton absorption cross section of the PNA molecule calculated with the aug-cc-pVDZ basis set for states A–E.

State A—which is of charge-transfer character—has lower cross sections than the other states. Moreover, this state has the largest differences between HF and DFT with increasing differences with the number of photons absorbed. Also, HF has the highest cross section for OPA and the lowest for TPA, 3PA and 4PA. The π → π* state B has the largest cross sections in all cases. Comparing the results from the different methods, we note in particular that the deviation between B3LYP on the one hand and HF/CAM-B3LYP on the other increases with increasing number of photons, whereas the results from HF and CAM-B3LYP remain similar. CAM-B3LYP and HF also behave similarly for state E, which is the other π → π* state considered. For states C and D—both of charge-transfer character—we observe that the differences between the methods can be an order of magnitude but that they do not increase with the number of photons absorbed.

To examine the basis set dependence of the absorption cross sections, we have performed a series of calculations on the PNA molecule with different Dunning-style basis sets and the CAM-B3LYP functional. The results are presented in Figure 6.

Figure 6

Multiphoton absorption behavior of the PNA molecule using different basis sets and the CAM-B3LYP functional.

The largest differences are observed between augmented and nonaugmented basis sets, in particular, for 2PA, 3PA, and 4PA. The ordering of the energy of the states does not depend much on the choice of basis set, the only exception being states B and C, which change order when changing from the nonaugmented to the augmented basis sets. For the augmented basis sets, the curves are almost indistinguishable in Figure 6. We therefore conclude that the use of the aug-cc-pVDZ basis set for MPA calculations is an excellent compromise between computational efficiency and accuracy of the results.

para-Nitroaminostilbene

As for para-nitroaniline, we start by discussing the nature of the lowest excited states followed by a discussion of the multiphoton absorption cross sections of different order. The excitation energies of the excited states are collected in Table 4.

Table 4

Excitation Energies (ΔE in eV) and Character of the Five Lowest States of PNAS Calculated with the aug-cc-pVDZ Basis Set at Different Levels of Theorya

HF		BLYP		B3LYP		CAM-B3LYP
ΔE		ΔE		ΔE		ΔE
3.871	BF	2.257	BF	2.783	BF	3.400	BF
4.888	AF	3.265	AF	3.752	AF	3.935	AF
4.955		3.333	BG	3.833	BG	4.378	BG
4.984	BH	3.510		4.110	CF	4.473	CF
5.231		3.523		4.162		4.572	BH

For the two-letter code referring to the character of the states we refer to the text.

The correlation of states obtained with different exchange-correlation functionals is more difficult than for PNA, in part because there is a large mixture of different orbital transitions in each electronic excitation and in part because the different functionals give very different results. We will limit our discussion to states that can be found for more than one method; see Table 4. Figure 7 shows the relevant occupied orbitals and the dominating virtual orbitals are shown in Figure 8.

Figure 7

Occupied orbitals dominating the excitations of PNAS. All orbitals have been plotted at a contour value of 0.05. Orbitals (a) and (b) are taken from a Hartree–Fock calculation while orbital (c) is from a CAM-B3LYP calculation. All orbitals were calculated using the aug-cc-pVDZ basis set. Orbital plots have been rendered using the MOLDEN program.[54]

Figure 8

Virtual orbitals dominating the excitations of PNAS. Orbitals (f), (g), and (i) have been rendered at a contour value of 0.05. Orbital (h) has been rendered at a contour value of 0.01. Orbital (f) was taken from a Hartree—Fock calculation while (g), (h), and (i) are from a CAM-B3LYP calculation. All orbitals were calculated using the aug-cc-pVDZ basis set. Orbital plots have been rendered using the MOLDEN program.[54]

The S1- and S2-states are the same for all methods. They are dominated by the orbital transitions bf and af (following the indices of the corresponding orbitals in Figures 7 and 8, respectively), respectively, and the states will therefore be denoted BF and AF in the following. Both states are of charge-transfer character. The bf transition is π → π* with charge moving from the amino- to the nitro-substituted aromatic ring. The AF state is n → π* with charge moving from the nitro group to the nitro-substituted ring. The BG-state is only found in the DFT calculations where it is always the S3-state, and it follows from the figures that it is a π → π* state without charge-transfer character. The CF-state is only found for the B3LYP and CAM-B3LYP functionals. Its main contributions are the cf and ai transitions, and thus, it can be considered to have some small charge-transfer character. Finally, the BH state has Rydberg character and is only found in the HF and CAM-B3LYP calculations.

Diagrams for the one-, two-, three-, and four-photon absorption cross sections of the different states discussed are shown in Figure 9.

Figure 9

Multiphoton absorption cross sections for different states of the PNAS molecule calculated with the aug-cc-pVDZ basis set.

In general, we note that the behavior of the different absorption properties is quite similar. The BF and BG states show rather similar absorption cross sections for the different methods. The BH state—a Rydberg state only found for HF and CAM-B3LYP—shows larger cross sections for CAM-B3LYP by 1 to 3 orders of magnitude. The AF state is found to be the state with the lowest absorption strength for all four absorption processes considered here, with HF always yielding the smallest cross sections. For the CF state, the cross section from CAMB3LYP is always smaller by several orders of magnitude than the one from B3LYP. None of the states show differences between the methods that increase with the number of photons absorbed.

Conclusion and Outlook

We have presented an open-ended recursive approach for the calculation of single residues of response functions using SCF-based theory. The approach is so far restricted to perturbations that do not affect the basis set. This new functionality has been used to calculate transition matrix elements for one-, two-, three-, and four-photon absorption. We have also presented a way to calculate multiphoton absorption cross sections from the transition matrix elements. This work therefore provides a generalization of the theoretical treatment of multiphoton absorption and enables the calculation of multiphoton absorption to arbitrary order.

In the application part of this article, we have presented the first theoretical treatment of four-photon absorption and we have provided a comparison with lower-order absorption properties as well as an assessement of different SCF-based methods and basis sets. We have found that the calculated multiphoton absorption cross sections are not very sensitive to the size of the basis set as long as a reasonably large basis set with diffuse functions is used. The choice of method (HF or DFT) and the choice of exchange-correlation functional significantly affect the calculated cross section with differences extending over several orders of magnitude. These differences increase with the number of photons absorbed only for some of the investigated states. Charge-transfer states—which are frequently the brightest states for molecules with strong multiphoton absorption—were not found to behave in a different way from other states in this respect. We conclude from our calculations on PNA and PNAS that the combination of the range-separated hybrid CAM-B3LYP density functional with the aug-cc-pVDZ basis set is a good computational prescription for the treatment of multiphoton absorption. This finding is also based on the finding from several studies showing the reliability of CAM-B3LYP for the description of two-photon absorption.[21,56] These studies feature a comparison between a CAM-B3LYP and a coupled cluster treatment of two-photon absorption. As we found the behavior of the different j-photon absorption properties to be reasonably similar, we consider these results to be transferable to 3PA and 4PA. Considering the results of this work, we would in general recommend to use range-separated hybrid functionals for the treatment of multiphoton absorption properties.

This work lays the foundation for computational chemistry to follow the exciting developments happening experimentally in the field of multiphoton spectroscopy.[11−15] Through computational modeling and a detailed understanding of the underlying quantum-mechanical effects determining multiphoton absorption cross sections, important insight into how to design molecules with large multiphoton absorption cross sections can be obtained. Work remains in being able to treat the effects of a surrounding chemical environment such as a solvent or a protein, as well as effects due to vibrational contributions for which the formalism needs to be extended also to perturbation-dependent basis sets. Work along these lines is in progress in our research group.