Simple Finite Field Nuclear Relaxation Method for Calculating Vibrational Contribution to Degenerate Four-wave Mixing

A simple extended finite field nuclear relaxation procedure for calculating vibrational contributions to degenerate four-wave mixing ͑also known as the intensity-dependent refractive index͒ is presented. As a by-product one also obtains the static vibrationally averaged linear polarizability, as well as the first and second hyperpolarizability. The methodology is validated by illustrative calculations on the water molecule. Further possible extensions are suggested.


I. INTRODUCTION
It is now well established that vibrational contributions to nonlinear optical ͑NLO͒ properties can be very significant [1][2][3][4] for most important NLO processes.About 15 years ago, Bishop and Kirtman presented a general perturbation theory ͑BK-PT͒ treatment, [5][6][7] which allows these contributions ͑as well as the vibrational linear polarizability͒ to be calculated for polyatomic molecules taking into account both mechanical and electrical anharmonicity.][12][13][14][15] For chemical systems of interest such as NLO materials, the requirement on the lowest electronic excitation energy is typically well satisfied.
In addition to BK-PT, alternative approaches continue to emerge.A variational method for linear polarizabilities, based on analytical linear response theory, 8,9 was presented very recently by Christiansen et al.For static and electrooptical hyperpolarizabilities ͑specified below͒ an earlier, finite field-nuclear relaxation ͑FF-NR͒ alternative formulation has proved valuable.The latter is related to BK-PT but avoids explicit calculation of high-order potential energy and electrical property derivatives, [16][17][18][19] which are otherwise required.7][18][19] Subsequent calculation of the electronic energy or dipole moment, linear polarizability, and first hyperpolarizability at the relaxed geometry yields the lowest-order ͑to be defined later͒ vibrational perturbation terms, [16][17][18] while all remaining vibrational terms are obtained from the zero-point vibrational average ͑ZPVA͒ of these properties. 19The ZPVA, in turn, can be obtained by applying variational procedures to solve the vibrational Schrödinger equation. 20,21perimental measurements are often carried out at optical frequencies much larger than the frequencies of the relevant vibrational excitations.In that event, one may invoke the infinite optical frequency approximation [22][23][24] ͑IOFA͒ in calculating vibrational NLO properties.Although not strictly required, 25 usually the FF-NR treatment is applied in conjunction with the IOFA.In that event, this method requires only the calculation of static electronic and ZPVA ͑hyper͒polarizabilities at the field induced geometries, which implies a major computational simplification.
The vibrational contributions to electro-optical NLO properties that can be computed by the FF-NR procedure include ␤͑− ; ,0͒ ͓dc-Pockels ͑dc-P͒ effect͔; ␥͑− ; ,0,0͒ ͓dc-Kerr ͑dc-K͒ effect͔; and ␥͑−2 ; , ,0͒ ͓electric-field-induced second harmonic generation ͑ESHG͔͒; as well as the static ␣, ␤, and ␥.An important purely dynamic NLO property that has, thus far, not proved accessible by the FF-NR approach is ␥͑− ; ,− , ͒ ͓degenerate four-wave mixing ͑DFWM͒, also known as the intensitydependent refractive index ͑IDRI͔͒.An approximate expression for the vibrational contribution to DFWM, based on the FF-NR approach, was presented some time ago. 16However, that expression is correct only through the first-order of BK-PT.Subsequent calculations have shown that second-order contributions from the ␥͑0;0,0,0͒ term in this approximate formula are typically quite large. 17,26,27For that reason, it has never been used and there are no published calculations of DFWM beyond the simple double harmonic approximation.Since DFWM is susceptible to large vibrational contributions and has a number of useful applications, it is of interest to remedy that situation.Some of the more noteworthy applications derive from optical bistability and self-focusing. 28,29he former, in particular, serves as the basis for devices used in optical memory elements or switches for optical communications and optical computing.
The purpose of this paper is to show how the FF-NR procedure may be extended to bring vibrational DFWM under the same rubric as the other NLO properties listed above.In addition, we obtain a simple way to calculate the static a͒ Electronic mail: kirtman@chem.ucsb.edu.
ZPVA of ␣, ␤, and ␥.As already noted, ␣ ZPVA and ␤ ZPVA are used to determine the three electro-optical NLO processes mentioned above, whereas ␥ ZPVA is required for our new procedure in order to isolate the pure vibrational DFWM.In Sec.II, we set out the methodology and, then, Sec.III provides illustrative calculations that serve to validate our treatment.Finally, in Sec.IV, we summarize and speculate about possible future developments.

