# Frequency-Dependent Quadratic Response Properties and Two-photon Absorption
from Relativistic Equation-of-Motion Coupled Cluster Theory

Xiang Yuan  [  Loïc Halbert  [  Lucas Visscher  l.visscher@ vu.nl [  André
Severo Pereira Gomes<EMAIL_ADDRESS>[

###### Abstract

We present the implementation of quadratic response theory based upon the
relativistic equation-of-motion coupled cluster method. We showcase our
implementation, whose generality allows us to consider both time-dependent and
time-independent electric and magnetic perturbations, by considering the
static and frequency-dependent hyperpolarizability of hydrogen halides (HX, X
= F-At), providing a comprehensive insight into their electronic response
characteristics. Additionally, we evaluated the Verdet constant for noble
gases Xe and Rn, and discussed the relative importance of relativistic and
electron correlation effects for these magneto-optical properties. Finally, we
calculate the two-photon absorption cross-sections of transition
($ns^{1}S_{0}\to(n+1)s^{1}S_{0}$) of Ga+, and I:n+, which are suggested as
candidates for new ion clocks. As our implementation allows for the use of
non-relativistic Hamiltonians as well, we have compared our EOM-QRCC results
to the QR-CC implementation in the DALTON code, and show that the differences
between CC and EOMCC response are in general smaller than 5% for the
properties considered. Collectively, the results underscore the versatility of
our implementation and its potential as a benchmark tool for other
approximated models such as density functional theory for higher-order
properties.

###### keywords:

American Chemical Society, LaTeX

Universite de Lille, CNRS, UMR 8523 − PhLAM] Physique des Lasers Atomes et
Molecules, Universite de Lille, F-59000 Lille, France [Vrije Universiteit
Amsterdam] Department of Chemistry and Pharmaceutical Sciences, Faculty of
Science, Vrije Universiteit Amsterdam, 1081 HV Amsterdam, The Netherlands
Universite de Lille, CNRS, UMR 8523 − PhLAM] Physique des Lasers Atomes et
Molecules, Universite de Lille, F-59000 Lille, France Vrije Universiteit
Amsterdam] Department of Chemistry and Pharmaceutical Sciences, Faculty of
Science, Vrije Universiteit Amsterdam, 1081 HV Amsterdam, The Netherlands
Universite de Lille CNRS, UMR 8523 − PhLAM] Physique des Lasers Atomes et
Molecules, Universite de Lille, F-59000 Lille, France IR,NMR,UV

##  1 Introduction

Nonlinear optical properties (NLO) of matter provide a wealth of information
on intra- and inter-molecular interactions and are therefore widely studied in
science and engineering1, 2, 3. NLO properties are also central to materials
and device design, with numerous important applications such as optical
devices for data transfer and storage. Among the materials being considered,
there is a growing interest in NLO properties of molecules containing heavier
elements, particularly in Lanthanide4, 5, 6, 7, 8, 9 and Actinide10, 11, 12
complexes, as they can offer a superior performance compared to molecules that
contain only light elements.

In order to compute and analyze molecular properties in the linear and non-
linear regime, one typically resorts to response theory13, 14, 15, 16. Within
this theory, the first-order nonlinear response is characterized by the
quadratic response function. Quadratic response functions have been
implemented for Hartree-Fock (HF) wave-functions17, 18 as well as at the
electron correlated level employing second-order Møller–Plesset perturbation
(MP2)19, multiconfigurational self-consistent field (MCSCF)20, 21, coupled
cluster (CC)22, 23, 24, and density functional theory (DFT)25, 26 reference
states. The common starting point of these developments has been the non-
relativistic molecular Hamiltonian.

For property calculations, the spin-orbit coupling operator can be added as
one of the perturbing properties. That will provide accurate results for
lighter elements, at the expense of needing to go one order higher in the
responses that are considered. As we move down the periodic table, we reach a
point where relativistic effects are too strong to be reliably treated as
perturbations. In this domain, it is therefore necessary to refine these
methods, ensuring relativistic effects are intrinsically accounted for by
employing a variationally stable relativistic Hamiltonian.

In the domain of relativistic quantum chemistry, to date, quadratic response
function derivations and implementations are primarily based on mean-field
models, such as HF27 and DFT28. To improve precision and establish benchmarks
for other models, in this manuscript we discuss the development of quadratic
response theory based on a relativistic equation-of-motion (EOM) coupled
cluster formulation (EOM-QRCC).

We showcase the generality and versatility of our implementation by examining
two molecular properties. First, we study the frequency-(in)dependent electric
first hyperpolarizability ($\beta$) as it can describe the nonlinear response
of a molecule to an applied electric field, which is significant for second-
harmonic generation29 associated with the design of optoelectronic devices and
can provide valuable insights into the intermolecular interaction30. For
instance, as discussed by Datta and Pati 31 $\beta$ is related to the weak
intermolecular forces such as dipolar interactions and hydrogen-bonding, thus
it is possible to control $\beta$ by modifying the interactions and accurate
calculations would be instrumental to provide insight into designing NLO
materials like $\pi$-conjugated molecular assemblies.

We consider magnetic circular birefringence, also known as the Faraday effect,
as the second property. One example of the interest in studying the Faraday
effect can be found in the observation by Savukov et al. 32 of the inverse
Faraday effect in the nuclear magnetic resonance (NMR) sample of liquid water
and liquid 129Xe , which has led to the suggestion that the nuclear spin-
induced optical rotation (NSOR) can provide a viable and potentially more
informative analog to the NMR chemical shift of traditional NMR detection.
There have been only a handful of theoretical investigations of this property,
however. For 129Xe, Ikäläinen et al. 33 performed non-relativistic (NR) time-
dependent Hartree-Fock (TDHF), time-dependent Density Functional Theory
(TDDFT), coupled cluster response, and relativistic TDHF, TDDFT calculations
on the Verdet constant and NSOR. In subsequent work, Cadène et al. 34
investigated the Verdet constant of 129Xe in both gas-phase experiments and
calculations derived from non-relativistic coupled cluster quadratic response
calculations (QR-CC), in which the relativistic effects were approximately
accounted for by employing relativistic effective core potentials (ECPs). With
our implementation, we shall complement these studies and in particular
investigate the relative importance of relativistic (scalar and spin-orbit
coupling) effects and electron correlation to these properties.

The characterization of Two-Photon Absorption (TPA) cross-sections, which can
be related to quadratic response theory, has also gained considerable
attention in different domains and is the third focus of our applications. TPA
was first predicted, using perturbation theory, by Göppert-Mayer 35 in 1931,
but not observed in experiments until the advent of the lasers that are
capable of delivering sufficiently high intensity. The main feature of TPA is
that it occurs with a probability depending quadratically on the incident
light intensity, which results in the TPA-based techniques offering better
spatial resolution than those based on one-photon absorption (OPA). In
materials science, materials with large TPA cross-sections enable applications
including drug delivery, photodynamic therapy, high-resolution, and optical
storage36. Moreover, TPA spectroscopy is also very useful as a research tool.
Concerning the different selection rules of TPA compared to OPA, TPA can
characterize the excited state in the spectrum in the case of OPA spectrum has
been large dispersions, particularly for the complex molecules, which contain
f electrons37, 38.

TPA is proportional to the imaginary part of the second-order
hyperpolarizability $\gamma$, which requires evaluation of the cubic response
function. However, under resonant conditions, it becomes possible to express
the TPA cross-sections in terms of the two-photon matrix39, which can be
obtained from the quadratic response of the reference state wave function.
With this strategy, the TPA cross-sections have been evaluated in various
standard models in quantum chemistry including Hartree-Fock21, MCSCF21, DFT40,
41, 42, and CC43, 44, 45. Moreover, in the last decades, the resonant
inelastic X-ray Scattering (RIXS)46, 47, a two-photon scattering process in
core excitation region, received considerable attention because of the
corresponding improvements in sensitivity and energy resolution48, 49, 50, 51,
which provides valuable information on the electronic structure of both
occupied and virtual states that are not easily accessible by the traditional
spectroscopies. Several approaches aimed at the description of RIXS spectra
for molecular systems based on non-relativistic or approximate relativistic
Hamiltonians have been proposed including algebraic diagrammatic construction
(ADC)52, MCSCF53, DFT54, and EOM-CC55, 56, 57, 58, 59, 60. However, in the
relativistic quantum chemistry field, the implementations of TPA cross-
sections are still scarce, owing to the additional complexity. An
implementation in the DIRAC program by Henriksson et al. 61 enabled pioneering
calculations of TPA cross-sections from the four-component Hartree-Fock
quadratic response theory. In this manuscript, we will focus on TPA for
valence processes and will investigate processes involving core electrons such
as RIXS in a subsequent publication.

Finally, we pay attention to methods that can lower computational costs. This
is of practical importance here since we utilize uncontracted basis sets with
adding many diffuse functions, which generate a large virtual orbital space in
CC calculations. The simplest and most often used method is the utilization of
the MP2 frozen natural orbitals (FNOs)62, 63, 64. While some authors have
pointed out the shortcomings of MP2FNOs for the calculation of linear response
properties65, 64, Surjuse et al. 66 recently suggested using MP2FNOs in EOM-CC
calculations can bring about reduce computational cost while retaining
sufficient accuracy for ionization energies. On the other hand, to the best of
our knowledge, there is no reference yet reporting the performance of MP2FNOs
on TPA calculations.

This manuscript is organized as follows: In Sec. 2, the EOM-CC quadratic
response theory and the corresponding two-photon absorption matrix formulation
are summarized. Section 3 is devoted to the details of the computations we
used to test the implementation. The calculations are presented and discussed
in Secs. 4. Finally, a brief summary of our findings is given in Sec. 5.

