# Cubic response#

## Second-order electric-dipole hyperpolarizability#

The second-order electric-dipole hyperpolarizability relates to the cubic response function accordin to

$\gamma_{\alpha\beta\gamma\delta}(-\omega_\sigma;\omega_1,\omega_2,\omega_3) = -\langle\langle \hat{\mu}_\alpha ;\hat{\mu}_\beta ,\hat{\mu}_\gamma ,\hat{\mu}_\delta \rangle\rangle_{\omega_1,\omega_2,\omega_3}$

where $$\omega_\sigma = \omega_1 + \omega_2 + \omega_3$$. In the limit of static fields, $$\gamma(0;0,0,0)$$ can also be determined from the fourth-order derivative of the energy with respect to electric field strengths.

### Pertinent optical processes#

For dynamic fields, nonlinear optical processes of interest relating to this molecular property include:

Optical process

Response function

Electric-Field Induced Second-Harmonic Generation

$$\gamma^\mathrm{ESHG}_{\mu\nu\eta\xi}(-2\omega;\omega,\omega,0) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi} \rangle\rangle_{\omega,\omega,0}$$

Third Harmonic Generation

$$\gamma^\mathrm{THG}_{\mu\nu\eta\xi}(-3\omega;\omega,\omega,\omega) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi} \rangle\rangle_{\omega,\omega,\omega}$$

Electro-optical or dc-Kerr effect

$$\gamma^{\rm dc-Kerr}_{\mu\nu\eta\xi}(-\omega;\omega,0,0) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi}\rangle\rangle_{\omega,0,0}$$

Optical or ac-Kerr effect

$$\gamma^\mathrm{ac-Kerr}_{\mu\nu\eta\xi}(-\omega_1;\omega_1,\omega_2,-\omega_2) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi}\rangle\rangle_{\omega_1,\omega_2,-\omega_2}$$

Intensity Dependent Refractive Index

$$\gamma^\mathrm{IDRI}_{\mu\nu\eta\xi}(-\omega;\omega,\omega,-\omega) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi}\rangle\rangle_{\omega,\omega,-\omega}$$

dc Optical Rectification

$$\gamma^\mathrm{dc-OR}_{\mu\nu\eta\xi}(0;\omega,-\omega,0) = - \langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta},\hat{\mu}_{\xi} \rangle\rangle_{\omega,-\omega,0}$$

We note that IDRI is also known under the name of degenerate four-wave mixing (DFWM).

### Isotropic quantities#

The experimentally measured quantities in isotropic fluids are usually the scalar components of the tensor $$\gamma$$ given by the isotropic average:

$\begin{eqnarray*} \gamma_{\parallel} & = & \frac{1}{15} (\gamma_{\xi\eta\eta\xi} +\gamma_{\xi\eta\xi\eta}+\gamma_{\xi\xi\eta\eta}), \\ \gamma_{\perp} & = & \frac{1}{15} (2\gamma_{\xi\eta\eta\xi} -\gamma_{\xi\xi\eta\eta}), \end{eqnarray*}$

where the same sequence of optical frequencies, not given explicitly above, is used for the laboratory axes and the molecular axes.

As for the first hyperpolarizability, the number of independent non-zero tensor elements depends on the optical process and on the symmetry of the molecule. For example, $$\gamma^{\rm THG}$$ is symmetric with respect to the second, third, and fourth indices, and this can be used to simplify the expression given for $$\gamma_{\perp}$$. For instance, the average value that should be compared with the experimental THG parallel component becomes

$\begin{eqnarray*} \gamma^{\rm THG}_{\|} = \frac{1}{5} \gamma^{\rm THG}_{\xi\xi\eta\eta}. \end{eqnarray*}$

Experimentally, if all fields have parallel polarization, one can measure the parallel components of the first and second hyperpolarizablities, which take into account the classical orientational averaging. In the case of ESHG, with the optical field polarized perpendicular to the static field, one measures the perpendicular components, and in the case of a dc-Kerr experiment, the differences between the parallel and perpendicular components.