II. THEORY
According to BK-PT, the vibrational frequencydependent second hyperpolarizability may be written in the general form where 1 , 2 , and 3 are the optical input frequencies and the output frequency = 1 + 2 + 3 .DFWM corresponds to the case ␥͑− ; ,− , ͒.The square bracket quantities on the right-hand side of Eq. ͑1͒ are functions of dipolar electrical property derivatives, harmonic vibrational frequencies, vibrational anharmonicity constants, and optical frequencies.
The symbols inside the square brackets indicate which electronic dipolar properties are involved.Although the electrical properties involved are the same for each NLO process, i.e., for each choice of 1 , 2 , and 3 , the specific formulas will vary.In the IOFA limit, however, the square bracket terms differ by, at most, a multiplicative constant.This constant may be zero so that some of the contributions vanish in that limit.Indeed, for DFWM, only the ͓␣ 2 ͔ term survives.The three remaining square brackets contain electric dipole derivatives.Thus, if the dipole derivative terms are subtracted out of the Hamiltonian when calculating some other NLO process by the FF-NR approach, then what remains will be the desired DFWM square bracket term.Since ͓␣ 2 ͔ does not contribute to ESHG in the IOFA limit, the one electro-optical process that could be employed for our purposes is the dc-K effect.This property is obtained from the electronic linear polarizability and ␣ ZPVA .However, it is simpler to use the static ␥, which is determined from the electronic and zero-point vibrational energy, ZPVE ͑or the electronic and ZPVA ͒.In that case, the ͓␣ 2 ͔ term is multiplied by 2 / 3 to get the correct DFWM result.So this is what we do.
There are two major steps involved in the FF-NR approach.The first is to calculate the relaxed geometry in a finite static field with all dipole derivative terms removed from the Hamiltonian.This goal can be easily achieved by subtracting F a a 01 a ͑Q , 0͒ from the electronic energy and F a a 11 a ͑Q , 0͒ from the ith component of the gradient vector.Here, and F a is the static electric field in the a direction.The quantity a 01 a ͑Q , 0͒, for example, is just the a component of the electronic dipole moment at the geometry corresponding to the normal coordinate displacement Q and at zero field.Subsequently, the nuclear relaxation contribution to DFWM, which is the zeroth-order or double harmonic approximation to ͓␣ 2 ͔, is obtained by calculating the electronic energy ͑or dipole moment͒ at the field-induced ͑relaxed͒ equilibrium geometry. 16,17This must be repeated at several different fields.The energy is then fitted to a power series in the field and the fourth-order term in that power series is the desired nuclear relaxation ␥ with electric dipole terms removed.
There is an alternative that is far more efficient for obtaining the relaxed geometries.It uses just the parameters a nm i,j,. ..,a,b. . .͑0 , 0͒, which are determined once and for all at the field-free equilibrium geometry.From our previous analysis, 18,19,26 it is easily seen that only the first and second dipole derivatives with respect to field-free normal modes Q i enter into the field-dependent equilibrium geometry calculation.Hence, we subtract from the energy and from ith element of the gradient vector.If one is interested solely in calculating the nuclear relaxation contribution to DFWM, then just the first terms of expressions ͑3͒ and ͑4͒ are required since the second terms contain anharmonic parameters.In order to calculate the higher-order BK-PT contributions to ͓␣ 2 ͔ a second major step is necessary.It involves evaluating the field-dependent ZPVE ͑or ZPVA ͒ at the relaxed geometry.1][32] They require fitting the potential energy surface ͑PES͒ at specified geometries Q.In order to obtain the exact DFWM ͑within the IOFA͒, it is necessary to subtract the dipolar contribution F a a 01 a ͑Q , 0͒ at each of these points.
An approximate treatment, which requires much less computational effort, utilizes only the leading fielddependent ZPVE term given ͑in a.u.͒ by where Q F is the field-dependent relaxed geometry.In the FF-NR procedure, this approximate ZPVE will give rise to the next highest nonvanishing BK-PT contributions beyond those included in the nuclear relaxation formula.For example, the second-order ͓␣ 2 ͔ term is determined in this manner ͑since the first-order term vanishes͒.Likewise, one obtains the second-order ͓␤͔, the third-order ͓ 2 ␣͔, and the fourth-order ͓ 4 ͔. 17Note that, again, we want to remove all dipolar contributions.Since a 20 ii ͑Q F , F͒ in Eq. ͑5͒ is obtained by diagonalizing the field-dependent Hessian, this implies that F a a 21 ij,a ͑Q F , 0͒ should be subtracted from the ͑i , j͒ element of the Hessian before calculating the approximate fielddependent ZPVE of Eq. ͑5͒.
The same general approach can be utilized to obtain the static ␣ ZPVA , ␤ ZPVA , and ␥ ZPVA ͑ ZPVA is already obtained in the standard FF-NR method͒.These quantities appear in the expression for ⌬E ZPVA as calculated by our new procedure, This expression immediately gives the static ␣ ZPVA and ␤ ZPVA , which are important in their own right and are also needed to extract higher-order vibrational contributions to the dc-P, dc-K, and ESHG processes, as described in Sec.I. Clearly, the static ␥ ZPVA is needed in order to extract the higher-order vibrational contributions to DFWM ͑denoted ␥ c-ZPVA here to correspond with our previous papers͒ from Eq. ͑6͒.Again, the first step in determining this quantity is to evaluate geometry relaxation due to the field.Besides removing all dipole derivatives from the Hamiltonian, now the linear polarizability derivatives must be deleted as well.Thus, the terms F a a 01 a ͑Q , 0͒ + F a 2 a 02 aa ͑Q , 0͒ and F a a 11 i,a ͑Q , 0͒ + F a 2 a 12 i,aa ͑Q , 0͒ should be removed from the electronic energy and ith component of the gradient vector, respectively, during the field-dependent geometry optimization.Again, there is a more efficient alternative formulation based on the equilibrium a nm i,j,. ..,a,b. . .͑0 , 0͒ parameters.Within this formulation is subtracted from the electronic energy and from ith element of the gradient vector.Finally, the second step again involves evaluating the ZPVE at the relaxed geometry with the dipolar and linear polarizability terms removed from the PES.If only the harmonic ZPVE is computed ͓i.e., Eq. ͑5͔͒, we must remove F a a 21 ij,a ͑Q F , 0͒ + F a 2 a 22 ij,aa ͑Q F , 0͒ from the ͑i , j͒ element of the fielddependent Hessian.
In summary, to evaluate the exact ␥ c-ZPVA ͑− ; , − , ͒ and exact static ␣ ZPVA , ␤ ZPVA , and ␥ ZPVA as efficiently as possible using our new FF-NR approach, the following quantities are required for the field-dependent geometry optimizations: a 11 i,a ͑0 , 0͒, a 21 ij,a ͑0 , 0͒, and a 12 i,aa ͑0 , 0͒.In addition, the ZPVE is obtained using a PES modified at each point Q by subtracting terms that involve the quantities a 01 a ͑Q , 0͒ and a 02 aa ͑Q , 0͒.If the harmonic approximation is employed for ZPVE, then only the reduced set of quantities: a 20 ii ͑Q F , F͒, a 21 ij,a ͑Q F , 0͒, and a 22 ij,aa ͑Q F , 0͒ is required.The latter approximation leads to vibrational NLO contributions through second order in BK-PT.In comparison, to evaluate the same properties ͑through second order͒ using the BK-PT method, it is necessary to calculate a 20 ii ͑0 , 0͒, a 30 ijk ͑0 , 0͒, ii,aaa ͑0 , 0͒, a 14 i,aaaa ͑0 , 0͒, and a 24 ii,aaaa ͑0 , 0͒.The cost of calcu-lating high-order normal coordinate derivatives, such as a 32 iij,aa ͑0 , 0͒ and a 40 iijk ͑0 , 0͒, increases dramatically with the number of atoms.For that reason, the BK-PT method is practical only for small molecules.Finally, although we have presented FF-NR expressions only for diagonal components, exactly analogous procedures can be applied for the nondiagonal components as well.