##  2 Theory

We base the theory on the time-averaged quasi-energy formalism, which has been
summarized in the landmark paper by Christiansen et al. 13. As the significant
part of the formalism to obtain the quadratic response functions is common to
that of linear response functions, and we have recently provided an extensive
discussion of the implementation details for linear response properties67, in
the current manuscript, we only focus on the equations related to quadratic
response.

The EOM-CC quadratic response function is expressed below:

$\begin{split}{}^{EOM}\langle\langle
X;Y,Z\rangle\rangle_{\omega_{Y},\omega_{Z}}=&\frac{1}{2}C^{\omega}P^{X,Y,Z}\\\
&[-\bar{\mathbb{t}}^{X}(\omega_{X})\mathbb{t}^{Y}(\omega_{Y})\bar{\mathbb{t}}^{0}\boldsymbol{\xi}^{Z}\\\
&+\mathbb{\bar{t}}^{X}(\omega_{X})^{EOM}\mathbb{A}^{Y}\mathbb{t}^{Z}(\omega_{Z})\\\
&-\bar{\mathbb{t}}^{0}\mathbb{t}^{Y}(\omega_{Y})\bar{\mathbb{t}}^{Z}(\omega_{Z})\boldsymbol{\xi}^{X}]\end{split}$
(1)

in terms of the EOM-CC property Jacobian matrix ${}^{EOM}\mathbb{A}^{X}$,

${}^{EOM}\mathbb{A}_{\mu\nu}^{X}=\bra{\mu}\Big{[}\bar{X},\ket{\nu}\bra{HF}\Big{]}\ket{HF}$
(2)

and similarity-transformed one-body property operators

$\bar{X}=e^{-\hat{T}}\hat{X}e^{\hat{T}}$ (3)

where $\mu$ and $\nu$ indicate excited Slater determinants (comprising single
and double excitations for the CCSD model).

The similarity-transformed CC Hamiltonian $\bar{\mathbb{H}}$ is not Hermitian,
which makes the left response amplitudes not just the complex conjugate of
their right counterparts. According to the 2n+1 and 2n+2 rules in perturbation
theory13, for obtaining the quadratic response, it is necessary to solve both
the left and right first-order response equations, given respectively by:

$(\bar{\mathbb{H}}-\omega_{X}\mathbb{I})\mathbb{t}^{X}=-\boldsymbol{\xi}^{X}$
(4)

and

${}^{EOM}\bar{\mathbb{t}}^{X}(\bar{\mathbb{H}}+\omega_{X}\mathbb{I})=-\boldsymbol{\eta}^{X}-\bar{\mathbb{t}}^{0}_{D}\boldsymbol{\xi}_{S}^{X}+(\bar{\mathbb{t}}^{0}\boldsymbol{\xi}^{X})\bar{\mathbb{t}}^{0}.$
(5)

The detailed working equations for the matrix elements of the different terms
in Eqs. 4 and 5 are given in our previous linear response work67, including
those for $\boldsymbol{\sigma}$ vectors and property gradients
$\boldsymbol{\xi}^{X}$. The working equations for new terms appearing in the
quadratic response functions, such as the ${}^{EOM}\mathbb{A}^{X}$ matrix are
presented in the supplementary information.

To define the TPA cross-section, we start from the sum-over states expression
for the two-photon transition matrix elements between the reference state
$\ket{0}$ and the target excited state $\ket{f}$52, 55:

$\begin{split}T_{XY}^{f0}(\omega)=\sum_{n}\left[\frac{\bra{f}\hat{X}\ket{n}\bra{n}\hat{Y}\ket{0}}{\omega_{n}-{(\omega+i\gamma})}+\frac{\bra{f}\hat{Y}\ket{n}\bra{n}\hat{X}\ket{0}}{\omega_{n}-(\omega^{\prime}+i\gamma)}\right]\end{split}$
(6)

where $\gamma$ is the damping factor representing the inverse lifetime. The
frequencies $\omega$ and $\omega^{\prime}$ represent the two external photons,
while $\omega_{f}$ corresponds to the excitation energy between reference
state $\ket{0}$ and the final excited state $\ket{f}$. For TPA the relation:

$\omega+\omega^{\prime}+\omega_{f}=0$ (7)

should be satisfied which means that for a given final state there is only one
independent variable. Within a response formulation, the EOM-CC right and left
frequency-dependent transition moments can be written as 55

$\begin{split}\mathrm{\textbf{Right:}}\qquad^{EOM}T_{XY}^{f0}(\omega)=&-\mathbb{L}_{f}[^{EOM}\mathbb{A}^{X}\mathbb{t}^{Y}(\omega+i\gamma)+^{EOM}\mathbb{A}^{Y}\mathbb{t}^{X}(\omega^{\prime}-i\gamma)\\\
&-(\bar{\mathbb{t}}_{0}\boldsymbol{\xi}^{X})\mathbb{t}^{Y}(\omega+i\gamma)-(\bar{\mathbb{t}}^{0}\boldsymbol{\xi}^{Y})\mathbb{t}^{X}(\omega^{\prime}-i\gamma)\\\
&-(\bar{\mathbb{t}}^{0}\mathbb{t}^{Y}(\omega+i\gamma))\boldsymbol{\xi}^{X}-(\bar{\mathbb{t}}^{0}\mathbb{t}^{X}(\omega^{\prime}-i\gamma))]\end{split}$
(8)
$\begin{split}\mathrm{\textbf{Left:}}\quad^{EOM}T_{XY}^{0f}(\omega)=&-[^{EOM}\bar{\mathbb{t}}^{X}(-\omega^{\prime}-i\gamma)^{EOM}\mathbb{A}^{Y}+^{EOM}\bar{\mathbb{t}}^{Y}(-\omega+i\gamma)^{EOM}\mathbb{A}^{X}\\\
&-(\bar{\mathbb{t}}^{0}\boldsymbol{\xi}^{X})^{EOM}\bar{\mathbb{t}}^{Y}(-\omega+i\gamma)-(\bar{\mathbb{t}}_{0}\boldsymbol{\xi}^{Y})^{EOM}\bar{\mathbb{t}}^{X}(-\omega^{\prime}-i\gamma)]\mathbb{R}_{f}\\\
&+(\bar{\mathbb{t}}_{0}\mathbb{R}_{f})[^{EOM}\bar{\mathbb{t}}^{Y}(-\omega+i\gamma)\boldsymbol{\xi}^{X}+^{EOM}\bar{\mathbb{t}}(-\omega^{\prime}-i\gamma)\boldsymbol{\xi}^{Y}]\end{split}$
(9)

where $\mathbb{R}_{f}$ and $\mathbb{L}_{f}$ are right and left target excited
states, respectively, obtained by solving EOM excitation energy (EOM-EE)
equations. With these left and right transition moments available, the total
scattering amplitudes can then be evaluated by the equation 13, 55:

$S_{XY,ZU}=T_{XY}^{0f}(\omega)T_{ZU}^{f0}(\omega)=\frac{1}{2}[T_{XY}^{0f}(\omega)T_{ZU}^{f0}(\omega)+(T_{ZU}^{0f}(\omega)T_{XY}^{f0}(\omega))^{*}]$
(10)

Finally, the TPA cross-section, $\delta_{TPA}$, is determined by the
components of scattering amplitudes matrix $\mathbf{S}$68:

$\delta_{TPA}=\frac{1}{15}\\{F\sum_{X,Y}S_{XX,YY}+G\sum_{X,Y}S_{XY,XY}+H\sum_{X,Y}S_{XY,YX}\\}$
(11)

The constants $F$, $G$, and $H$ depend on the polarization of the incident
light. In this work, $F=G=H=1$ is selected to represent parallel linearly
polarized light. Moreover, we set up the frequency of the external field as
half of the excitation energy of the target state
($\omega^{\prime}=\omega=-\omega_{f}/2$).

##  3 Computational details

All EOM-CC quadratic response and two-photon absorption calculations were
carried out with development versions (revisions 0923e70dd0, fe4351caf9 and
380df6b) of the DIRAC code69, 70, employing the uncontracted triply-augmented
valence triple zeta Dyall basis set (defined as t-aug-dyall.v3z in inputs) for
heavy elements (In, I, At, Xe, Rn) 71, 72, and an equivalent triply-augmented
uncontracted Dunning basis set (defined as t-aug-cc-pVTZ in inputs) for light
elements (H, F, Cl, Ga, Br) 73, 74, 75. We utilized the exact two-component
(X2C)76 relativistic Hamiltonian, and in some cases, to show the effect of
relativity explicitly, we also provide results using the non-relativistic
Hamiltonian77, 78(as activated by the .Levy-Leblond keyword). To study the
effect of electron correlation, we performed quadratic-response and two-photon
absorption calculations based on mean-field methods such as Hartree-Fock (HF)
and density-functional theory (employing the B3LYP79 density functional
approximation). The relativistic and non-relativistic calculations have been
carried out with the Gaussian type80 and point charge nucleus model,
respectively.

In what follows, we shall use the term orbital as shorthand for both spinors
and spin-orbitals, depending on the Hamiltonian used in the calculation.

In our calculations for heavy elements (HI, HAt, Xe, and Rn), we have profited
from the components of an ongoing implementation in ExaCorr of the Cholesky-
decomposition approach81, 82, 83 to reduce the memory footprint of our
calculations in the step to transform two-electron integrals from AO to MO
basis, with thresholds of 10-9 (Xe and Rn), and 10-4 (HI and HAt), the latter
is looser than the one employed in our previous work; we have carried out
benchmark calculations on selected systems to verify this change did not
significantly alter our results.