In the ESHG experiment a laser beam passes through the sample in a static electric field and a weak, collinear, frequency-doubled beam is detected.
Absolute values for the hyperpolarizabilities cannot be extracted; the signal from the sample is compared to that of a known buffer gas (ultimately helium, for which there are accurate theoretical values) or a solid. In analogy with the derivation for $$\tilde{\alpha}$$, the classical thermal averaging yields

$\tilde{\gamma}^{\rm ESHG} = \left( \tilde{\beta}_{ZZZ}^{F}/{F} \right)_{{F}{\rightarrow}0} = \langle \gamma_{ZZZZ} \rangle_\Omega + \mu_z \langle \beta_{ZZZ} \rangle_\Omega /3kT.$

The molecular hyperpolarizabilities $$\gamma^{\rm ESHG}$$ and $$\beta^{\rm SHG}$$ can be obtained in this experiment. For non-centrosymmetric molecules a series of measurements over a range of temperatures has to be performed, whereas for centrosymmetric molecules $$\beta = 0$$ and thus a single measurement at one temperature is sufficient.

The majority of existing experimental data on hyperpolarizabilities are derived from ESHG and dc-Kerr measurements. The dc-Kerr effect differs from the other nonlinear optical processes as it allows for absolute measurements without the need for a reference measurement. The measured molar Kerr constant is

$\begin{eqnarray*} A^{\rm dc-Kerr} &=& \frac{N_A}{81\epsilon_0} \left[ \frac{3}{2} (\langle \gamma_{ZZZZ} \rangle_\Omega-\langle \gamma_{ZXXZ} \rangle_\Omega) +\mu_z (\langle \beta_{ZZZ} \rangle_\Omega - \langle \beta_{XXZ} \rangle_\Omega ) /kT \right. \\ &+& \left. \frac{3}{10kT} \left( \alpha_{\alpha\beta}\alpha_{\alpha\beta}^0 - \langle \alpha_{ZZ} \rangle_\Omega \langle \alpha_{ZZ}^0 \rangle_\Omega + \frac{1}{kT} \mu_z^2 ( \alpha_{zz} - \langle \alpha_{ZZ} \rangle_\Omega ) \right) \right], \end{eqnarray*}$

where $$N_A$$ is Avogadro’s number, and it involves both dynamic hyperpolarizabilities ($$\gamma^{\rm dc-Kerr}$$ and $$\beta^{\rm EOPE}$$) as well as static (superscript 0) and dynamic polarizabilities.

The computational requirements for the second hyperpolarizability are more or less the same as for the first hyperpolarizability. Although the number of studies that analyze the importance of electron correlation effects (not including here DFT) is rather limited, the available results confirm in general the findings made for the first hyperpolarizabilities.

The atomic units of $$\gamma$$ are e$$^4$$ a$$_0^4$$ E$$_h^{-3}$$. 1 au of $$\gamma$$ corresponds to 6.23538 $$\times 10^{-65}$$ C$$^4$$ m$$^4$$ J$$^{-3}$$ (SI) or to 5.0367 $$\times 10^{-40}$$ statvolt$$^{-2}$$ cm$$^5$$ (esu).

### Dispersion#

Similar to Cauchy expansions, expansions may be applied to describe the dispersions of hyperpolarizabilities

$\gamma_{\parallel}(-\omega_\sigma;\omega_1,\omega_2,\omega_3) = \gamma_\parallel(0;0,0,0) \Big[ 1 + {\cal A}^{\gamma}(\omega_\sigma^2+\omega_1^2+\omega_2^2+\omega_3^2) +{\cal{O}}(\omega_i^4) \Big]$

Since the expansion coefficient, $${\cal A}^{\gamma}$$, is frequency independent in itself, dispersions are transferable in between different nonlinear optical processes.

## In VeloxChem#