III. ILLUSTRATIVE EXAMPLE
The primary goal of this section is to validate our new methodology rather than to calculate benchmark values for comparison with experiment or previous literature results.It is noteworthy, nevertheless, that we report vibrational contributions to DFWM calculated for the first time beyond the double harmonic approximation ͑denoted here by nr as in previous usage͒.Our calculations were done on the water molecule at the HF and MP2 levels with the d-aug-cc-pVTZ basis set.For the sake of simplicity, we used the approximation of Eq. ͑5͒, and only the diagonal component of the property ͑in the direction parallel to the dipole moment͒ was considered.Table I compares our new FF-NR method with BK-PT for the nr and P c-ZPVA ͑second-order͒ contribution to ␥ zzzz ͑− ; ,− , ͒ →ϱ , as well as the ZPVA for the static linear polarizability ͑␣ zz ͑0;0͒͒, first hyperpolarizability ͑␤ zzz ͑0;0,0͒͒, and second hyperpolarizability ͑␥ zzzz ͑0;0,0,0͒͒.It is clear from the table that there is very good agreement between the FF-NR and BK-PT results.The small differences for P c-ZPVA are of the same magnitude as the numerical error involved in calculating this property.
The data in Table I support the well-established conclusion that electron correlation can have a very strong effect on the vibrational contribution.Usually the importance of vibrational contributions decreases when the number of static fields involved in the property decreases.The property ␥͑− ; ,− , ͒ →ϱ is an exception because there is a cancellation of + and − in energy denominators of the sumover-states formula, which simulates a static field.Indeed, at the MP2 level the recently published 21 IOFA ratio ͑␥ nr + ␥ c-ZPVA ͒ / ␥ zzzz e ͑0;0,0,0͒ is 0.04 for the dc-K process ͑two static fields͒ and 0.01 for ESHG ͑one static field͒, whereas the corresponding IDRI ratio found here is 0.09 ͑although the Sadlej basis 33 was used in Ref.21, rather than d-aug-cc-pVTZ͒.The corresponding ratio for the static property is larger yet ͑0.21͒, but that is not always the case.Finally, as an indicator for the convergence of the perturba- tion treatment of anharmonicity we usually use the ratio P c-ZPVA / P nr . 20In agreement with our previous published results for the static and electrooptical ͑hyper͒polarizabilities of water, 20,21 the ratio obtained here for DFWM is less than 0.18, which indicates that the perturbation series is converging rapidly.