The molecular structures employed in all calculations have been taken from the
literature: from Huber 84 for HX (X=F, Cl, Br, I), and from Gomes and Visscher
85 for HAt. The internuclear distances employed are thus H-F (0.91680 Å), H-Cl
(1.27455 Å), H-Br (1.41443 Å), H-I (1.60916 Å), and H–At (1.722 Å).

In the calculations, the size of the correlated virtual spaces in the coupled
cluster is truncated by discarding orbitals with energies above 5 a.u. For the
occupied orbitals, we correlate only valence electrons.

##  4 Results and discussion

###  4.1 First hyperpolarizability of HX(X=F, Cl, Br, I, At)

To demonstrate our implementation we first apply it to calculate the parallel
component of the static first hyperpolarizability ($\beta_{||}$)30 of the
hydrogen halide molecules.

$\beta_{||}=\frac{1}{5}\sum_{i=x,y,z}(\beta_{iiz}+\beta_{izi}+\beta_{zii})$
(12)

Each component is defined by the equation16:

$\beta_{ijk}(-\omega_{\sigma};\omega_{1},\omega_{2})=\sum
P_{-\sigma,1,2}\sum_{n,m}\frac{\bra{0}\hat{\mu}_{i}\ket{n}\bra{n}\hat{\mu}_{j}\ket{m}\bra{m}\hat{\mu}_{k}\ket{0}}{(\omega_{n0}-\omega_{\sigma})(\omega_{m0}-\omega_{2})}$
(13)

where $\hat{\mu}_{i}$ are Cartesian components of the electric dipole
operators, and $\sum P_{-\sigma,1,2}$ indicates the sum of six terms by
permuting the pairs ($i,-\omega_{\sigma}$), ($j,\omega_{1}$),
($k,\omega_{2}$).

Before proceeding with the calculation, it’s crucial to select an appropriate
basis set and establish the correlation space. Our study evaluates the impact
of the basis set and correlation space on the hyperpolarizability of HF
molecules. The results are presented in Table 1 where they are compared with
results from the DALTON program86, 23 and experimental data.

An analysis of the first three rows reveals that both diffuse functions and
polarization functions significantly influence the calculation of
hyperpolarizability, as is well-known in the literature23, 87. For example,
when utilizing the doubly-augmented d-aug-cc-pVDZ basis set, the result is
only 58% of the value obtained with the augmented s-aug-cc-pVDZ basis set.
Conversely, the effect of the correlation space is relatively minor. By
comparing the results of the third and fourth rows, it’s evident that
correlating all virtual orbitals enhances the value by merely around 1%.

In the fourth row, we observe that our calculation, when using the non-
relativistic Hamiltonian, matches the DALTON value (-7.3385 a.u.) precisely
and this serves as a validation of our implementation. Furthermore, based on
the DALTON results, the disparity between EOM-QRCC and QR-CC is approximately
4.5%. This deviation stems from the absence of size extensivity in the
transition moments of the EOM model. This inconsistency between EOM and CC was
previously highlighted in research on linear response properties (see Yuan et
al. 67 and references therein), and we plan to delve deeper into this topic by
studying a wider array of molecules for both linear and quadratic response
properties in follow-up work.

The DALTON QR-CC results reveal a discrepancy of 3 a.u. between the calculated
(-7.94 a.u.) and the experimental value (-10.88 a.u.). Beyond the limitations
of the basis set, which might benefit from the inclusion of more diffuse and
polarization functions, it is essential not to overlook the vibrational
effects. These effects, amounting to 1.2 a.u., were discussed at the Hartree-
Fock level by Bishop and Norman 88.

Comparison between the fourth, fifth, and sixth rows of the QR-CC calculations
reveals that adding more diffuse functions (from d-aug-cc-pVTZ to t-aug-cc-
pVTZ) improves accuracy more significantly than incorporating additional
polarization functions (from d-aug-cc-pVTZ to d-aug-cc-pVQZ). Given that the
QZ calculations are notably more resource-intensive than TZ ones, we will
employ the t-aug-cc-pVTZ basis set for the subsequent calculations on heavier
elements.

Table 1: Benchmark of basis sets and correlation virtual orbital space for QR-
CC calculations of the static $\beta_{||}$ (a.u.) of the HF molecule

Basis |  NR-EOM |  NR-EOMa |  NR-QR-CCa |  Exp89
---|---|---|---|---
s-aug-ccpVDZb |  -9.4232 |  |  |
d-aug-ccpVDZb |  -5.5463 |  |  |
d-aug-ccpVTZb |  -7.2677 |  |  |
d-aug-ccpVTZc |  -7.3385 |  -7.3385 |  -7.6718 |
d-aug-ccpVQZc |  |  |  -7.7521 |
t-aug-ccpVTZc |  |  |  -7.9385 |
|  |  |  |  -10.88

  * a

Calculations were performed using the DALTON program

  * b

Truncating the virtual orbital space at 5 a.u.

  * c

Correlating all virtual orbitals

In Table 2, the static hyperpolarizability of hydrogen halides molecules (from
F to At) is displayed, in which we show the Hartree-Fock, B3LYP, and EOM-CC
results for both the non-relativistic and the X2C Hamiltonian.

At Hartree-Fock level, $\beta_{zxx}$, $\beta_{zzz}$, and $\beta_{||}$ all
generally exhibit an upward trend from HF to HAt in both relativistic and non-
relativistic calculations. This pattern is also discernible in the correlated
calculations, though the exact values vary slightly. Both CC and B3LYP results
indicate that electron correlation tends to decrease the value of
$\beta_{zxx}$ for all species on the series. Nevertheless, for the
$\beta_{zzz}$, CC results indicate an increase in value due to electron
correlation for all molecules except HF, whereas B3LYP shows the opposite
pattern. This divergence between CC and B3LYP leads to discrepancies in the
final $\beta_{||}$ value. For the heavier molecules, B3LYP values deviate
considerably from both the Hartree-Fock and the CC ones.

Accounting for relativistic effects is, as expected, absolutely essential for
systems containing heavier elements. For instance, both Hartree-Fock and CC
models reveal that for HAt, the non-relativistic outcomes are approximately
half of their relativistic counterparts. With the exception of HF, all three
models consistently suggest that relativistic effects augment the $\beta_{||}$
values for all molecules. Furthermore, the effects of relativity on
$\beta_{zxx}$ are much larger than that of $\beta_{zzz}$.

To understand the different magnitude of relativistic correction for
$\beta_{zxx}$ and $\beta_{zzz}$, we can look at equation 13. Through spin-
orbit coupling (SOC), relativity modifies the response function by shifting
the location of the poles of the response functions and by introducing
additional transition channels. Without SOC, the ground state is in all cases
of pure ${}^{1}\Sigma_{0}^{+}$ symmetry. When introducing SOC, the
${}^{1}\Sigma$ designation is no longer strictly valid, and allowed
transitions are only characterized by the remaining quantum numbers, thereby
yielding $0^{+}\to 0^{+}$, and $0^{+}\to 1$ transitions. The former is
evidently connected to the $z$-component of the transition dipole moment,
while the latter corresponds to the $x$ and $y$-components. Specifically, the
$\beta_{zzz}$ component only permits contributions from transitions to $0^{+}$
states, originating from ${}^{3}\Pi_{0+}$, and ${}^{3}\Sigma_{0+}$ states. On
the other hand, for the $\beta_{zxx}$ component, transitions to the $1$ state,
potentially emerging from ${}^{3}\Pi_{1}$, ${}^{3}\Sigma_{1}$ and
${}^{1}\Pi_{1}$ states, are also allowed, with an enhanced contribution from
the singlet states ${}^{1}\Pi$.

Table 2: Static hyperpolarizability (a.u.) of hydrogen halides HX(X=F, Cl, Br,
I, At)

| HFa | HFb | B3LYPa | B3LYPb | CCa | CCb
---|---|---|---|---|---|---
$\beta_{zxx}$
HF | -0.5091 | -0.5100 | -1.4916 | -1.5565 | -1.5276 | -1.5286
HCl | 2.3017 | 2.3672 | -0.1127 | 0.0117 | -0.0858 | -0.0157
HBr | 5.2994 | 5.9155 | 2.8120 | 3.8513 | 2.6604 | 3.4213
HI | 10.7728 | 13.6493 | 4.8719 | 10.2006 | 6.7316 | 9.8900
HAt | 16.0854 | 38.3140 | 8.7719 | 42.5584 | 11.1347 | 32.2297
$\beta_{zzz}$
HF | -8.3950 | -8.4157 | -9.8253 | -9.8476 | -9.4973 | -9.4983
HCl | -11.4505 | -11.4435 | -13.1441 | -13.0321 | -10.2439 | -10.1522
HBr | -11.0481 | -11.1427 | -11.4062 | -10.8785 | -6.4784 | -5.9790
HI | -2.9448 | -3.9634 | -6.0444 | -5.3245 | 5.1093 | 4.9752
HAt | 5.3664 | 5.5476 | 1.3711 | 5.2236 | 16.9624 | 11.4539
$\beta_{||}$
HF | -5.6479 | -5.6614 | -7.6852 | -7.7763 | -7.5315 | -7.5333
HCl | -4.1082 | -4.0255 | -8.0218 | -7.8052 | -6.2493 | -6.1101
HBr | -0.2696 | 0.4130 | -3.4693 | -1.9055 | -0.6946 | 0.5181
HI | 11.1606 | 14.0011 | 2.2196 | 9.0460 | 11.1434 | 14.8531
HAt | 22.5224 | 49.3054 | 11.3490 | 54.2043 | 23.5391 | 45.5480

  * a