As an illustration, let us calculate $$\gamma_{zzzz}(-\omega; \omega, -\omega, \omega)$$ using damped response theory.

import veloxchem as vlx


First, we determine the reference state of the system.

mol_str = """3

O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(molecule, basis)


Second, we calculate the needed cubic response functions with use of the CubicResponseDriver class.

crf_drv = vlx.CubicResponseDriver()

crf_drv.a_operator = "electric dipole"
crf_drv.b_operator = "electric dipole"
crf_drv.c_operator = "electric dipole"
crf_drv.d_operator = "electric dipole"

# other available operators
#lrf.b_operator = "electric dipole"
#lrf.b_operator = "linear momentum"
#lrf.b_operator = "angular momentum"
#lrf.b_operator = "magnetic dipole"

crf_drv.a_component = "z"
crf_drv.b_component = "z"
crf_drv.c_component = "z"
crf_drv.d_component = "z"

crf_drv.b_frequencies = [0.0, 0.0656]
crf_drv.c_frequencies = [0.0, -0.0656]
crf_drv.d_frequencies = [0.0, 0.0656]

crf_drv.damping = 0.004556  # 1000 cm-1

crf_results = crf_drv.compute(molecule, basis, scf_results)

Hide code cell output

Cubic Response Driver Setup
=============================

ERI Screening Threshold         : 1.0e-12
Convergance Threshold           : 1.0e-04
Max. Number of Iterations       : 150
Damping Parameter               : 4.556000e-03
Exchange-Correlation Functional : B3LYP
Molecular Grid Level            : 4



Complex Response Solver Setup
===============================

Number of Frequencies           : 1
Max. Number of Iterations       : 150
Convergence Threshold           : 1.0e-04
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Threshold         : 1.0e-12
Exchange-Correlation Functional : B3LYP
Molecular Grid Level            : 4


* Info * Molecular grid with 40720 points generated in 0.02 sec.

* Info * Processing Fock builds... (batch size: 5)
* Info *   batch 1/1


* Info * 2 gerade trial vectors in reduced space
* Info * 3 ungerade trial vectors in reduced space

* Info * 75.47 kB of memory used for subspace procedure on the master node
* Info * 2.58 GB of memory available for the solver on the master node

*** Iteration:   1 * Residuals (Max,Min): 8.43e-01 and 8.28e-01


* Info * Processing Fock builds... (batch size: 6)
* Info *   batch 1/1


* Info * 5 gerade trial vectors in reduced space
* Info * 6 ungerade trial vectors in reduced space

* Info * 91.79 kB of memory used for subspace procedure on the master node
* Info * 2.58 GB of memory available for the solver on the master node

*** Iteration:   2 * Residuals (Max,Min): 9.39e-02 and 9.19e-02


* Info * Processing Fock builds... (batch size: 4)
* Info *   batch 1/1


* Info * 7 gerade trial vectors in reduced space
* Info * 8 ungerade trial vectors in reduced space

* Info * 102.67 kB of memory used for subspace procedure on the master node
* Info * 2.58 GB of memory available for the solver on the master node

*** Iteration:   3 * Residuals (Max,Min): 9.72e-03 and 9.46e-03


* Info * Processing Fock builds... (batch size: 6)
* Info *   batch 1/1


* Info * 10 gerade trial vectors in reduced space
* Info * 11 ungerade trial vectors in reduced space

* Info * 118.99 kB of memory used for subspace procedure on the master node
* Info * 2.58 GB of memory available for the solver on the master node

*** Iteration:   4 * Residuals (Max,Min): 3.69e-04 and 3.61e-04


* Info * Processing Fock builds... (batch size: 4)
* Info *   batch 1/1


* Info * 12 gerade trial vectors in reduced space
* Info * 13 ungerade trial vectors in reduced space

* Info * 129.87 kB of memory used for subspace procedure on the master node
* Info * 2.58 GB of memory available for the solver on the master node

*** Iteration:   5 * Residuals (Max,Min): 1.18e-05 and 1.17e-05


               *** Complex response converged in 5 iterations. Time: 1.28 sec


* Info * Molecular grid with 40720 points generated in 0.02 sec.

Fock Matrix Computation
=========================

* Info * Processing Fock builds... (batch size: 4)
* Info *   batch 1/4

* Info *   batch 2/4

* Info *   batch 3/4

* Info *   batch 4/4


* Info * Time spent in Fock matrices: 2.91 sec



Complex Response Solver Setup
===============================

Number of Frequencies           : 1
Max. Number of Iterations       : 150
Convergence Threshold           : 1.0e-04
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Threshold         : 1.0e-12
Exchange-Correlation Functional : B3LYP
Molecular Grid Level            : 4


* Info * Molecular grid with 40720 points generated in 0.02 sec.

* Info * Processing Fock builds... (batch size: 8)
* Info *   batch 1/1


* Info * 4 gerade trial vectors in reduced space
* Info * 4 ungerade trial vectors in reduced space

* Info * 67.66 kB of memory used for subspace procedure on the master node
* Info * 2.59 GB of memory available for the solver on the master node

*** Iteration:   1 * Residuals (Max,Min): 9.76e-02 and 8.18e-02


* Info * Processing Fock builds... (batch size: 8)
* Info *   batch 1/1


* Info * 8 gerade trial vectors in reduced space
* Info * 8 ungerade trial vectors in reduced space

* Info * 89.42 kB of memory used for subspace procedure on the master node
* Info * 2.59 GB of memory available for the solver on the master node

*** Iteration:   2 * Residuals (Max,Min): 6.61e-03 and 5.82e-03


* Info * Processing Fock builds... (batch size: 8)
* Info *   batch 1/1


* Info * 12 gerade trial vectors in reduced space
* Info * 12 ungerade trial vectors in reduced space

* Info * 111.18 kB of memory used for subspace procedure on the master node
* Info * 2.59 GB of memory available for the solver on the master node

*** Iteration:   3 * Residuals (Max,Min): 4.85e-04 and 4.28e-04


* Info * Processing Fock builds... (batch size: 8)
* Info *   batch 1/1


* Info * 16 gerade trial vectors in reduced space
* Info * 16 ungerade trial vectors in reduced space

* Info * 132.94 kB of memory used for subspace procedure on the master node
* Info * 2.59 GB of memory available for the solver on the master node

*** Iteration:   4 * Residuals (Max,Min): 1.67e-05 and 1.55e-05


               *** Complex response converged in 4 iterations. Time: 1.24 sec


                                                 Fock Matrix Computation
=========================

* Info * Processing Fock builds... (batch size: 4)
* Info *   batch 1/1


* Info * Time spent in Fock matrices: 0.62 sec



Cubic response function: << z;z,z,z >>  (0.0,0.0,0.0)
=======================================================

Real             Imaginary
----------------------------------------------------
CRF              -998.63439127           0.00000000j

Cubic response function: << z;z,z,z >>  (0.0656,-0.0656,0.0656)
=================================================================

Real             Imaginary
----------------------------------------------------
CRF             -1152.42872036         -12.04492771j

*** Time spent in cubic response calculation: 6.13 sec ***



Third, we print the requested molecular properties.

print(f"{'w1':>8s}{'w2':>8s}{'w3':>8s}{'gamma_zzzz':>24s}")
print("-"*48)
for key in crf_results.keys():
print(f"{key:8.4f}{key:8.4f}{key:8.4f}{-crf_results[key]:24.4f}")

      w1      w2      w3              gamma_zzzz
------------------------------------------------
0.0000  0.0000  0.0000        998.6344-0.0000j
0.0656 -0.0656  0.0656      1152.4287+12.0449j


Note

Taking into account the three perturbation operators, all of the form $$\hat{V}^\omega = - \hat{\mu}_\alpha F_\alpha^\omega$$, we need to change the sign of the response function values to get results for the second-order hyperpolarizability.