IV. SUMMARY AND FUTURE OUTLOOK
We have presented a simple extension of our FF-NR procedure for calculating vibrational nonlinear optical contributions to the DFWM ͑or IDRI͒ process.The basic idea is to subtract out unwanted electric dipole terms from our previous formulation for ͑static͒ properties obtained from the zero point vibrational energy E ZPVA .This new methodology also leads to a simple procedure for evaluating the static ␣ ZPVA , ␤ ZPVA , and ␥ ZPVA .The first two of these can be utilized in our previous treatment of the dc-P, dc-K, and ESHG properties; the last of these is used in extracting the vibrational DFWM from the E ZPVA calculation.Our new procedure is validated by application to the water molecule using the lowest-order approximation for E ZPVA .
Future plans call for generalizing the extended FF-NR procedure to calculate dc-P, dc-K, and ESHG vibrational contributions from E ZPVA .Currently, we obtain these properties from the static ␣ ZPVA and ␤ ZPVA , which are computationally more expensive to determine than E ZPVA .In addition, we want to incorporate accurate methods for obtaining the static E ZPVA , such as vibrational CI, [34][35][36] within our procedure.The possibility that an extended FF-NR approach can be adapted for vibrational two photon absorption 37 will also be considered.

TABLE I .
Comparison between BK-PT and our new FF-NR method for the diagonal component ͑along the direction parallel to ͒ of various properties calculated for water at the HF/ and MP2/d-aug-cc-pVTZ levels using the approximation of Eq. ͑5͒.All quantities are in a.u.