non-relativistic calculation using the Levy-Leblond Hamiltonian

  * b

Relativistic calculation using the X2C Hamiltonian

We now turn to the frequency-dependence of the hyperpolarizability, focusing
on the impacts of relativistic and electron correlation effects. Using both CC
and Hartree-Fock methods, we assess the hyperpolarizability of the single-
frequency optical processes -the second harmonic generation (SHG) of HI, in
which we have $\omega_{\sigma}=2\omega_{1}=2\omega_{2}=2\omega$. Figure 1
displays the result of frequency ranging from 0.0 to 0.115 a.u. To further
interpret these curves, we also calculate the excitation energy for the lowest
five electronic states by diagonalizing $\bar{\mathbb{H}}$. These results are
summarized in Table 3.

Overall the dispersion curves for CC and HF are qualitatively similar for both
non-relativistic and relativistic calculations, though the values for CC
values are larger than the HF ones, which is due to larger excitation energy
of the lowest lying dipole allowed states, that can be observed in Table 3.

One can clearly find singularities in the relativistic results, located at
0.0958 a.u for HF, and located at 0.1000 a.u. for CC. According to equation
13, in SHG, a singularity should appear twice on the curve: the first pole
corresponds to half of the excitation energy, while the second aligns with the
full excitation energy. Observing this pattern, we pinpoint the singularities
in our curves to the first pole associated with the $a^{3}\Pi_{0+}$ state
(with excitation energies for HF and CC being 0.1916 a.u. and 0.2001 a.u.,
respectively).

In the absence of SOC, the transition to $a^{3}\Pi_{0+}$ is spin-forbidden.
This observation aligns with the singularities appearing exclusively in the
relativistic calculations. On the other hand, the pole related to the
transition to $a^{3}\Pi_{1}$, is not observed, despite its permissibility with
SOC. This is attributed to our focus on the $\beta_{zxx}$ and $\beta_{zzz}$
components. One can find that in equation 13, the transition dipole moment
$\bra{0}\hat{\mu}_{z}\ket{n}$ is zero for all $\ket{n}$=$\ket{1}$ states.

Figure 1: Frequency-dependent hyperpolarizability of HI (Left: $\beta_{zxx}$, Right: $\beta_{zzz}$). The red and black vertical lines are half of the excitation energy of $a^{3}\Pi_{0+}$ state for HF(0.0958 a.u.) and EOM-CC (0.1001 a.u.), respectively.  Table 3: Excitation energy (a.u.) of HI for the lowest five states State |  X2C-EOM-CC |  X2C-HF |  NR-EOM-CC |  NR-HF
---|---|---|---|---
$a^{3}\Pi_{2}$ |  0.1773 |  0.1670 |  0.1894 |  0.1775
$a^{3}\Pi_{1}$ |  0.1831 |  0.1743 |  0.1894 |  0.1775
$a^{3}\Pi_{0-}$ |  0.1979 |  0.1861 |  0.1894 |  0.1775
$a^{3}\Pi_{0+}$ |  0.2001 |  0.1916 |  0.1894 |  0.1775
$A^{1}\Pi_{1}$ |  0.2108 |  0.2120 |  0.2111 |  0.2189

###  4.2 Verdet constant of Xe and Rn

In the current section, we show a calculation of the Verdet constant as an
illustrative example of the use of our implementation for a mixed electric-
magnetic property. The Verdet constant is evaluated with the following
frequency-dependent quadratic response function90, 27:

$V(\omega)=\omega\frac{eN\epsilon_{xyz}}{24c_{0}\epsilon_{0}m_{e}}\text{Im}\langle\langle\hat{\mu}_{x};\hat{\mu}_{y},\hat{m}_{z}\rangle\rangle_{\omega,0}$
(14)

with $N$ the number density of the gas, $e$ the elementary charge, $m_{e}$ the
electron mass, $c_{0}$ the speed of light in vacuo, and $\hat{m}_{z}$ is
magnetic dipole moment operator.

We calculate the Verdet constant at three different laser wavelengths (589 nm,
694.3 nm, and 1064 nm) for Xe and Rn and list the results in Table 4. For Xe,
compared to the experiment, the HF value shows sizeable relative errors of
about 10%. The relativistic effect increases the value and reduces the error
to 5% at the HF level. It’s evident that the scalar relativistic results (with
the SFDC Hamiltonian) closely align with the NR-HF values, but deviate
significantly from the DC and X2C results. One can note the scalar
relativistic effects decrease Verdet’s constant value while SOC moves the
results in the opposite direction, but more strongly. This suggests that a
major portion of the relativistic correction originates from the spin-orbit
coupling, and considering only scalar relativistic effects may lead to an
underestimation of results. Additionally, upon investigating the influence of
the Gaunt interaction, we determine it to further increase the value of the
Verdet constants, but much more modestly (0.5%).

The effect of electron correlation is also to increase the Verdet constant,
but the higher the degree of electron correlation recovered, the less
important the increase. If we compare the DALTON NR-CCSD results obtained with
a truncated correlation in space with those in which all occupied and virtual
orbitals are included in the calculation, we observe a 2% difference, with the
latter calculation showing smaller values. On the other hand, comparing CCSD
and CC391 results with the truncated correlation space, we observe a small
increase in the Verdet constant for CC3. Relative to the QR-CC results, our
non-relativistic EOM calculation seems to overestimate the correlation effect
by about 1%, due to the non-extensivity issue discussed for the
hyperpolarizabilities (and in ref 67). Even though this overestimation causes
the X2C-EOM value to be significantly larger than experimental results, we
anticipate that, given the downward trend in the QR-CC results upon improving
the correlation space discussed above (approximately -0.08 (10-3 rad/(T m) at
the 1064 nm wavelength), the X2C-EOM value with a complete orbital space are
expected to come closer to the experimental values.

When we examine the Rn, we find the relativistic effect to be substantial for
both HF and CC. For example, at wavelength 589 nm, the non-relativistic EOM
value is 18.90 (10-3 rad/(T m)), but the relativistic EOM value is 25.91 (10-3
rad/(T m)). Even in the absence of experimental data for reference, such a
pronounced correction underscores the importance of accounting for
relativistic effects. Beyond this amplified relativistic effect, other
observations for Rn align with those for Xe. This includes the dominance of
the relativistic effect by spin-orbit coupling, the marginal influence of the
Gaunt interaction, and a comparable magnitude of difference between EOM and
QR-CC.

It’s also worth noting, as reported in the supplementary material of Ref 33,
that the performance of B3LYP is somewhat poor. It tends to overestimate the
values in both non-relativistic and relativistic calculations, with errors
approaching 18% for Xe compared to the experiment.

Method |  $\lambda$(589 nm) |  $\lambda$(694.3 nm) |  $\lambda$(1064 nm)
---|---|---|---
Xe
NR-HF | 11.18 | 7.88 | 3.26
SFDC-HFd | 11.12 | 7.83 | 3.24
X2C-HF | 11.61 | 8.17 | 3.37
DC-HFb | 11.61 | 8.17 | 3.37
DCG-HFc | 11.66 | 8.20 | 3.39
NR-B3LYP | 13.62 | 9.56 | 3.93
X2C-B3LYP | 14.60 | 10.21 | 4.18
NR-EOM | 12.55 | 8.83 | 3.65
X2C-EOM | 13.10 | 9.21 | 3.79
NR-CCSDa | 12.39 | 8.72 | 3.60
NR-CCSDa * | 12.11 | 8.53 | 3.52
NR-CC3a | 12.46 | 8.77 | 3.62
Exp | 12.30e |  | 3.56$\pm$0.10f
Rn
NR-HF | 16.80 | 11.80 | 4.86
SFDC-HFd | 16.56 | 11.60 | 4.76
X2C-HF | 23.13 | 16.02 | 6.48
DC-HFb | 23.14 | 16.03 | 6.48
DCG-HFc | 23.25 | 16.10 | 6.51
NR-B3LYP | 20.08 | 14.04 | 5.74
X2C-B3LYP | 30.38 | 20.83 | 8.31
NR-EOM | 18.90 | 13.26 | 5.44
X2C-EOM | 25.91 | 17.92 | 7.22
NR-CCSDa | 18.65 | 13.08 | 5.37

  * a

Calculations were performed using the DALTON program

  * *

Include all occupied and virtual orbitals

  * b

Dirac-Coulomb Hamiltonian92

  * c

Dirac-Coulomb plus Gaunt interaction Hamiltonian

  * d

Dirac-Coulomb without spin-orbit coupling Hamiltonian93

  * e

Reference94

  * f

Reference34

Table 4: Verdet Constant V($\omega$) [in 10-3 rad/(T m)] for Gaseous Xe and Rn
at Different Laser Wavelengths

###  4.3 Two-photon absorption cross-sections

Finally, we consider the two-photon absorption cross-sections. As a showcase
for our implementation, we highlight the two-photon transitions for group IIIB
divalent ions, namely, (Ga+, and In+). Such systems have been discussed in the
literature95, 96 for their potential use as an atomic clock. As an initial
step in this exploration, we focus on spin-allowed transitions, allowing for
comparative analyses with other non-relativistic programs. We intend to
address spin-forbidden transitions in subsequent work.

In Table 5, we calculate the TPA cross-sections and the corresponding
excitation energy of the target states for Ga+, and In+ by HF and CC. Both
methods indicate that relativistic effects reduce TPA cross-sections. Notably
for In+, the X2C value is roughly 60% of its non-relativistic counterpart.
When comparing HF and CC, we find that electron correlation further diminishes
the TPA cross-sections. For Ga+, the electron correlation effect is slightly
larger than that of relativity, whereas, in the case of In+, relativistic
effects show larger contributions to the final TPA cross-sections than
electron correlation.

The observed trend, wherein both electron correlation and relativistic effects
reduce the TPA, can be understood by examining their effect on the excitation
energies. We find that both factors increase the excitation energy of the
target state ($(4s5s)^{1}S_{0}$ for Ga+, $(5s6s)^{1}S_{0}$ for In+) and from
equation 6, a higher excitation energy in the denominator usually results in a
smaller transition amplitude.

In non-relativistic calculation, our EOM results closely match the TPA results
derived from QR-CC in DALTON, with only about a 3% discrepancy, in line with
the discrepancies observed for the hyperpolarizabilties.

Systems(transitions) | NR-HF | X2C-HF | DC-HF | NR-EOMa | X2C-EOMa | NR-CCb
---|---|---|---|---|---|---
Two-photon absorption cross-sections
Ga+(4s-5s) | 3211.01 | 2710.27 | 2704.53 | 2535.17 | 2128.33 | 2455.97
In+(5s-6s) | 7964.78 | 4851.59 | 4831.34 | 6203.59 | 3745.02 | 6022.11
Excitation energy
Ga+(4s-5s) | 0.4382 | 0.4472 | 0.4473 | 0.4718 | 0.4812 | 0.4718
In+(5s-6s) | 0.3760 | 0.3983 | 0.3986 | 0.4042 | 0.4269 | 0.4042

  * a

Correlate both (n-1)d and (n)s shell of total 12 electrons

  * b

Calculations were performed using the DALTON program

Table 5: Two-photon absorption cross-sections $\delta$ (a.u.) and excitation
energy (a.u.) of the target states for Ga+, and In+

We test MP2FNO-based EOM-CC energy and TPA for Ga+ and display the results in
Table 6. The relativistic MP2FNOs are generated based upon our previous
implementation in DIRAC97. From the first three columns, it’s evident that
when using standard selection schemes, we cannot achieve reasonable TPA values
in comparison to the original canonical orbital results, even with a threshold
of 1.0d-6 for the occupation number where we retrieve more than 99.99%
correlation energy for the ground state. We find the excitation energy of the
TPA target state (4s5s)${}^{1}S_{0}$ is markedly overestimated in FNO
calculations, especially when using threshold of 1.0d-4 and 1.0d-5.

At the MP2 level, the 5s orbital has a small contribution to the correlation
energy in the ground state. As a result, when we obtain the MP2 density matrix
and natural orbitals, the occupation numbers for natural orbitals primarily
influenced by 5s orbitals are exceedingly small. These orbitals are therefore
omitted by the selection scheme based purely on the threshold of occupation
numbers. While the 5s orbital is not be particularly important for the ground
state, it plays a significant role in the excited state under consideration.

Recognizing this, we’ve adjusted the selection scheme in the MP2FNO
implementation. Besides selecting natural orbitals with occupation numbers
exceeding the threshold, we also incorporate all doubly-degenerate orbitals
with occupation numbers below the threshold (in atomic systems, that
corresponds to s1/2 and p1/2 orbitals). We provide a more in-depth discussion
on this point in the supplementary information, where we include the
excitation energy for the eight lowest states and the energy of virtual
orbitals.

With this expanded natural orbital space, we revisit the CC and TPA
calculations, presenting the outcomes in the 4th to 6th columns of Table 6. We
note a marked improvement in the excitation energy of the target state. Even
at a threshold of 1.0d-4, the discrepancy when compared to full canonical
results is around 1%. For the TPA, we also see more accurate results. For
example, at a threshold of 1.0d-6, the error drops from 57% to 30% by
correlating only two additional doubly-degenerate orbitals. We also detect a
consistent trend of approaching the canonical orbitals results when going from
1.0d-4 to 1.0d-6.

There remains a 30% discrepancy between truncated FNOs and canonical orbital
results. From equation 6, to achieve precise scattering amplitudes, we cannot
overlook the impact on the transition dipole moment, which is highly sensitive
to diffuse orbitals. Conversely, these diffuse orbitals will have low
occupation numbers in the MP2 calculations for the ground state and will
therefore be removed from the correlating space, even with the slightly
modified procedure we used.

FNOa | FNOb | FNOc | FNOd | FNOe | FNOf | Canonical | Expg
---|---|---|---|---|---|---|---
Number of correlated virtual orbitals
25 | 35 | 48 | 29 | 38 | 50 | 82 |
CCSD ground state correlation energy
-0.2371 | -0.2396 | -0.2399 | -0.2373 | -0.2397 | -0.2399 | -0.2399 |
Excitation energy of the target state (4s5s)${}^{1}S$
0.5864 | 0.6414 | 0.5192 | 0.4861 | 0.4838 | 0.4835 | 0.4812 | 0.4860
Two-photon absorption cross-sections
0.19976 |  $<$1.0d-15 | 3353.67 | 200.06 | 1300.01 | 1491.51 | 2128.33 |

  * a

FNOs with the threshold of occupation number 1.0d-4

  * b

FNOs with the threshold of occupation number 1.0d-5

  * c

FNOs with the threshold of occupation number 1.0d-6

  * d

FNOs with the threshold of occupation number 1.0d-4 plus doubly-degenerate
orbitals

  * e

FNOs with the threshold of occupation number 1.0d-5 plus doubly-degenerate
orbitals

  * f

FNOs with the threshold of occupation number 1.0d-6 plus doubly-degenerate
orbitals

  * g

Results from NIST

Table 6: Performance of MP2FNOs on Two-photon absorption cross-sections
$\delta$ (a.u.) and excitation energy (a.u.) of the target states for Ga+

##  5 Conclusion

In this work, we implement the relativistic Equation-of-Motion Coupled Cluster
method to study the molecular quadratic response properties and two-photon
absorption cross-sections. This implementation is accomplished in the GPU-
accelerated coupled cluster module of DIRAC (ExaCorr), extending our previous
linear response coupled cluster code 67 to solve both left and right response
equations.

We have validated the implementation by assessing the purely electric
properties such as static and frequency-dependent first hyperpolarizability of
six hydrogen halide molecules from HF to HAt. Using a non-relativistic
Hamiltonian, our code exactly reproduces the EOM-based quadratic response
properties implemented in the DALTON code. Compared to Hartree-Fock and B3LYP
linear response calculations, our relativistic EOM quadratic response
calculation shows the significance of both relativistic effect and electron
correlation.

We have also investigated the Verdet constant, a mixed electric-magnetic
property, for Xe and Rn with different Hamiltonian and correlation models.
Both correlation and spin-orbit coupling are found to have pronounced effects.
Compared to NR-QR-CC, our NR-EOM calculation overestimated the results by
roughly 1%. While the X2C-EOM calculation deviates from the experimental value
more than its non-relativistic counterpart, we find this to be due to error
cancellation in the treatment of electron correlation, and we estimate that
using larger correlating spaces should bring our X2C results more in line with
experiment. We note such calculations are currently not feasible with our
single-node code due to memory limits in computational resources at our
disposal.

On the other hand, as consistent with previous works, we also observe the
performance of B3LYP on the Verdet constant is poor with an error of 18%
compared to the experiment for Xe. Such large deviations suggest it is
important to investigate these properties with the more reliable coupled
cluster method.

At last, in our study of the two-photon absorption in Ga+, and In+, we find
relativistic and electron correlation effects both decrease the corresponding
TPA cross-sections. We analyzed the results by evaluating the excitation
energy of the target state and found that both effects increase the excitation
energies.

It is worth noting that most calculations are limited in size since the
quadratic response properties usually require more diffuse functions, which is
challenging for the memory requirement in the current single-node
implementation. There is an imperative to develop algorithms that can lower
computational costs. In the current work, we utilize the MP2FNOs to reduce the
virtual orbital space in TPA calculations and we find for the low-lying
states, MP2FNO can effectively decrease the calculation cost while maintaining
accuracy.

For higher states such as the target state in the two-photon transition
considered, the bias of MP2FNOs towards the ground state may remove diffuse
orbitals which will be important for excited states. A better way to consider
the influence of these diffuse orbitals is to take account of the excited
state in a more sophisticated manner, such as introducing the corresponding
natural transition orbitals98. Exploring this further is among our future
research objectives. Another natural development is to extend the current code
to use libraries tailored for distributed memory computing architectures, such
as the ExaTENSOR library, something which we are currently pursuing.

This research used resources of the Oak Ridge Leadership Computing Facility,
which is a DOE Office of Science User Facility supported under Contract DE-
AC05-00OR22725 (allocations CHM160, CHM191 and CHP109). XY, LH and ASPG
acknowledge funding from projects Labex CaPPA (Grant No. ANR-11-LABX-0005-01)
and CompRIXS (Grant Nos. ANR-19CE29-0019 and DFG JA 2329/6-1), the I-SITE ULNE
project OVERSEE and MESONM International Associated Laboratory (LAI) (Grant
No. ANR-16-IDEX-0004), as well support from the French national supercomputing
facilities (Grant Nos. DARI A0130801859, A0110801859, project grand challenge
AdAstra GDA2210).

The data (input/output) corresponding to the calculations of this paper are
available at the Zenodo repository under DOI: 10.5281/zenodo.8341323.

## References

  * Cronstrand et al. 2005 Cronstrand, P.; Luo, Y.; Ågren, H. _Advances in Quantum Chemistry_ ; Elsevier, 2005; Vol. 50; pp 1–21
  * Papadopoulos et al. 2006 Papadopoulos, M. G., Sadlej, A. J., Leszczynski, J., Eds. _Non-linear optical properties of matter: from molecules to condensed phases_ ; Challenges and advances in computational chemistry and physics 1; Springer: Dordrecht, 2006; OCLC: 72145327
  * Barron 2009 Barron, L. D. _Molecular light scattering and optical activity_ ; Cambridge University Press, 2009
  * Sénéchal et al. 2004 Sénéchal, K.; Toupet, L.; Ledoux, I.; Zyss, J.; Le Bozec, H.; Maury, O. First lanthanide dipolar complexes for second-order nonlinear optics. _Chem. Commun._ 2004, 2180–2181
  * Law et al. 2010 Law, G.-L.; Wong, K.-L.; Lau, K.-K.; Lap, S.-t.; Tanner, P. A.; Kuo, F.; Wong, W.-T. Nonlinear optical activity in dipolar organic–lanthanide complexes. _J. Mater. Chem._ 2010, _20_ , 4074
  * Andraud and Maury 2009 Andraud, C.; Maury, O. Lanthanide Complexes for Nonlinear Optics: From Fundamental Aspects to Applications. _Eur. J. Inorg. Chem._ 2009, _2009_ , 4357–4371
  * Tancrez et al. 2005 Tancrez, N.; Feuvrie, C.; Ledoux, I.; Zyss, J.; Toupet, L.; Le Bozec, H.; Maury, O. Lanthanide Complexes for Second Order Nonlinear Optics: Evidence for the Direct Contribution of f Electrons to the Quadratic Hyperpolarizability. _J. Am. Chem. Soc._ 2005, _127_ , 13474–13475
  * Sénéchal-David et al. 2006 Sénéchal-David, K.; Hemeryck, A.; Tancrez, N.; Toupet, L.; Williams, J. A. G.; Ledoux, I.; Zyss, J.; Boucekkine, A.; Guégan, J.-P.; Le Bozec, H.; Maury, O. Synthesis, Structural Studies, Theoretical Calculations, and Linear and Nonlinear Optical Properties of Terpyridyl Lanthanide Complexes: New Evidence for the Contribution of f Electrons to the NLO Activity. _J. Am. Chem. Soc._ 2006, _128_ , 12243–12255
  * Valore et al. 2010 Valore, A.; Cariati, E.; Righetto, S.; Roberto, D.; Tessore, F.; Ugo, R.; Fragalà, I. L.; Fragalà, M. E.; Malandrino, G.; De Angelis, F.; Belpassi, L.; Ledoux-Rak, I.; Thi, K. H.; Zyss, J. Fluorinated $\beta$-Diketonate Diglyme Lanthanide Complexes as New Second-Order Nonlinear Optical Chromophores: The Role of f Electrons in the Dipolar and Octupolar Contribution to Quadratic Hyperpolarizability. _J. Am. Chem. Soc._ 2010, _132_ , 4966–4970
  * Klepov et al. 2014 Klepov, V. V.; Serezhkina, L. B.; Vologzhanina, A. V.; Pushkin, D. V.; Sergeeva, O. A.; Stefanovich, S. Y.; Serezhkin, V. N. Tris(acrylato)uranylates as a scaffold for NLO materials. _Inorg. Chem. Commun_ 2014, _46_ , 5–8
  * Wang et al. 2010 Wang, S.; Alekseev, E. V.; Ling, J.; Liu, G.; Depmeier, W.; Albrecht-Schmitt, T. E. Polarity and Chirality in Uranyl Borates: Insights into Understanding the Vitrification of Nuclear Waste and the Development of Nonlinear Optical Materials. _Chem. Mater._ 2010, _22_ , 2155–2163
  * Serezhkin et al. 2017 Serezhkin, V. N.; Grigoriev, M. S.; Abdulmyanov, A. R.; Fedoseev, A. M.; Savchenkov, A. V.; Stefanovich, S. Y.; Serezhkina, L. B. Syntheses, Crystal Structures, and Nonlinear Optical Activity of Cs ${}_{\textrm{2}}$ Ba[AnO ${}_{\textrm{2}}$ (C ${}_{\textrm{2}}$ H ${}_{\textrm{5}}$ COO) ${}_{\textrm{3}}$ ] ${}_{\textrm{4}}$ (An = U, Np, Pu) and Unprecedented Octanuclear Complex Units in KR ${}_{\textrm{2}}$ (H ${}_{\textrm{2}}$ O) ${}_{\textrm{8}}$ [UO ${}_{\textrm{2}}$ (C ${}_{\textrm{2}}$ H ${}_{\textrm{5}}$ COO) ${}_{\textrm{3}}$ ] ${}_{\textrm{5}}$ (R = Sr, Ba). _Inorg. Chem._ 2017, _56_ , 7151–7160
  * Christiansen et al. 1998 Christiansen, O.; Jørgensen, P.; Hättig, C. Response functions from Fourier component variational perturbation theory applied to a time-averaged quasienergy. _Int. J. Quantum Chem._ 1998, _68_ , 1–52
  * Helgaker et al. 2012 Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.; Ruud, K. Recent Advances in Wave Function-Based Methods of Molecular-Property Calculations. _Chem. Rev._ 2012, _112_ , 543–631
  * Norman 2011 Norman, P. A perspective on nonresonant and resonant electronic response theory for time-dependent molecular properties. _Phys. Chem. Chem. Phys._ 2011, _13_ , 20519
  * Norman et al. 2018 Norman, P.; Ruud, K.; Saue, T. _Principles and practices of molecular properties: theory, modeling, and simulations_ ; Wiley: Hoboken, NJ, 2018
  * Sekino and Bartlett 1986 Sekino, H.; Bartlett, R. J. Frequency dependent nonlinear optical properties of molecules. _J. Chem. Phys._ 1986, _85_ , 976–989
  * Rice et al. 1990 Rice, J. E.; Amos, R. D.; Colwell, S. M.; Handy, N. C.; Sanz, J. Frequency dependent hyperpolarizabilities with application to formaldehyde and methyl fluoride. _J. Chem. Phys._ 1990, _93_ , 8828–8839
  * Rice and Handy 1992 Rice, J. E.; Handy, N. C. The calculation of frequency-dependent hyperpolarizabilities including electron correlation effects. _Int. J. Quantum Chem._ 1992, _43_ , 91–118
  * Olsen and Jørgensen 1985 Olsen, J.; Jørgensen, P. Linear and nonlinear response functions for an exact state and for an MCSCF state. _J. Chem. Phys._ 1985, _82_ , 3235–3264
  * Hettema et al. 1992 Hettema, H.; Jensen, H. J. A.; Jørgensen, P.; Olsen, J. Quadratic response functions for a multiconfigurational self‐consistent field wave function. _J. Chem. Phys._ 1992, _97_ , 1174–1190
  * Rozyczko and Bartlett 1997 Rozyczko, P.; Bartlett, R. J. Frequency dependent equation-of-motion coupled cluster hyperpolarizabilities: Resolution of the discrepancy between theory and experiment for HF? _J. Chem. Phys._ 1997, _107_ , 10823–10826
  * Hättig et al. 1997 Hättig, C.; Christiansen, O.; Koch, H.; Jørgensen, P. Frequency-dependent first hyperpolarizabilities using coupled cluster quadratic response theory. _Chem. Phys. Lett_ 1997, _269_ , 428–434
  * Gauss et al. 1998 Gauss, J.; Christiansen, O.; Stanton, J. F. Triple excitation effects in coupled-cluster calculations of frequency-dependent hyperpolarizabilities. _Chem. Phys. Lett_ 1998, _296_ , 117–124
  * Sałek et al. 2002 Sałek, P.; Vahtras, O.; Helgaker, T.; Ågren, H. Density-functional theory of linear and nonlinear time-dependent molecular properties. _J. Chem. Phys._ 2002, _117_ , 9630–9645
  * Schipper et al. 2000 Schipper, P. R. T.; Gritsenko, O. V.; Van Gisbergen, S. J. A.; Baerends, E. J. Molecular calculations of excitation energies and (hyper)polarizabilities with a statistical average of orbital model exchange-correlation potentials. _J. Chem. Phys._ 2000, _112_ , 1344–1352
  * Norman and Jensen 2004 Norman, P.; Jensen, H. J. A. Quadratic response functions in the time-dependent four-component Hartree-Fock approximation. _J. Chem. Phys._ 2004, _121_ , 6145–6154
  * Henriksson et al. 2008 Henriksson, J.; Saue, T.; Norman, P. Quadratic response functions in the relativistic four-component Kohn-Sham approximation. _J. Chem. Phys._ 2008, _128_ , 024105
  * Burland et al. 1994 Burland, D. M.; Miller, R. D.; Walsh, C. A. Second-order nonlinearity in poled-polymer systems. _Chem. Rev._ 1994, _94_ , 31–75
  * Shelton and Rice 1994 Shelton, D. P.; Rice, J. E. Measurements and calculations of the hyperpolarizabilities of atoms and small molecules in the gas phase. _Chem. Rev._ 1994, _94_ , 3–29
  * Datta and Pati 2006 Datta, A.; Pati, S. K. Dipolar interactions and hydrogen bonding in supramolecular aggregates: understanding cooperative phenomena for 1st hyperpolarizability. _Chem. Soc. Rev._ 2006, _35_ , 1305
  * Savukov et al. 2006 Savukov, I. M.; Lee, S.-K.; Romalis, M. V. Optical detection of liquid-state NMR. _Nature_ 2006, _442_ , 1021–1024
  * Ikäläinen et al. 2012 Ikäläinen, S.; Lantto, P.; Vaara, J. Fully Relativistic Calculations of Faraday and Nuclear Spin-Induced Optical Rotation in Xenon. _J. Chem. Theory Comput._ 2012, _8_ , 91–98
  * Cadène et al. 2015 Cadène, A.; Fouché, M.; Rivère, A.; Battesti, R.; Coriani, S.; Rizzo, A.; Rizzo, C. Circular and linear magnetic birefringences in xenon at $\lambda$ = 1064 nm. _J. Chem. Phys._ 2015, _142_ , 124313
  * Göppert-Mayer 1931 Göppert-Mayer, M. Über Elementarakte mit zwei Quantensprüngen. _Ann. Phys._ 1931, _401_ , 273–294
  * Lee et al. 2021 Lee, C.; Xu, E. Z.; Liu, Y.; Teitelboim, A.; Yao, K.; Fernandez-Bravo, A.; Kotulska, A. M.; Nam, S. H.; Suh, Y. D.; Bednarkiewicz, A.; Cohen, B. E.; Chan, E. M.; Schuck, P. J. Giant nonlinear optical responses from photon-avalanching nanoparticles. _Nature_ 2021, _589_ , 230–235
  * Barker et al. 1987 Barker, T. J.; Denning, R. G.; Thorne, J. R. G. Applications of two-photon spectroscopy to inorganic compounds. 1. Spectrum and electronic structure of dicesium tetrachlorodioxouranate. _Inorg. Chem._ 1987, _26_ , 1721–1732
  * Barker et al. 1992 Barker, T. J.; Denning, R. G.; Thorne, J. R. G. Applications of two-photon spectroscopy to inorganic compounds. 2. Spectrum and electronic structure of cesium uranyl nitrate, CsUO2(NO3)3. _Inorg. Chem._ 1992, _31_ , 1344–1353
  * Hättig et al. 1998 Hättig, C.; Christiansen, O.; Jørgensen, P. Multiphoton transition moments and absorption cross sections in coupled cluster response theory employing variational transition moment functionals. _J. Chem. Phys._ 1998, _108_ , 8331–8354
  * Sałek et al. 2003 Sałek, P.; Vahtras, O.; Guo, J.; Luo, Y.; Helgaker, T.; Ågren, H. Calculations of two-photon absorption cross sections by means of density-functional theory. _Chem. Phys. Lett_ 2003, _374_ , 446–452
  * Frediani et al. 2005 Frediani, L.; Rinkevicius, Z.; Ågren, H. Two-photon absorption in solution by means of time-dependent density-functional theory and the polarizable continuum model. _J. Chem. Phys._ 2005, _122_ , 244104
  * Nayyar et al. 2013 Nayyar, I. H.; Masunov, A. E.; Tretiak, S. Comparison of TD-DFT Methods for the Calculation of Two-Photon Absorption Spectra of Oligophenylvinylenes. _J. Phys. Chem. C_ 2013, _117_ , 18170–18189
  * Hättig et al. 1998 Hättig, C.; Christiansen, O.; Jørgensen, P. Coupled cluster response calculations of two-photon transition probability rate constants for helium, neon and argon. _J. Chem. Phys._ 1998, _108_ , 8355–8359
  * Friese et al. 2012 Friese, D. H.; Hättig, C.; Ruud, K. Calculation of two-photon absorption strengths with the approximate coupled cluster singles and doubles model CC2 using the resolution-of-identity approximation. _Phys. Chem. Chem. Phys._ 2012, _14_ , 1175–1184
  * Nanda and Krylov 2015 Nanda, K. D.; Krylov, A. I. Two-photon absorption cross sections within equation-of-motion coupled-cluster formalism using resolution-of-the-identity and Cholesky decomposition representations: Theory, implementation, and benchmarks. _J. Chem. Phys._ 2015, _142_ , 064118
  * De Groot 2001 De Groot, F. High-Resolution X-ray Emission and X-ray Absorption Spectroscopy. _Chem. Rev._ 2001, _101_ , 1779–1808
  * Ament et al. 2011 Ament, L. J. P.; van Veenendaal, M.; Devereaux, T. P.; Hill, J. P.; van den Brink, J. Resonant inelastic x-ray scattering studies of elementary excitations. _Rev. Mod. Phys._ 2011, _83_ , 705–767
  * Fuchs et al. 2008 Fuchs, O.; Zharnikov, M.; Weinhardt, L.; Blum, M.; Weigand, M.; Zubavichus, Y.; Bär, M.; Maier, F.; Denlinger, J. D.; Heske, C.; Grunze, M.; Umbach, E. Isotope and Temperature Effects in Liquid Water Probed by X-Ray Absorption and Resonant X-Ray Emission Spectroscopy. _Phys. Rev. Lett._ 2008, _100_ , 027801
  * Hennies et al. 2010 Hennies, F.; Pietzsch, A.; Berglund, M.; Föhlisch, A.; Schmitt, T.; Strocov, V.; Karlsson, H. O.; Andersson, J.; Rubensson, J.-E. Resonant Inelastic Scattering Spectra of Free Molecules with Vibrational Resolution. _Phys. Rev. Lett._ 2010, _104_ , 193002
  * Kunnus et al. 2016 Kunnus, K.; Zhang, W.; Delcey, M. G.; Pinjari, R. V.; Miedema, P. S.; Schreck, S.; Quevedo, W.; Schröder, H.; Föhlisch, A.; Gaffney, K. J.; Lundberg, M.; Odelius, M.; Wernet, P. Viewing the Valence Electronic Structure of Ferric and Ferrous Hexacyanide in Solution from the Fe and Cyanide Perspectives. _J. Phys. Chem. B_ 2016, _120_ , 7182–7194
  * Kjellsson et al. 2020 Kjellsson, L. et al. Resonant Inelastic X-Ray Scattering Reveals Hidden Local Transitions of the Aqueous OH Radical. _Phys. Rev. Lett._ 2020, _124_ , 236001
  * Rehn et al. 2017 Rehn, D. R.; Dreuw, A.; Norman, P. Resonant Inelastic X-ray Scattering Amplitudes and Cross Sections in the Algebraic Diagrammatic Construction/Intermediate State Representation (ADC/ISR) Approach. _J. Chem. Theory Comput._ 2017, _13_ , 5552–5559
  * Josefsson et al. 2012 Josefsson, I.; Kunnus, K.; Schreck, S.; Föhlisch, A.; De Groot, F.; Wernet, P.; Odelius, M. Ab Initio Calculations of X-ray Spectra: Atomic Multiplet and Molecular Orbital Effects in a Multiconfigurational SCF Approach to the L-Edge Spectra of Transition Metal Complexes. _J. Phys. Chem. Lett._ 2012, _3_ , 3565–3570
  * Nascimento et al. 2021 Nascimento, D. R.; Biasin, E.; Poulter, B. I.; Khalil, M.; Sokaras, D.; Govind, N. Resonant Inelastic X-ray Scattering Calculations of Transition Metal Complexes Within a Simplified Time-Dependent Density Functional Theory Framework. _J. Chem. Theory Comput._ 2021, _17_ , 3031–3038
  * Faber and Coriani 2019 Faber, R.; Coriani, S. Resonant Inelastic X-ray Scattering and Nonesonant X-ray Emission Spectra from Coupled-Cluster (Damped) Response Theory. _J. Chem. Theory Comput._ 2019, _15_ , 520–528
  * Schnack-Petersen et al. 2023 Schnack-Petersen, A. K.; Moitra, T.; Folkestad, S. D.; Coriani, S. New Implementation of an Equation-of-Motion Coupled-Cluster Damped-Response Framework with Illustrative Applications to Resonant Inelastic X-ray Scattering. _J. Phys. Chem. A_ 2023, _127_ , 1775–1793, arXiv:2211.12215 [physics]
  * Skomorowski and Krylov 2021 Skomorowski, W.; Krylov, A. I. Feshbach–Fano approach for calculation of Auger decay rates using equation-of-motion coupled-cluster wave functions. I. Theory and implementation. _J. Chem. Phys._ 2021, _154_ , 084124
  * Skomorowski and Krylov 2021 Skomorowski, W.; Krylov, A. I. Feshbach–Fano approach for calculation of Auger decay rates using equation-of-motion coupled-cluster wave functions. II. Numerical examples and benchmarks. _J. Chem. Phys._ 2021, _154_ , 084125
  * Skeidsvoll et al. 2022 Skeidsvoll, A. S.; Moitra, T.; Balbi, A.; Paul, A. C.; Coriani, S.; Koch, H. Simulating weak-field attosecond processes with a Lanczos reduced basis approach to time-dependent equation-of-motion coupled-cluster theory. _Phys. Rev. A_ 2022, _105_ , 023103
  * Ranga and Dutta 2021 Ranga, S.; Dutta, A. K. A Core–Valence Separated Similarity Transformed EOM-CCSD Method for Core-Excitation Spectra. _J. Chem. Theory Comput._ 2021, _17_ , 7428–7446
  * Henriksson et al. 2005 Henriksson, J.; Norman, P.; Jensen, H. J. A. Two-photon absorption in the relativistic four-component Hartree–Fock approximation. _J. Chem. Phys._ 2005, _122_ , 114106
  * Taube and Bartlett 2005 Taube, A. G.; Bartlett, R. J. Frozen Natural Orbitals: Systematic Basis Set Truncation for Coupled-Cluster Theory. _Collect. Czech. Chem. Commun._ 2005, _70_ , 837–850
  * Taube and Bartlett 2008 Taube, A. G.; Bartlett, R. J. Frozen natural orbital coupled-cluster theory: Forces and application to decomposition of nitroethane. _J. Chem. Phys._ 2008, _128_ , 164101
  * Crawford et al. 2019 Crawford, T. D.; Kumar, A.; Bazanté, A. P.; Di Remigio, R. Reduced‐scaling coupled cluster response theory: Challenges and opportunities. _WIREs Comput Mol Sci_ 2019, _9_
  * Kumar and Crawford 2017 Kumar, A.; Crawford, T. D. Frozen Virtual Natural Orbitals for Coupled-Cluster Linear-Response Theory. _J. Phys. Chem. A_ 2017, _121_ , 708–716
  * Surjuse et al. 2022 Surjuse, K.; Chamoli, S.; Nayak, M. K.; Dutta, A. K. A low-cost four-component relativistic equation of motion coupled cluster method based on frozen natural spinors: Theory, implementation, and benchmark. _J. Chem. Phys._ 2022, _157_ , 204106
  * Yuan et al. 2023 Yuan, X.; Halbert, L.; Pototschnig, J.; Papadopoulos, A.; Coriani, S.; Visscher, L.; Gomes, A. S. P. Formulation and Implementation of Frequency-Dependent Linear Response Properties with Relativistic Coupled Cluster Theory for GPU-accelerated Computer Architectures. 2023; http://arxiv.org/abs/2307.14296, arXiv:2307.14296 [physics]
  * McClain 1971 McClain, W. M. Excited State Symmetry Assignment Through Polarized Two-Photon Absorption Studies of Fluids. _J. Chem. Phys._ 1971, _55_ , 2789–2796
  * Saue et al. 2020 Saue, T.; Bast, R.; Gomes, A. S. P.; Jensen, H. J. A.; Visscher, L.; Aucar, I. A.; Di Remigio, R.; Dyall, K. G.; Eliav, E.; Fasshauer, E., et al. The DIRAC code for relativistic molecular calculations. _J. Chem. Phys_ 2020, _152_ , 204104
  * 70 DIRAC, a relativistic ab initio electronic structure program, Release DIRAC23 (2023), written by R. Bast, A. S. P. Gomes, T. Saue and L. Visscher and H. J. Aa. Jensen, with contributions from I. A. Aucar, V. Bakken, C. Chibueze, J. Creutzberg, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, L. Halbert, E. D. Hedegård, T. Helgaker, B. Helmich–Paris, J. Henriksson, M. van Horn, M. Iliaš, Ch. R. Jacob, S. Knecht, S. Komorovský, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, N. H. List, H. S. Nataraj, M. K. Nayak, P. Norman, A. Nyvang, G. Olejniczak, J. Olsen, J. M. H. Olsen, A. Papadopoulos, Y. C. Park, J. K. Pedersen, M. Pernpointner, J. V. Pototschnig, R. di Remigio, M. Repisky, K. Ruud, P. Sałek, B. Schimmelpfennig, B. Senjean, A. Shee, J. Sikkema, A. Sunaga, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, M. L. Vidal, S. Villaume, O. Visser, T. Winther, S. Yamamoto and X. Yuan (available at https://doi.org/10.5281/zenodo.7670749, see also https://www.diracprogram.org)
  * Dyall et al. 2022 Dyall, K. G.; Tecmer, P.; Sunaga, A. Diffuse Basis Functions for Relativistic s and d Block Gaussian Basis Sets. _J. Chem. Theory. Comput_ 2022,
  * Dyall 2006 Dyall, K. G. Relativistic quadruple-zeta and revised triple-zeta and double-zeta basis sets for the 4p, 5p, and 6p elements. _Theor. Chem. Acc_ 2006, _115_ , 441–447
  * Kendall et al. 1992 Kendall, R. A.; Dunning Jr, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. _J. Chem. Phys_ 1992, _96_ , 6796–6806
  * Woon and Dunning Jr 1993 Woon, D. E.; Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. _J. Chem. Phys_ 1993, _98_ , 1358–1371
  * Wilson et al. 1999 Wilson, A. K.; Woon, D. E.; Peterson, K. A.; Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. IX. The atoms gallium through krypton. _J. Chem. Phys_ 1999, _110_ , 7667–7676
  * Iliaš and Saue 2007 Iliaš, M.; Saue, T. An infinite-order two-component relativistic Hamiltonian by a simple one-step transformation. _J. Chem. Phys_ 2007, _126_ , 064102
  * Lévy-Leblond 1967 Lévy-Leblond, J.-M. Nonrelativistic particles and wave equations. _Commun. Math. Phys_ 1967, _6_ , 286–311
  * Visscher and Saue 2000 Visscher, L.; Saue, T. Approximate relativistic electronic structure methods based on the quaternion modified Dirac equation. _J. Chem. Phys_ 2000, _113_ , 3996–4002
  * Becke 1993 Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. _J. Chem. Phys_ 1993, _98_ , 1372–1377
  * Visscher and Dyall 1997 Visscher, L.; Dyall, K. G. Dirac–Fock atomic electronic structure calculations using different nuclear charge distributions. _Atom. Data. Nucl. Data_ 1997, _67_ , 207–224
  * Beebe and Linderberg 1977 Beebe, N. H.; Linderberg, J. Simplifications in the generation and transformation of two-electron integrals in molecular calculations. _Int. J. Quantum. Chem_ 1977, _12_ , 683–705
  * Koch et al. 2003 Koch, H.; Sánchez de Merás, A.; Pedersen, T. B. Reduced scaling in electronic structure calculations using Cholesky decompositions. _J. Chem. Phys_ 2003, _118_ , 9481–9484
  * Aquilante et al. 2007 Aquilante, F.; Pedersen, T. B.; Lindh, R. Low-cost evaluation of the exchange Fock matrix from Cholesky and density fitting representations of the electron repulsion integrals. _J. Chem. Phys_ 2007, _126_ , 194106
  * Huber 1979 Huber, K.-P. _Molecular spectra and molecular structure_ 1979,
  * Gomes and Visscher 2004 Gomes, A. S. P.; Visscher, L. The influence of core correlation on the spectroscopic constants of HAt. _Chem. Phys. Lett._ 2004, _399_ , 1
  * Aidas et al. 2014 Aidas, K. et al. The Dalton quantum chemistry program system: The Dalton program. _WIREs Comput Mol Sci_ 2014, _4_ , 269–284
  * Rizzo et al. 2002 Rizzo, A.; Coriani, S.; Fernández, B.; Christiansen, O. A coupled cluster response study of the electric dipole polarizability, first and second hyperpolarizabilities of HCl. _Phys. Chem. Chem. Phys._ 2002, _4_ , 2884–2890
  * Bishop and Norman 1999 Bishop, D. M.; Norman, P. Effects of vibration on the polarizability and the first and second hyperpolarizabilities of HF, HCl, and HBr. _J. Chem. Phys._ 1999, _111_ , 3042–3050
  * Dudley and Ward 1985 Dudley, J.; Ward, J. Measurements of second- and third-order nonlinear polarizabilities for HF and HCI. _J. Chem. Phys._ 1985, _82_ , 4673–3677
  * Coriani et al. 1997 Coriani, S.; Hättig, C.; Jørgensen, P.; Halkier, A.; Rizzo, A. Coupled cluster calculations of Verdet constants. _Chem. Phys. Lett_ 1997, _281_ , 445–451
  * Pecul et al. 2006 Pecul, M.; Pawłowski, F.; Jørgensen, P.; Köhn, A.; Hättig, C. High-order correlation effects on dynamic hyperpolarizabilities and their geometric derivatives: A comparison with density functional results. _J. Chem. Phys._ 2006, _124_ , 114101
  * Visscher 1997 Visscher, L. Approximate molecular relativistic Dirac-Coulomb calculations using a simple Coulombic correction. _Theor. Chem. Acc._ 1997, _98_ , 68–70
  * Dyall 1994 Dyall, K. G. An exact separation of the spin-free and spin-dependent terms of the Dirac–Coulomb–Breit Hamiltonian. _J. Chem. Phys._ 1994, _100_ , 2118–2127
  * Ingersoll and Liebenberg 1956 Ingersoll, L. R.; Liebenberg, D. H. Faraday Effect in Gases and Vapors. II$\ast$. _J. Opt. Soc. Am._ 1956, _46_ , 538–542
  * Zuhrianda et al. 2012 Zuhrianda, Z.; Safronova, M. S.; Kozlov, M. G. Anomalously small blackbody radiation shift in the Tl + frequency standard. _Phys. Rev. A_ 2012, _85_ , 022513
  * Safronova et al. 2011 Safronova, M. S.; Kozlov, M. G.; Clark, C. W. Precision Calculation of Blackbody Radiation Shifts for Optical Frequency Metrology. _Phys. Rev. Lett._ 2011, _107_ , 143006
  * Yuan et al. 2022 Yuan, X.; Visscher, L.; Gomes, A. S. P. Assessing MP2 frozen natural orbitals in relativistic correlated electronic structure calculations. _J. Chem. Phys._ 2022, _156_ , 224108
  * Höfener and Klopper 2017 Höfener, S.; Klopper, W. Natural transition orbitals for the calculation of correlation and excitation energies. _Chem. Phys. Lett_ 2017, _679_ , 52–59