Second-Order Multi-Reference Algebraic Diagrammatic Construction Theory for Photoelectron Spectra of Strongly Correlated Systems

TL;DR

This paper introduces MR-ADC(2), a second-order multi-reference algebraic diagrammatic construction method that efficiently simulates photoelectron spectra of strongly correlated systems, capturing ionization energies and intensities across all orbitals.

Contribution

The paper develops a new second-order MR-ADC(2) method that avoids large density matrices and offers lower computational scaling for photoelectron spectrum simulations.

Findings

Successfully applied to small molecules, carbon dimer, and hydrogen chains.

Provides comprehensive ionization information for all orbitals.

Achieves lower computational cost compared to traditional methods.

Abstract

We present a second-order formulation of multi-reference algebraic diagrammatic construction theory [Sokolov, A. Yu. J. Chem. Phys. 2018, 149, 204113] for simulating photoelectron spectra of strongly correlated systems (MR-ADC(2)). The MR-ADC(2) method uses second-order multi-reference perturbation theory (MRPT2) to efficiently obtain ionization energies and intensities for many photoelectron transitions in a single computation. In contrast to conventional MRPT2 methods, MR-ADC(2) provides information about ionization of electrons in all orbitals (i.e., core and active) and allows to compute transition intensities in straightforward and efficient way. Although equations of MR-ADC(2) depend on four-particle reduced density matrices, we demonstrate that computation of these large matrices can be completely avoided without introducing any approximations. The resulting MR-ADC(2) implementation has a lower computational scaling compared to conventional MRPT2 methods. We present results of MR-ADC(2) for photoelectron spectra of small molecules, carbon dimer, and equally-spaced hydrogen chains (H10 and H30) and outline directions for future developments.

Figures

b](#S2.F3.sf2)). Conveniently, the $\mathbf{S}^{-1/2}$ matrix can be constructed together with the $(\mathbf{S^{[i]}})^{-1/2}$ matrices used for solution of the amplitude equations ([Section 3.2](#S3.SS2)). As an example, we consider non-zero elements of $\mathbf{S}$ for $a_{ij}^{x}\ket{\Psi_{0}}$ that have the form $S_{ijx,ijy}=\braket{\Psi_{0}}{a_{y}^{ij}a_{ij}^{x}}{\Psi_{0}}=\braket{\Psi_{0}}{a_{y}a^{\dagger}_{x}}{\Psi_{0}}$. These elements are equal to the $\mathbf{S^{[+1]}}$ matrix elements $S_{ijay,ijax}^{[+1]}=\braket{\Psi_{0}}{a_{ay}^{ij}a_{ij}^{ax}}{\Psi_{0}}=\braket{\Psi_{0}}{a_{y}a^

b](#S5.F4.sf2)). MR-ADC(2) shows performance similar to that for equilibrium geometries ([Table 2](#S5.T2)), with $\Delta_{\mathrm{MAE}}$ (0.50 eV) and $\Delta_{\mathrm{STD}}$ (0.36 eV) smaller than the corresponding errors for the single-reference methods. The best agreement with FCI is again shown by QD-NEVPT2 with $\Delta_{\mathrm{MAE}}$ = 0.18 eV and $\Delta_{\mathrm{STD}}$ = 0.23 eV. As for the equilibrium geometries, the MR-ADC(2) ionization energies for stretched geometries are systematically overestimated relative to FCI, reproducing energy spacings between ionized states within 0.1 eV

a](#S5.F8.sf1) shows MR-ADC(2) DOS computed using the STO-6G and cc-pVDZ basis sets. For both geometries, DOS computed using MR-ADC(2) with the STO-6G basis closely resembles LDOS of the same system from the DMRG study of Ronca *et al*.[147](#bib.bib147) [Figure 8b](#S5.F8.sf2) plots contributions to MR-ADC(2)/STO-6G DOS from core and active orbitals separately. For the compressed chain ($r$ = 1.4 $a_{0}$), contributions from active orbitals dominate the low-energy part of the spectrum, whereas, for equilibrium geometry ($r$ = 1.8 $a_{0}$), core and active orbitals have similar contributions t

Tables

Table 1. Table 1: Size-consistency errors of the MR-ADC(2) ionization energies ( Δ ​ Ω Δ Ω \Delta\Omega , eV) and spectroscopic factors ( Δ ​ P Δ 𝑃 \Delta P ) for the \ce (H2O)2 and \ce (HF)2 systems composed of two identical monomers separated by 10000 Å (aug-cc-pVDZ basis set). For \ce H2O, r e subscript 𝑟 𝑒 r_{e} = r 𝑟 r ( \ce O-H) = 1.0 Å and ∠ ∠ \angle (H–O–H) = 104.5 ∘ . For \ce HF, r e subscript 𝑟 𝑒 r_{e} = 0.917 Å . The (4e, 4o) and (6e, 5o) active spaces were used for the \ce H2O and \ce HF monomer CASSCF reference wavefunctions. For dimers, (8e, 8o) and (12e, 10o) active spaces were used, respectively. The number of CASCI ionized states was set to 10 and 20, for monomers and dimers, respectively.

System	State	$\Delta \Omega$	$\Delta P$
\ce (H2O)2 ($r_e$)	$1b_1$	$-$2.4 $\times$ ${10}^{-5}$	6.0 $\times$ ${10}^{-8}$
	$3a_1$	9.4 $\times$ ${10}^{-6}$	$-$2.7 $\times$ ${10}^{-7}$
	$1b_2$	2.8 $\times$ ${10}^{-6}$	2.5 $\times$ ${10}^{-7}$
\ce (H2O)2 ($2r_e$)	$1b_1$	$-$1.3 $\times$ ${10}^{-5}$	3.4 $\times$ ${10}^{-7}$
	$3a_1$	1.5 $\times$ ${10}^{-5}$	2.3 $\times$ ${10}^{-6}$
	$1b_2$	1.3 $\times$ ${10}^{-5}$	1.8 $\times$ ${10}^{-6}$
\ce (HF)2 ($r_e$)	$1\pi$	1.2 $\times$ ${10}^{-4}$	$-$1.0 $\times$ ${10}^{-6}$
	$3\sigma$	5.4 $\times$ ${10}^{-5}$	2.4 $\times$ ${10}^{-6}$
\ce (HF)2 ($2r_e$)	$1\pi$	1.0 $\times$ ${10}^{-4}$	$-$4.5 $\times$ ${10}^{-6}$
	$3\sigma$	5.9 $\times$ ${10}^{-5}$	$-$9.1 $\times$ ${10}^{-7}$

Table 2. Table 2: Computed vertical ionization energies ( Ω Ω \Omega , eV) and spectroscopic factors ( P 𝑃 P ) of molecules with equilibrium geometries. See Section   4 for active spaces used in the reference CASSCF computations, structural parameters, and basis sets. Also shown are mean absolute errors ( Δ MAE subscript Δ MAE \Delta_{\mathrm{MAE}} ) and standard deviations ( Δ STD subscript Δ STD \Delta_{\mathrm{STD}} ) of the results, relative to FCI.

System	State	SR-ADC(2)		SR-ADC(3)		SR-ADC(3+)a		MR-ADC(2)		EOM-CCSD	QD-NEVPT2	FCI
		$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$\Omega$	$\Omega$
HF	$1\pi$	14.41	0.89	16.79	0.93	16.41	0.93	16.35	0.93	15.85	16.00	16.07
	$3\sigma$	18.69	0.90	20.65	0.94	20.30	0.94	20.38	0.94	19.88	20.04	20.06
F2	$1{\pi }_g$	13.90	0.87	16.03	0.89	15.87	0.90	16.55	0.88	15.40	15.38	15.64
	$1{\pi }_u$	17.06	0.84	19.25	0.75	19.11	0.81	19.86	0.80	18.77	18.58	18.83
	$3{\sigma }_g$	20.25	0.89	21.26	0.89	21.01	0.88	22.08	0.87	21.16	20.88	21.15
CO	$5\sigma$	13.78	0.91	13.57	0.90	13.80	0.89	14.07	0.92	13.99	13.53	13.74
	$1\pi$	16.24	0.89	17.16	0.90	16.88	0.90	17.38	0.90	16.93	16.75	16.90
	$4\sigma$	18.28	0.85	20.46	0.76	20.10	0.79	20.15	0.85	19.67	19.48	19.56
N2	$3{\sigma }_g$	14.79	0.88	15.42	0.91	15.60	0.91	15.76	0.91	15.43	15.21	15.30
	$1{\pi }_u$	16.98	0.91	16.60	0.92	16.77	0.92	17.33	0.92	17.11	16.75	16.83
	$2{\sigma }_u$	17.96	0.85	18.79	0.82	18.93	0.82	19.00	0.83	18.71	18.44	18.50
H2O	$1b_1$	11.23	0.89	12.99	0.92	12.78	0.92	12.74	0.93	12.38	12.55	12.53
	$3a_1$	13.53	0.89	15.28	0.92	15.08	0.93	15.07	0.93	14.66	14.85	14.81
	$1b_2$	17.95	0.90	19.34	0.93	19.16	0.93	19.28	0.94	18.89	19.05	18.98
CS	$7\sigma$	10.99	0.86	10.99	0.85	11.33	0.85	11.59	0.85	11.36	10.95	11.13
	$2\pi$	12.84	0.91	12.67	0.90	12.66	0.90	13.43	0.91	12.94	12.74	12.83
	$6\sigma$	16.88	0.85	15.53	0.18	15.51	0.19	16.83	0.40	17.02	15.83	15.88
H2CO	$2b_2$	9.46	0.87	11.11	0.91	10.87	0.91	11.23	0.92	10.62	10.28	10.72
	$1b_1$	13.73	0.88	14.54	0.88	14.30	0.88	15.14	0.90	14.47	14.07	14.48
	$5a_1$	14.62	0.86	16.61	0.90	16.20	0.90	16.70	0.90	15.95	15.64	16.01
	$1b_2$	16.67	0.88	17.04	0.69	17.32	0.65	17.76	0.88	17.21	16.50	16.86
C2H4	$1b_{1u}$	10.14	0.91	10.47	0.91	10.46	0.91	11.01	0.90	10.58	10.41	10.58
	$1b_{1g}$	12.79	0.91	13.22	0.91	13.19	0.91	13.75	0.92	13.22	13.05	13.21
	$3a_g$	13.78	0.89	14.34	0.91	14.36	0.91	14.74	0.89	14.31	14.12	14.25
	$1b_{2u}$	16.13	0.87	16.50	0.74	16.49	0.79	17.10	0.84	16.61	16.35	16.45
${\Delta }_{\mathrm{MAE}}$		0.83		0.30		0.21		0.56		0.17	0.17
${\Delta }_{\mathrm{STD}}$		0.68		0.32		0.22		0.23		0.28	0.14

Table 3. Table 3: Computed vertical ionization energies ( Ω Ω \Omega , eV) and spectroscopic factors ( P 𝑃 P ) of molecules with stretched geometries. See Section   4 for active spaces used in the reference CASSCF computations, structural parameters, and basis sets. Also shown are mean absolute errors ( Δ MAE subscript Δ MAE \Delta_{\mathrm{MAE}} ) and standard deviations ( Δ STD subscript Δ STD \Delta_{\mathrm{STD}} ) of the results, relative to FCI.

System	State	SR-ADC(2)		SR-ADC(3)		MR-ADC(2)		EOM-CCSD	QD-NEVPT2	FCI
		$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$\Omega$
HF	$1\pi$	9.84	0.77	16.15	0.84	13.86	0.60	13.67	13.61	13.65
	$3\sigma$	13.30	0.84	14.68	0.76	14.98	0.73	14.76	14.83	14.84
F2	$1{\pi }_g$	10.63	0.64	17.55	0.88	18.12	0.74	16.86	16.81	17.13
	$1{\pi }_u$	10.66	0.64	17.69	0.89	18.16	0.82	16.95	16.87	17.19
N2	$3{\sigma }_g$	15.70	0.63	$-$2.60	1.69	14.00	0.69	14.36	13.06	13.38
	$1{\pi }_u$	17.50	0.55	$-$5.24	2.16	14.17	0.51	14.77	13.21	13.49
H2O	$1b_1$	6.53	0.71	12.24	0.66	11.31	0.64	10.65	10.99	11.07
	$3a_1$	10.49	0.75	12.78	0.67	13.22	0.67	12.69	12.99	13.02
	$1b_2$	11.18	0.75	13.01	0.72	13.78	0.71	13.26	13.53	13.56
H2CO	$2b_2$	10.65	0.85	8.31	0.21	11.51	0.39	9.85	10.24	10.37
	$1b_1$	10.69	0.86	8.35	0.22	11.21	0.48	9.66	10.38	10.55
	$5a_1$	10.60	0.91	10.97	0.88	13.16	0.57	10.97	12.32	13.16
C2H4	$1b_{1u}$	9.37	0.76	6.87	0.83	9.69	0.53	9.41	9.26	9.25
	$3a_g$	11.38	0.79	8.74	0.91	11.36	0.73	11.17	10.94	10.93
${\Delta }_{\mathrm{MAE}}$		2.70		3.66		0.50		0.56	0.18
${\Delta }_{\mathrm{STD}}$		3.10		6.28		0.36		0.81	0.23

Table 4. Table 4: Vertical ionization energies ( Ω Ω \Omega , eV) and spectroscopic factors ( P 𝑃 P ) of carbon dimer with r 𝑟 r \ce (C-C) = 1.2425 Å computed using the aug-cc-pVDZ basis set. For MR-ADC(2) and QD-NEVPT2, the CASSCF reference wavefunction was computed using the (8e, 12o) active space.

Configuration	State	SR-ADC(3)		MR-ADC(2)		QD-NEVPT2	FCI
		$\Omega$	$P$	$\Omega$	$P$	$\Omega$	$\Omega$
${\left(2{\sigma }_u\right)}^2{\left(1{\pi }_u\right)}^3{\left(3{\sigma }_g\right)}^0$	$1^2{\Pi }_u$	11.69	0.9215	12.50	0.8986	12.28	12.34
${\left(2{\sigma }_u\right)}^2{\left(1{\pi }_u\right)}^2{\left(3{\sigma }_g\right)}^1$	$1^2{\Delta }_g$	11.17	0.0002	14.31	0.0002	13.92	13.94
${\left(2{\sigma }_u\right)}^2{\left(1{\pi }_u\right)}^2{\left(3{\sigma }_g\right)}^1$	$1^2{\Sigma }_g^{-}$	a	a	14.55	0.0000	14.12	14.15
${\left(2{\sigma }_u\right)}^2{\left(1{\pi }_u\right)}^2{\left(3{\sigma }_g\right)}^1$	$1^2{\Sigma }_g^{+}$	11.43	0.0004	14.60	0.0047	14.26	14.29
${\left(2{\sigma }_u\right)}^1{\left(1{\pi }_u\right)}^4{\left(3{\sigma }_g\right)}^0$	$1^2{\Sigma }_u^{+}$	13.95	0.8738	15.33	0.7190	15.07	15.09
${\left(2{\sigma }_u\right)}^2{\left(1{\pi }_u\right)}^1{\left(3{\sigma }_g\right)}^2$	$2^2{\Pi }_u$	a	a	14.77	0.0183	15.35	15.43

Equations

$$G_{\mu \nu }(\omega )$$

$$=⟨\Psi |q_{\mu }(\omega -H+E{)}^{-1}{q}_{\nu }^{†}|\Psi ⟩$$

$$\pm ⟨\Psi |q_{\nu }^{†}(\omega +H-E{)}^{-1}{q}_{\mu }|\Psi ⟩$$

$$|\Psi ⟩$$

$$T_k$$

$$H^{(0)}$$

$$H_{act}$$

$$C$$

$$f_p^q$$

$$\mathbf{G}(\omega )$$

$${\mathbf{G}}_{\pm }(\omega )$$

$${\mathbf{M}}_{\pm }{\mathbf{Y}}_{\pm }={\mathbf{S}}_{\pm }{\mathbf{Y}}_{\pm }{\mathbf{\Omega }}_{\pm }$$

$${\mathbf{X}}_{\pm }={\mathbf{T}}_{\pm }{\mathbf{S}}_{\pm }^{-1/2}{\mathbf{Y}}_{\pm }$$

$${\mathbf{G}}_{\pm }(\omega )$$

$$T(\omega )$$

$$M_{\mu \nu }^{(n)}$$

$$T_{\mu \nu }^{(n)}$$

$$S_{\mu \nu }^{(n)}$$

$${\tilde{H}}^{(0)}$$

$${\tilde{H}}^{(1)}$$

$${\tilde{H}}^{(2)}$$

$${\tilde{q}}_{\mu }^{(0)}$$

$${\tilde{q}}_{\mu }^{(1)}$$

$${\tilde{q}}_{\mu }^{(2)}$$

$$\{\text{A}|\text{B}\}$$

$$\{\text{A}|\tilde{\mathcal{H}}|\text{B}\}$$

$$\mathbf{M}$$

$$+\{{\mathbf{h}}^{(0)†}|{\tilde{\mathcal{H}}}^{(1)}|{\mathbf{h}}^{(0)†}\}+\{{\mathbf{h}}^{(0)†}|{\tilde{\mathcal{H}}}^{(0)}|{\mathbf{h}}^{(1)†}\}$$

$$+\{{\mathbf{h}}^{(1)†}|{\tilde{\mathcal{H}}}^{(1)}|{\mathbf{h}}^{(0)†}\}+\{{\mathbf{h}}^{(1)†}|{\tilde{\mathcal{H}}}^{(0)}|{\mathbf{h}}^{(1)†}\}$$

$$+\{{\mathbf{h}}^{(0)†}|{\tilde{\mathcal{H}}}^{(1)}|{\mathbf{h}}^{(1)†}\}+\{{\mathbf{h}}^{(0)†}|{\tilde{\mathcal{H}}}^{(2)}|{\mathbf{h}}^{(0)†}\}$$

$$\mathbf{T}$$

$$+\{{\tilde{\mathbf{q}}}^{(1)}|{\mathbf{h}}^{(1)†}\}+\{{\tilde{\mathbf{q}}}^{(2)}|{\mathbf{h}}^{(0)†}\}$$

$$\mathbf{S}$$

$$+\{{\mathbf{h}}^{(0)†}|{\mathbf{h}}^{(1)†}\}+\{{\mathbf{h}}^{(1)†}|{\mathbf{h}}^{(1)†}\}$$

$$T^{(k)}={\mathbf{t}}^{(\mathbf{k})}\mathbf{\tau }=\sum _{\mu }{t}_{\mu }^{(k)}{\tau }_{\mu }$$

$${\mathbf{t}}^{(1)}=$$

$$t_{ij}^{xy(1)};t_{xy}^{ab(1)};t_{ix}^{ay(1)};t_{ix}^{yz(1)};t_{xy}^{az(1)}\}$$

$$\mathbf{\tau }=$$

$$a_{ij}^{xy};a_{xy}^{ab};a_{ix}^{ay};a_{ix}^{yz};a_{xy}^{az}\}$$

$$⟨{\Psi }_0⟩{\tau }_{\mu }^{†}{\tilde{H}}^{(1)}{\Psi }_0=0$$

$${\mathbf{H}}^{(0)}{\mathbf{t}}^{(1)}=-{\mathbf{V}}^{(1)}$$

$$H_{\mu \nu }^{(0)}$$

Full text

1 Introduction

Recently, there has been a significant progress in increasing tractability of strong electron correlation problem. New methods enable computations of systems with a large number of strongly correlated electrons in the ground or excited electronic states.1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 These approaches usually start by computing a multi-configurational wavefunction that describes strong correlation in a subset of frontier (active) molecular orbitals with near-degeneracies.15, 16, 17 The remaining (dynamic) correlation effects outside of the active orbitals are usually captured by multi-reference perturbation theory (MRPT),18, 19, 20, 21, 22, 23, 24, 25, 26, 27 configuration interaction,28, 29, 30, 31, 32 or coupled cluster (CC) methods.33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50 In particular, low-order MRPT methods have been very successful at computing accurate energies of large strongly correlated systems, due to their relatively low computation cost and ability to treat large active spaces with up to $\sim$ 30 orbitals.51, 52, 53, 54, 55, 56, 57, 58

Despite significant advances, application of conventional MRPT methods to a wider range of problems, such as simulating excited-state or spectroscopic properties, is hindered by a number of limitations. For example, computation of transition intensities in MRPT is not straightforward due to complexity of the underlying response equations.59 Another limitation is that MRPT methods do not describe electronic transitions involving orbitals outside active space that are important for simulating broadband spectra or core-level excitations in X-ray spectroscopies. Furthermore, for computations involving many electronic states of the same symmetry, MRPT methods rely on using state-averaged reference wavefunctions, which introduce dependence of their results on the number of states and weights used in state-averaging. This motivates the development of new efficient multi-reference theories that are not bound by these limitations.

We have recently proposed a multi-reference formulation of algebraic diagrammatic construction theory (MR-ADC) for simulating spectroscopic properties of strongly correlated systems.60 MR-ADC is a generalization of the conventional (single-reference) ADC theory proposed by Schirmer in 1982.61 Rather than computing energies and wavefunctions of individual electronic states, in MR-ADC excitation energies and transition intensities are directly obtained from poles and residues of a retarded propagator approximated using multi-reference perturbation theory. In contrast to conventional MRPT, MR-ADC describes electronic transitions involving all orbitals (i.e., core, active, and external), enables simulations of various spectroscopic processes (e.g., ionization or two-photon excitation), and provides direct access to spectral properties. In this regard, MR-ADC is related to multi-reference propagator theories,62, 63, 64, 65, 66, 67, 68, 69, 70 but has an advantage of a Hermitian eigenvalue problem and including dynamic correlation effects beyond single excitations. For electronic excitations, MR-ADC can also be considered as a low-cost alternative to multi-reference equation-of-motion (MR-EOM) theories, such as MR-EOM-CC,41, 42, 43 and internally-contracted linear-response theories, such as ic-MRCC.71

In this work, we present a second-order formulation of MR-ADC (MR-ADC(2)) for photoelectron spectra of multi-reference systems. We begin by describing the derivation of MR-ADC(2) (Section 2) and discuss details of its implementation (Section 3), demonstrating that it has a lower computational scaling with the number of active orbitals compared to conventional MRPT methods. Next, we describe computational details (Section 4) and test the performance of MR-ADC(2) for computing photoelectron energies and transition intensities of small molecules, carbon dimer, as well as equally-spaced hydrogen chains \ceH10 and \ceH30 (Section 5). Finally, we present our conclusions (Section 6) and outline future developments.

2 Theory

2.1 Multi-Reference Algebraic Diagrammatic Construction Theory (MR-ADC)

We begin with a brief overview of MR-ADC. In Ref. 60, we have described the derivation of MR-ADC using the formalism of effective Liouvillean theory.72 Here, we only summarize the main results. Our starting point is a general expression for the retarded propagator73, 74 that describes response of a many-electron system to an external perturbation with frequency $\omega$:

[TABLE]

Here, $G_{\mu\nu}^{+}(\omega)$ and $G_{\mu\nu}^{-}(\omega)$ are the forward and backward components of the propagator, $\ket{\Psi}$ and $E$ are the eigenfunction and eigenvalue of the electronic Hamiltonian $H$, and the frequency $\omega\equiv\omega^{\prime}+i\eta$ is written in terms of its real component ($\omega^{\prime}$) and an infinitesimal imaginary number ($i\eta$). Depending on the form of operators $q^{\dagger}_{\nu}$, the propagator $G_{\mu\nu}(\omega)$ can describe various spectroscopic processes. Choosing $q^{\dagger}_{\nu}=a^{\dagger}_{p}a_{q}-\braket{\Psi}{a^{\dagger}_{p}a_{q}}{\Psi}$, where $a^{\dagger}_{p}$ and $a_{p}$ are the usual creation and annihilation operators, corresponds to polarization propagator that provides information about electronic excitations in optical (e.g., UV/Vis) spectroscopy. Alternatively, a propagator with $q^{\dagger}_{\nu}=a^{\dagger}_{p}$ describes electron attachment and ionization processes. The number of creation and annihilation operators in $q^{\dagger}_{\nu}$ (odd or even) determines the sign ($+$ or $-$) of the second term in Section 2.1.

Evaluation of the exact propagator is very expensive computationally. For this reason, many approximate methods75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 61, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 62, 63, 64, 65, 66, 67, 68, 69 have been developed to compute $G_{\mu\nu}(\omega)$ for realistic systems. A common assumption in most of these approaches is that the eigenfunction $\ket{\Psi}$ can be well approximated by a single Slater determinant. Although this assumption significantly simplifies the underlying equations, such single-reference methods do not provide reliable results when strong correlation is important and the wavefunction $\ket{\Psi}$ becomes multi-configurational.

To efficiently and accurately compute $G_{\mu\nu}(\omega)$ for strongly correlated systems, in MR-ADC we consider an expansion of Section 2.1 using multi-reference perturbation theory, where the zeroth-order (reference) wavefunction $\ket{\Psi_{0}}$ is obtained by solving the complete active space configuration interaction (CASCI) or self-consistent field (CASSCF) variational problem in a set of active molecular orbitals (Figure 1). The eigenfunction $\ket{\Psi}$ is related to $\ket{\Psi_{0}}$ via a unitary transformation44, 45, 46, 47, 48, 27, 49

[TABLE]

where $T$ generates all internally-contracted excitations between core, active, and external orbitals (see Figure 1 for orbital index notation). Defining the zeroth-order Hamiltonian to be the Dyall Hamiltonian107, 24, 25, 26

[TABLE]

expressed in the basis of diagonal core and external generalized Fock operators ($f_{{i}}^{{j}}=\varepsilon_{i}\delta_{{i}}^{{j}}$, $f_{{a}}^{{b}}=\varepsilon_{a}\delta_{{a}}^{{b}}$), we expand the propagator in Section 2.1 in perturbative series with respect to the perturbation $V=H-H^{(0)}$:

[TABLE]

Truncating Eq. 8 at the $n$th order in perturbation theory corresponds to the propagator of the MR-ADC(n) approximation.

An important property of MR-ADC (along with that of its single-reference variant)72 is that the forward and backward components of the propagator in Section 2.1 are decoupled and, thus, perturbative expansion (8) can be performed for $G_{\mu\nu}^{+}(\omega)$ and $G_{\mu\nu}^{-}(\omega)$ separately. The MR-ADC(n) $G_{\mu\nu}^{+}(\omega)$ and $G_{\mu\nu}^{-}(\omega)$ contributions are expressed in the matrix form

[TABLE]

where $\mathbf{M}_{\pm}$, $\mathbf{T}_{\pm}$, and $\mathbf{S}_{\pm}$ are the effective Liouvillean, transition moment, and overlap matrices, respectively, each evaluated up to $n$th order in perturbation theory. The $\mathbf{M}_{\pm}$ matrix contains information about transition energies, which are obtained by solving the Hermitian generalized eigenvalue problem

[TABLE]

where $\boldsymbol{\Omega}_{\pm}$ is a diagonal matrix of eigenvalues. The eigenvectors $\mathbf{Y}_{\pm}$ are used to compute spectroscopic amplitudes

[TABLE]

which are related to transition intensities. Combining the eigenvalues $\boldsymbol{\Omega}_{\pm}$ and spectroscopic amplitudes $\mathbf{X}_{\pm}$, we obtain expressions for the MR-ADC(n) propagator and spectral function

[TABLE]

2.2 Second-Order MR-ADC for Ionization Energies and Spectra

2.2.1 Overview

In this work, we consider the MR-ADC(2) approximation for photoelectron spectra, which incorporates all contributions to $\mathbf{G}(\omega)$ up to the second order in perturbation theory. A propagator of choice for the description of electron ionization processes is the backward component of the one-particle Green’s function $\mathbf{G}_{-}(\omega)$, which can be defined by specifying $q^{\dagger}_{\nu}=a^{\dagger}_{p}$ in the second term of Section 2.1. To simplify our notation, we will drop the subscript $-$ everywhere in the equations. Thus, matrices $\mathbf{M}$, $\mathbf{T}$, and $\mathbf{S}$ will refer to the components of $\mathbf{G}_{-}(\omega)$ in Eq. 9. Following the effective Liouvillean approach,72, 60 we express the $n$th-order MR-ADC matrices as:

[TABLE]

where $[\ldots]$ and $[\ldots]_{+}$ denote commutator and anticommutator, respectively. In Eqs. 14, 15 and 16, $\tilde{H}^{(k)}$ and $\tilde{q}_{\mu}^{(k)}$ are the $k$th-order contributions to the effective Hamiltonian $\tilde{H}=e^{-A}He^{A}$ and observable $\tilde{q}_{\mu}=e^{-A}q_{\mu}e^{A}$ operators. These contributions can be obtained by expanding $\tilde{H}$ and $\tilde{q}_{\mu}$ using the Baker–Campbell–Hausdorff (BCH) formula and collecting terms at the $k$th order. The low-order components of these operators have the form

[TABLE]

where $A^{(k)}=T^{(k)}-T^{(k){\dagger}}$ as shown in Eq. 2. The operators $h_{\mu}^{(k){\dagger}}$ compose the $k$th-order ionization operator manifold that is used to construct a set of internally-contracted (ionized) basis states $\ket{\Psi_{\mu}^{(k)}}=h_{\mu}^{(k){\dagger}}\ket{\Psi_{0}}$ necessary for representing the eigenstates in Eq. 10.

Introducing shorthand notations72 for the matrix elements of arbitrary operator sets $\mathbf{A}=\{A_{\mu}\}$ and $\mathbf{B}=\{B_{\mu}\}$

[TABLE]

we express contributions to the MR-ADC(2) matrices in the following form

[TABLE]

Computing matrix elements in Sections 2.2.1, 2.2.1 and 2.2.1 requires solving for amplitudes of the excitation operators ($T^{(1)}$ and $T^{(2)}$) and determining the ionization operator manifolds ($h_{\mu}^{(k){\dagger}}$, $k=0,1$).

2.2.2 Amplitudes of the Excitation Operators

To solve for amplitudes of the $T^{(k)}$ $(k=1,2)$ operators, we express these operators in a general form

[TABLE]

where $t_{\mu}^{(k)}$ are the $k$th-order coefficients and $\tau_{\mu}$ are the corresponding excitation operators (Eq. 3). The first-order operator $T^{(1)}$ includes up to two-body terms ($T^{(1)}=T^{(1)}_{1}+T^{(1)}_{2}$) parametrized using three classes of single excitation and eight classes of double excitation amplitudes

[TABLE]

Defining $a_{q}^{p}\equiv a^{\dagger}_{p}a_{q}$ and $a_{rs}^{pq}\equiv a^{\dagger}_{p}a^{\dagger}_{q}a_{s}a_{r}$, the corresponding excitation operators are

[TABLE]

To compute $\mathbf{t^{(1)}}$, we consider a system of projected linear equations

[TABLE]

Using the definition of $\tilde{H}^{(1)}$ from Eq. 18, this system of equations can be expressed in the matrix form60

[TABLE]

where the zeroth-order Hamiltonian and perturbation matrix elements are defined as

[TABLE]

and $E_{0}$ is the zeroth-order (reference) energy. Eq. 32 is identical to equation that defines the first-order wavefunction in the standard Rayleigh–Schrödinger perturbation theory. Since $H^{(0)}$ is the Dyall Hamiltonian, the first-order MR-ADC reference wavefunction $\ket{\Psi^{(1)}}=T^{(1)}\ket{\Psi_{0}}$ is equivalent to the first-order wavefunction in internally-contracted second-order $N$-electron valence perturbation theory (NEVPT2).24, 25, 26 Importantly, this suggests that solutions of Eq. 31 do not suffer from intruder-state problems, provided that $\ket{\Psi_{0}}$ is the ground-state reference wavefunction. The $\mathbf{t^{(1)}}$ amplitudes can be used to compute the second-order correlation correction to the reference energy

[TABLE]

which is equivalent to the NEVPT2 correlation energy. We note that Eqs. 32 and 35 have been recently derived in the context of perturbation expansion of internally-contracted multi-reference coupled cluster theory.108

Evaluating the MR-ADC(2) matrices in Sections 2.2.1 and 2.2.1 also requires semi-internal amplitudes of the second-order excitation operator $T^{(2)}$

[TABLE]

These parameters are obtained by solving the second-order linear equations

[TABLE]

where the matrix elements of $\mathbf{V^{(2)}}$ are defined as

[TABLE]

Eq. 37 is analogous to the first-order Eq. 32 with r.h.s. modified by the second-order matrix $\mathbf{V^{(2)}}$ and, thus, can be solved in a similar way. In practice, only a small number of terms in Sections 2.2.1 and 2.2.1 depend on the $\mathbf{t^{(2)}}$ amplitudes and their contributions have a very small effect on the ionization energies and spectral intensities. We will discuss solution of the first- and second-order amplitude equations in more detail in Section 3.2.

2.2.3 Ionization Operator Manifolds

To determine the ionization operators $h_{\mu}^{(k){\dagger}}$ ($k=0,1$), we use the fact that these operators must satisfy two requirements:72, 60 (i) at the $k$th order, the particle-hole rank of $h_{\mu}^{(k){\dagger}}$ must not exceed that of $\tilde{q}_{\mu}^{(k){\dagger}}$ or $\tilde{q}_{\mu}^{(k)}$ for the forward or backward components of the propagator, respectively; (ii) $h_{\mu}^{(k){\dagger}}$ must fulfill the vacuum annihilation condition (VAC)75, 76, 77, 78 with respect to the reference state, i.e. $h_{\mu}^{(k)}\ket{\Psi_{0}}=0$, which ensures decoupling of the forward and backward components of the propagator in Section 2.1.72, 60 To obtain $h_{\mu}^{(0){\dagger}}$, we recall that $\tilde{q}_{\mu}^{(0)}=a_{p}$, where the annihilation operator can be of three different types: $a_{i}$, $a_{x}$, or $a_{a}$ (core, active, or external). Out of these three classes, only the core operator $a_{i}$ satisfies VAC with respect to $\ket{\Psi_{0}}$ ($a^{\dagger}_{i}\ket{\Psi_{0}}=0$) and, thus, can be added to $h_{\mu}^{(0){\dagger}}$. Since $\ket{\Psi_{0}}$ does not contain electrons in the active space, the external operator $a_{a}$ is redundant ($a_{a}\ket{\Psi_{0}}=0$) and cannot be included in $h_{\mu}^{(0){\dagger}}$. Although the active-space operator $a_{x}$ does not fulfill VAC ($a^{\dagger}_{x}\ket{\Psi_{0}}\neq 0$), it can be expanded60 in the form $a_{x}=\sum_{I}Z^{\dagger}_{I}c_{I,x}$, where $Z^{\dagger}_{I}$ is a complete set of active-space eigenoperators,109, 110, 111 defined as:

[TABLE]

Here, $\ket{\Psi_{I}^{N-1}}$ are the CASCI states of the ionized system with $N-1$ electrons computed using the active space and one-electron basis of the reference state $\ket{\Psi_{0}}$. We note that in the context of propagator theory the configurational operators $Z^{\dagger}_{I}$ were first used by Freed and Yeager109 and have two important properties: they are linearly-independent and include all types of active-only ionization operators ($a_{x}$, $a^{\dagger}_{x}a_{y}a_{z}$, $\ldots$). Incidentally, these operators also satisfy VAC with respect to $\ket{\Psi_{0}}$ and can be added to $h_{\mu}^{(0){\dagger}}$. Although we have assumed that the set of operators $Z^{\dagger}_{I}$ is complete, only a subset of these operators corresponding to CASCI states in the spectral region of interest need to be included in practice. We summarize that the MR-ADC(2) zeroth-order manifold $h_{\mu}^{(0){\dagger}}$ consists of two sets of operators:

[TABLE]

Following a similar strategy, we determine that the first-order operators $h_{\mu}^{(1){\dagger}}$ have a general form $a_{qr}^{p}\equiv a^{\dagger}_{p}a_{r}a_{q}$ and can be further divided into five classes

[TABLE]

describing ionization in the core or active spaces accompanied by core-active, active-external, or core-external single excitations, as shown in Figure 2. The all-active operators $a_{zy}^{x}$ do not appear in $\mathbf{h}^{(1){\dagger}}$, since they are already included in the $\mathbf{h}^{(0){\dagger}}$ manifold by the $Z^{\dagger}_{I}$ operators.

Figure 3 illustrates perturbative structure of the MR-ADC(2) effective Liouvillean ($\mathbf{M}$) and overlap ($\mathbf{S}$) matrices. The $\{\mathbf{h}^{(0){\dagger}}|\tilde{\mathcal{H}}^{(k)}|\mathbf{h}^{(0){\dagger}}\}$ block of the $\mathbf{M}$ matrix includes all contributions up to $k$ $=$ 2, while the coupling block $\{\mathbf{h}^{(1){\dagger}}|\tilde{\mathcal{H}}^{(k)}|\mathbf{h}^{(0){\dagger}}\}$ is evaluated to first order, as given by Section 2.2.1. In the manifold of first-order ionized states, the $\{\mathbf{h}^{(1){\dagger}}|\tilde{\mathcal{H}}^{(0)}|\mathbf{h}^{(1){\dagger}}\}$ sector is block-diagonal with non-zero elements for the $h_{\mu}^{(1){\dagger}}$ excitations from the same class (Eq. 41). Overall, the general perturbative structure of the MR-ADC(2) matrices closely resembles that of non-Dyson SR-ADC(2)93, 94, 95 and the two methods become equivalent in the limit of single-determinant reference wavefunction $\ket{\Psi_{0}}$.

3 Implementation

3.1 General Algorithm

In this section, we describe a general algorithm of our MR-ADC(2) implementation for complete active space (CAS) reference wavefunctions. Although in this work we always employ the ground-state CASSCF wavefunction of a neutral system as a reference, in MR-ADC other choices of reference orbitals are possible (e.g., Hartree-Fock, state-averaged, or unrestricted natural orbitals).112 The main steps of the MR-ADC(2) algorithm are summarized below:

Choose active space, compute the reference orbitals and CAS wavefunction $\ket{\Psi_{0}}$ for the neutral system with $N$ electrons. 2. 2.

Using reference orbitals, compute the CASCI energies $E_{I}^{N-1}$ and wavefunctions $\ket{\Psi_{I}^{N-1}}$ for $N_{\mathrm{CI}}$ lowest-energy states of the ionized system with $(N-1)$ electrons. 3. 3.

Compute active-space reduced density matrices (RDMs) for the reference state $\ket{\Psi_{0}}$, transition RDMs between $\ket{\Psi_{0}}$ and ionized states $\ket{\Psi_{I}^{N-1}}$, and transition RDMs between two ionized states $\ket{\Psi_{I}^{N-1}}$. 4. 4.

Solve linear amplitude Eqs. 32 and 37 to compute $\mathbf{t^{(1)}}$ and $\mathbf{t^{(2)}}$. 5. 5.

Solve the generalized eigenvalue problem (10) to obtain ionization energies $\boldsymbol{\Omega}$. 6. 6.

Compute spectroscopic amplitudes (11) and (if necessary) spectral function (13).

As discussed in Section 2.2, the number of active-space ionized states ($N_{\mathrm{CI}}$) should be sufficiently large to include all important CASCI states in the spectral region of interest. Implementation of the algorithm outlined above requires derivation of equations for contributions to the M, T, and S matrices (Sections 2.2.1, 2.2.1 and 2.2.1). Although most of these contributions have compact expressions, matrix elements of the second-order effective Hamiltonian (e.g., $\{\mathbf{h}^{(0){\dagger}}|\tilde{\mathcal{H}}^{(2)}|\mathbf{h}^{(0){\dagger}}\}$) are very complicated containing $\sim$ 250-300 terms for each matrix block. Such algebraic complexity is a common feature of many internally-contracted multi-reference theories.113, 41, 31, 59, 54

To speed up tedious derivation and implementation of MR-ADC(2), we have developed a Python program that automatically generates equations and code for arbitrary-order MR-ADC(n) approximation. Our code generator is a modified version of the SecondQuantizationAlgebra (SQA) program developed by Neuscamman and co-workers.113 We use SQA to define and normal-order all active-space creation and annihilation operators in Sections 2.2.1, 2.2.1 and 2.2.1 with respect to the physical vacuum. Next, we additionally normal-order core creation and annihilation operators relative to the Fermi vacuum and evaluate expectation values with respect to the active-space states $\ket{\Psi_{0}}$ and $\ket{\Psi_{I}^{N-1}}$. The resulting equations, written as contractions of the one- and two-electron integrals, $\mathbf{t^{(1)}}$ and $\mathbf{t^{(2)}}$ amplitudes, and RDMs, are used to generate code and can be implemented using any available tensor contraction engine. We present working equations for all matrix elements in Sections 2.2.1, 2.2.1 and 2.2.1 in the Supporting Information.

In Sections 3.2, 3.3 and 3.4, we provide more details about the solution of amplitude equations, efficient computation of terms that depend on high-order RDMs, and solution of the generalized eigenvalue problem.

3.2 Amplitude Equations

General form of the first- and second-order amplitude equations has been discussed in Section 2.2.2. Since the Dyall Hamiltonian (Eq. 4) does not contain terms that couple excitations outside of the active space, its matrix representation $\mathbf{H^{(0)}}$ (Eq. 33) is block-diagonal and the amplitude equations (32) and (37) can be solved for each block separately. Using the standard notation for classifying excitations adopted in N-electron valence perturbation theory,24, 25, 26 operators $\boldsymbol{\tau}$ in Section 2.2.2 are split into eight groups $\boldsymbol{\tau}^{\mathbf{[i]}}$ ($i$ $\in$ $\{0;+1;-1;+2;-2;+1^{\prime};-1^{\prime};0^{\prime}\}$), where $i$ is the number of electrons added to ($i>0$) or removed from ($i<0$) active space upon excitation. The operator classes with $i$ $\in$ $\{+1^{\prime};-1^{\prime};0^{\prime}\}$ are used to represent three coupled sets of single and semi-internal double excitations: $\boldsymbol{\tau}^{\mathbf{[+1^{\prime}]}}=\{a_{i}^{x};\ a_{ix}^{yz}\}$, $\boldsymbol{\tau}^{\mathbf{[-1^{\prime}]}}=\{a_{x}^{a};\ a_{xy}^{az}\}$, and $\boldsymbol{\tau}^{\mathbf{[0^{\prime}]}}=\{a_{i}^{a};\ a_{ix}^{ay}\}$.

Separating the $\mathbf{H^{(0)}}$, $\mathbf{t^{(1)}}$, and $\mathbf{V^{(1)}}$ matrices in Eq. 32 into blocks according to excitation classes $\boldsymbol{\tau}^{\mathbf{[i]}}$ (denoted as $\mathbf{K^{[i]}}$, $\mathbf{t^{[i](1)}}$, and $\mathbf{V^{[i](1)}}$, respectively), we express the first-order amplitude equations in the following form

[TABLE]

To solve Eq. 42 for each excitation class, we consider the generalized eigenvalue problem for the matrix $\mathbf{K^{[i]}}$

[TABLE]

which allows to obtain expression for the first-order amplitudes60

[TABLE]

where $K_{\mu\nu}^{[i]}=\braket{\Psi_{0}}{\tau^{[i]{\dagger}}_{\mu}(H^{(0)}-E_{0})\tau^{[i]}_{\nu}}{\Psi_{0}}$, $S^{[i]}_{\mu\nu}=\braket{\Psi_{0}}{\tau^{[i]{\dagger}}_{\mu}\tau^{[i]}_{\nu}}{\Psi_{0}}$, and $\mathbf{\tilde{Z}^{[i]}}=(\mathbf{S^{[i]}})^{1/2}\,\mathbf{Z^{[i]}}$. Computing the $\mathbf{t^{[i](1)}}$ amplitudes in Eq. 44 requires diagonalizing $\mathbf{K^{[i]}}$ and $\mathbf{S^{[i]}}$ and removing linear dependencies corresponding to eigenvectors of $\mathbf{S^{[i]}}$ with small eigenvalues. Since the matrix elements $K_{\mu\nu}^{[i]}$ and $S^{[i]}_{\mu\nu}$ are zero when the operators $\tau^{[i]{\dagger}}_{\mu}$ and $\tau^{[i]}_{\nu}$ do not share the same core and external indices, diagonalization of $\mathbf{K^{[i]}}$ and $\mathbf{S^{[i]}}$ can be performed very efficiently. For the semi-internal amplitudes $\mathbf{t^{[i](1)}}$ ($i$ $\in$ $\{+1^{\prime};-1^{\prime};0^{\prime}\}$), removing redundancies in the overlap matrix may introduce small size-consistency errors of the MR-ADC energies due to the appearance of disconnected terms in the amplitude equations that become non-zero when linear dependencies are eliminated.60, 114 To restore full size-consistency of the MR-ADC energies, we use the approach developed by Hanauer and Köhn115 that removes the disconnected terms by transforming the excitation operators $\boldsymbol{\tau}^{\mathbf{[i]}}$ ($i$ $\in$ $\{+1^{\prime};-1^{\prime};0^{\prime}\}$) to a generalized normal-ordered form. We will demonstrate size-consistency of the MR-ADC(2) ionization energies in Section 5.1.

We use Eq. 44 to compute $\mathbf{t^{[i](1)}}$ for all double ($i$ $\in$ $\{0;+1;-1;+2;-2\}$) and one class of semi-internal ($i$ $=$ $0^{\prime}$) excitations. For the $\mathbf{t^{[+1^{\prime}](1)}}$ and $\mathbf{t^{[-1^{\prime}](1)}}$ amplitudes, diagonalization of $\mathbf{K^{[+1^{\prime}]}}$ and $\mathbf{K^{[-1^{\prime}]}}$ requires the four-particle reduced density matrix (4-RDM) of the reference state $\ket{\Psi_{0}}$, which is expensive to compute and store in memory for large active spaces (see Section 3.3 for details). To avoid computation of 4-RDM, we evaluate $\mathbf{t^{[+1^{\prime}](1)}}$ and $\mathbf{t^{[-1^{\prime}](1)}}$ using imaginary-time algorithm developed in Ref. 60, which employs a Laplace transform116, 57 to evaluate the operator resolvent $(H^{(0)}-E_{0})^{-1}$ without explicitly constructing and inverting the $\mathbf{K^{[+1^{\prime}]}}$ and $\mathbf{K^{[-1^{\prime}]}}$ matrices.

The second-order amplitude equations (37) need to be solved only for the semi-internal amplitudes $\mathbf{t^{[+1^{\prime}](2)}}$, $\mathbf{t^{[-1^{\prime}](2)}}$, and $\mathbf{t^{[0^{\prime}](2)}}$ (Eq. 36). Among these, only $\mathbf{t^{[+1^{\prime}](2)}}$ enter equations for the M matrix, while all three sets of semi-internal amplitudes are necessary to compute the T matrix elements. The second-order amplitudes can be obtained in a similar way as their first-order counterparts $\mathbf{t^{[i](1)}}$, i.e. by expressing $\mathbf{t^{[i](2)}}$ in the form of Eq. 44 (with $\mathbf{V^{[i](1)}}$ replaced by $\mathbf{V^{[i](2)}}$ defined in Eq. 38) or using the imaginary-time algorithm. Although solving the second-order equations is straightforward, matrix elements of the perturbation operator $\mathbf{V^{[i](2)}}$ contain $\sim$ 600 terms and are rather tedious to evaluate. On the other hand, since the primary role of $\mathbf{t^{[i](2)}}$ ($i$ $\in$ $\{+1^{\prime};-1^{\prime};0^{\prime}\}$) is to describe relaxation of the orbitals, their contributions are expected to have a small effect on the results of the MR-ADC(2) method that already incorporates orbital relaxation via the first-order amplitudes $\mathbf{t^{[i](1)}}$ and ionization operators $\mathbf{h}^{(1){\dagger}}$. To test this, we considered an approximation where we neglect contributions of $\mathbf{t^{[+1^{\prime}](2)}}$ and $\mathbf{t^{[-1^{\prime}](2)}}$ and approximate $\mathbf{t^{[0^{\prime}](2)}}$ by setting $t_{ix}^{ay(2)}\approx 0$ and neglecting all terms that depend on active-space RDMs in $\mathbf{V^{[0^{\prime}](2)}}$ to obtain $t_{i}^{a(2)}$ (see the Supporting Information). The resulting amplitude equations ensure that MR-ADC(2) is equivalent to SR-ADC(2) in the single-reference limit. As demonstrated in the Supporting Information, approximating the $\mathbf{t^{(2)}}$ terms has a very small effect on the MR-ADC(2) results with errors of $\leq$ 0.005 eV and $\leq$ $3\times 10^{-4}$ in ionization energies and spectroscopic factors, respectively. For this reason, we adopted this approximation in our implementation of MR-ADC(2).

3.3 Avoiding High-Order Reduced Density Matrices

As other internally-contracted multi-reference perturbation theories, MR-ADC(2) contains terms that depend on high-order reduced density matrices (e.g., 4-RDM) in its equations. In this section, we will demonstrate that these terms can be efficiently evaluated without computing and storing 4-RDMs in memory. There are two sources of high-order RDMs in the MR-ADC(2) equations: (i) $\mathbf{t^{(1)}}$ and $\mathbf{t^{(2)}}$ amplitude equations and (ii) second-order contributions to the effective Liouvillean matrix $\mathbf{M}$. As discussed in Section 3.2, using the imaginary-time algorithm60 allows to completely avoid computation of 4-RDM in the amplitude equations.

For the $\mathbf{M}$ matrix, 4-RDMs appear in expectation values of the second-order effective Hamiltonian $\tilde{\mathcal{H}}^{(2)}$ with respect to the reference ($\braket{\Psi_{0}}{\tilde{\mathcal{H}}^{(2)}}{\Psi_{0}}$) and ionized ($\braket{\Psi_{I}^{N-1}}{\tilde{\mathcal{H}}^{(2)}}{\Psi_{J}^{N-1}}$) wavefunctions. In particular, the latter matrix elements depend on transition 4-RDMs between all CASCI ionized states (e.g., $\braket{\Psi_{I}^{N-1}}{a^{\dagger}_{w}a^{\dagger}_{x}a^{\dagger}_{y}a^{\dagger}_{z}a_{z^{\prime}}a_{y^{\prime}}a_{x^{\prime}}a_{w^{\prime}}}{\Psi_{J}^{N-1}}$), which have a high $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}^{2}N^{8}_{\mathrm{act}})$ computational scaling, where $N_{\mathrm{det}}$ is the dimension of CAS Hilbert space, $N_{\mathrm{CI}}$ is the number of CASCI ionized states, and $N_{\mathrm{act}}$ is the number of active orbitals. To demonstrate how to avoid computation of 4-RDMs, we consider one of the contributions to the $\braket{\Psi_{I}^{N-1}}{\tilde{\mathcal{H}}^{(2)}}{\Psi_{J}^{N-1}}$ matrix elements

[TABLE]

where we use shorthand notation for the first-order amplitudes $t_{xy}^{az(1)}\equiv t_{xy}^{az}$ and CASCI states $\ket{\Psi_{I}^{N-1}}\equiv\ket{\Psi_{I}}$. Changing the order of creation and annihilation operators, we express Eq. 45 in the following form

[TABLE]

where the remaining terms involve contractions of transition 2- and 3-RDMs. Computing intermediate states

[TABLE]

we evaluate the first term in Eq. 46 using a compact expression

[TABLE]

Using Eqs. 47, 48 and 49 allows us to significantly lower the cost of computing transition 4-RDM terms from $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}^{2}N^{8}_{\mathrm{act}})$ to $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}^{2}N^{3}_{\mathrm{act}}N_{\mathrm{ext}})$, where $N_{\mathrm{ext}}$ is the number of external orbitals. We use the same technique to efficiently evaluate all 4-RDM terms that appear in the $\braket{\Psi_{I}^{N-1}}{\tilde{\mathcal{H}}^{(2)}}{\Psi_{J}^{N-1}}$ and $\braket{\Psi_{0}}{\tilde{\mathcal{H}}^{(2)}}{\Psi_{0}}$ matrix elements. We note that similar techniques have been used to avoid computation of 4-RDM in implementations of complete active space second-order perturbation theory (CASPT2) and NEVPT2 in combination with matrix product state wavefunctions.116, 117, 57

The $\mathbf{M}$ matrix elements also depend on transition RDMs of the form $\braket{\Psi_{0}}{a^{\dagger}_{w}a^{\dagger}_{x}a^{\dagger}_{y}a^{\dagger}_{z}a_{z^{\prime}}a_{y^{\prime}}a_{x^{\prime}}}{\Psi_{I}^{N-1}}$, which we denote as 3.5-RDMs. These RDMs contribute to the second-order matrix elements $\braket{\Psi_{0}}{a^{\dagger}_{i}\tilde{\mathcal{H}}^{(2)}}{\Psi_{I}^{N-1}}$, as well as some elements of the first-order off-diagonal blocks $\{\mathbf{h}^{(1){\dagger}}|\tilde{\mathcal{H}}^{(1)}|\mathbf{h}^{(0){\dagger}}\}$ and $\{\mathbf{h}^{(0){\dagger}}|\tilde{\mathcal{H}}^{(1)}|\mathbf{h}^{(1){\dagger}}\}$ in Section 2.2.1. For example, a 3.5-RDM contribution to $\braket{\Psi_{0}}{a^{\dagger}_{i}\tilde{\mathcal{H}}^{(2)}}{\Psi_{I}^{N-1}}$ has a form

[TABLE]

To evaluate this term, we reorder creation and annihilation operators, contract $v^{xy}_{zw}$ and $t^{vu^{\prime}}_{aw^{\prime}}$ with $a^{\dagger}_{x}a^{\dagger}_{y}a_{w}\rvert\Psi_{0}\rangle$ and $a^{\dagger}_{v}a^{\dagger}_{u^{\prime}}a_{w^{\prime}}\rvert\Psi_{I}\rangle$ to form intermediate states ($\rvert v^{z}\rangle$ and $\rvert t_{I}^{a}\rangle$), and contract $t^{iz}_{au}$ with their overlap matrix element ($\langle v^{z}\lvert a^{\dagger}_{u}\rvert t^{a}_{I}\rangle$). As in the case of 4-RDM, using intermediate states allows to completely avoid computation and storage of 3.5-RDMs for all terms of the $\mathbf{M}$ matrix, lowering computational scaling from $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}N^{7}_{\mathrm{act}})$ to $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}N^{2}_{\mathrm{act}}N_{\mathrm{ext}})$.

Combining efficient algorithms for the solution of amplitude equations and evaluation of high-order RDM terms, our MR-ADC(2) implementation has $\mathcal{O}(N_{\mathrm{det}}N^{2}_{\mathrm{CI}}N^{6}_{\mathrm{act}})$ computational scaling, which is significantly lower than the $\mathcal{O}(N_{\mathrm{det}}N^{8}_{\mathrm{act}})$ scaling of the conventional multi-reference perturbation theories (e.g., CASPT2 or NEVPT2) with the number of active orbitals. Although the scaling of our current MR-ADC(2) algorithm originates from computing transition 3-RDMs ($\braket{\Psi_{I}^{N-1}}{a^{\dagger}_{w}a^{\dagger}_{x}a^{\dagger}_{y}a_{y^{\prime}}a_{x^{\prime}}a_{w^{\prime}}}{\Psi_{J}^{N-1}}$) for all ionized states, we note that using intermediate states the computational cost can be further lowered to $\mathcal{O}(N_{\mathrm{det}}N_{\mathrm{CI}}N^{6}_{\mathrm{act}})$. We did not take advantage of it in our present implementation.

3.4 Solution of the Generalized Eigenvalue Problem

Finally, we briefly discuss solution of the MR-ADC(2) generalized eigenvalue problem in Eq. 10. Since the $\mathbf{M}$ and $\mathbf{S}$ matrices are computed in the non-orthogonal basis of internally-contracted ionized states, we transform the eigenvalue equation to the symmetrically-orthogonalized form

[TABLE]

where $\mathbf{\tilde{M}}=\mathbf{S}^{-1/2}\mathbf{M}\mathbf{S}^{-1/2}$ and $\mathbf{\tilde{Y}}=\mathbf{S}^{1/2}\mathbf{Y}$. Here, the overlap matrix $\mathbf{S}$ contains four non-diagonal blocks corresponding to ionized states $\ket{\Psi_{\mu}}=\{a_{ij}^{x}\ket{\Psi_{0}}$; $a_{ix}^{a}\ket{\Psi_{0}}$; $a_{xy}^{a}\ket{\Psi_{0}}$; $a_{i}\ket{\Psi_{0}}$; $a_{ix}^{y}\ket{\Psi_{0}}\}$ (Figure 3b). Conveniently, the $\mathbf{S}^{-1/2}$ matrix can be constructed together with the $(\mathbf{S^{[i]}})^{-1/2}$ matrices used for solution of the amplitude equations (Section 3.2). As an example, we consider non-zero elements of $\mathbf{S}$ for $a_{ij}^{x}\ket{\Psi_{0}}$ that have the form $S_{ijx,ijy}=\braket{\Psi_{0}}{a_{y}^{ij}a_{ij}^{x}}{\Psi_{0}}=\braket{\Psi_{0}}{a_{y}a^{\dagger}_{x}}{\Psi_{0}}$. These elements are equal to the $\mathbf{S^{[+1]}}$ matrix elements $S_{ijay,ijax}^{[+1]}=\braket{\Psi_{0}}{a_{ay}^{ij}a_{ij}^{ax}}{\Psi_{0}}=\braket{\Psi_{0}}{a_{y}a^{\dagger}_{x}}{\Psi_{0}}$. Thus, by diagonalizing the density matrix $\braket{\Psi_{0}}{a_{y}a^{\dagger}_{x}}{\Psi_{0}}$ and removing linearly-dependent eigenvectors corresponding to small eigenvalues ($<$ $\eta_{d}$, where $\eta_{d}$ is a user-defined truncation parameter), we simultaneously obtain elements of $(\mathbf{S^{[+1]}})^{-1/2}$ and $\mathbf{S}^{-1/2}$ for the $a_{ij}^{x}\ket{\Psi_{0}}$ ionized wavefunctions. Similarly, we construct $(\mathbf{S^{[-1]}})^{-1/2}$ and $(\mathbf{S^{[-2]}})^{-1/2}$ together with $\mathbf{S}^{-1/2}$ for $a_{ix}^{a}\ket{\Psi_{0}}$ and $a_{xy}^{a}\ket{\Psi_{0}}$, respectively.

For the $a_{ij}^{x}\ket{\Psi_{0}}$, $a_{ix}^{a}\ket{\Psi_{0}}$, and $a_{xy}^{a}\ket{\Psi_{0}}$ states, numerical instabilities due to linear dependencies are completely eliminated when using small truncation parameters ($\eta_{d}$ $\sim$ $10^{-10}$). Except for very small active spaces ($N_{\mathrm{act}}<6$), orthogonalization of these ionized states does not require discarding any eigenvectors of the overlap matrix. The zeroth-order $a_{i}\ket{\Psi_{0}}$ and first-order $a_{ix}^{y}\ket{\Psi_{0}}$ ionized states exhibit much stronger linear dependencies in their overlap matrix. To remove these linear dependencies, we project out $a_{i}\ket{\Psi_{0}}$ from $a_{ix}^{y}\ket{\Psi_{0}}$ using the projection approach developed by Hanauer and Köhn114 and subsequently orthogonalize $a_{ix}^{y}\ket{\Psi_{0}}$ between each other. Importantly, this ensures that the zeroth-order states $a_{i}\ket{\Psi_{0}}$, which are already orthogonal, are not affected by removing redundancies in the first-order $a_{ix}^{y}\ket{\Psi_{0}}$ ionization manifold. To discard linearly-dependent eigenvectors of the $a_{ix}^{y}\ket{\Psi_{0}}$ overlap matrix, we use a larger truncation parameter ($\eta_{s}$ $\sim$ $10^{-6}$) than the one used for other ionized states ($\eta_{d}$).

We solve the eigenvalue problem (51) using a multi-root implementation of the Davidson algorithm,118, 119 which avoids storing the full M and S matrices, significantly reducing the memory requirements. Since the second-order block $\{\mathbf{h}^{(0){\dagger}}|\tilde{\mathcal{H}}^{(2)}|\mathbf{h}^{(0){\dagger}}\}$ of M is small (with $(N_{\mathrm{CI}}+N_{\mathrm{act}})^{2}$ elements) and its computation is the most time-consuming step of the MR-ADC(2) implementation, we precompute this block, store it memory, and use it for the efficient evaluation of matrix-vector products in the Davidson procedure.

4 Computational Details

We implemented MR-ADC(2) for photoelectron spectra in our pilot code Prism, which was interfaced with Pyscf120 to obtain integrals and CASCI/CASSCF reference wavefunctions. Our implementation follows the general algorithm outlined in Section 3.1. All MR-ADC(2) computations used the CASSCF reference wavefunctions with molecular orbitals optimized for the ground electronic state of each (neutral) system. To remove linear dependencies in the solution of amplitude equations and generalized eigenvalue problem, we truncated eigenvectors of the overlap matrices using two parameters: $\eta_{s}$ = $10^{-6}$ and $\eta_{d}$ = $10^{-10}$ (see Section 3.4 for details). The $\eta_{s}$ parameter was used to orthogonalize the $a_{ix}^{y}\ket{\Psi_{0}}$ ionized states and to compute the semi-internal $\mathbf{t^{[i](1)}}$ ($i$ $\in$ $\{+1^{\prime};-1^{\prime};0^{\prime}\}$) amplitudes (Section 3.2), while $\eta_{d}$ was employed for other amplitudes and ionized states. To efficiently compute $\mathbf{t^{[+1^{\prime}](1)}}$ and $\mathbf{t^{[-1^{\prime}](1)}}$, our implementation used imaginary-time algorithm,60, 116, 57 where propagation in imaginary time was performed using the embedded Runge-Kutta method that automatically determines time step based on the accuracy parameter $\Delta_{it}$.121 In all computations, we used $\Delta_{it}$ = $10^{-7}$ $E_{\mathrm{h}}$, which allows to obtain very accurate amplitudes and reference NEVPT2 correlation energy. All MR-ADC(2) results were converged with respect to the number of CASCI ionized states ($N_{\mathrm{CI}}$). For most of the systems employed in this study, using $N_{\mathrm{CI}}$ = 20 was enough to obtain well-converged results.

We benchmarked the accuracy of MR-ADC(2) for a set of small molecules (\ceHF, \ceF2, \ceCO, \ceN2, \ceH2O, \ceCS, \ceH2CO, and \ceC2H4), carbon dimer (\ceC2), and hydrogen chains (\ceH10 and \ceH30). For small molecules, equilibrium and stretched geometries were considered. The equilibrium structures were taken from Ref. 94. For diatomic molecules, the stretched geometries were obtained by increasing the bond length by a factor of two. For the \ceH2O, \ceH2CO, and \ceC2H4 stretched geometries, we doubled the \ceO-H, \ceC-O, and \ceC-C bond distances, respectively. The \ceC-C bond length in \ceC2 was set to 1.2425 Å, which is very close to its equilibrium geometry. Unless noted otherwise, all computations employed the aug-cc-pVDZ basis set.122 For \ceH2CO and \ceC2H4, the cc-pVDZ basis set was used for the hydrogen atoms, as employed in Ref. 94. We denote active spaces used in CASCI/CASSCF as ($n$e, $m$o), where $n$ is the number of active electrons and $m$ is the number of active orbitals. Active spaces of small molecules included 10 orbitals with $n$ = 8, 14, 10, 10, 8, 10, 12, and 10 active electrons for \ceHF, \ceF2, \ceCO, \ceN2, \ceH2O, \ceCS, \ceH2CO, and \ceC2H4, respectively. For \ceC2, the (8e, 12o) active space was used. For the hydrogen chains, we employed the (10e, 10o) active space.

The MR-ADC(2) results were compared to results of single-reference non-Dyson ADC methods (SR-ADC(2) and SR-ADC(3)),93, 94, 95 equation-of-motion coupled cluster theory for ionization energies with single and double excitations (EOM-CCSD),123, 124, 83 quasi-degenerate strongly-contracted second-order N-electron valence perturbation theory (QD-NEVPT2),26 as well as full configuration interaction (FCI). All methods employed the same geometries and basis sets as those used for MR-ADC(2). SR-ADC(2) and SR-ADC(3) were implemented by our group as a module in the development version of Pyscf. The FCI results were computed using the semistochastic heat-bath configuration interaction algorithm (SHCI) implemented in the Dice program.12, 13, 14 The SHCI electronic energies were extrapolated using a linear fit according to procedure described in Ref. 14. We estimate that errors of the computed SHCI energy differences relative to FCI do not exceed 0.03 eV. For \ceH2CO and \ceC2H4, the $1s$ atomic orbitals of carbon and oxygen were not correlated in the SHCI computations. For all other methods, all electrons were correlated in all computations. The EOM-CCSD and QD-NEVPT2 results were obtained using Q-Chem125 and Orca,126 respectively. For the ground state of each neutral system, QD-NEVPT2 used the same active spaces and CASSCF reference wavefunctions as those employed in MR-ADC(2). The QD-NEVPT2 computations of ionized states used the state-averaged CASSCF reference wavefunctions, where state-averaging included four electronic states for each abelian subgroup irreducible representation of the full symmetry point group.

Intensities of photoelectron transitions were characterized by computing spectroscopic factors

[TABLE]

where $X_{p,\mu}$ are elements of the spectroscopic amplitude matrix $\mathbf{X}_{\pm}$ defined in Eq. 11. Spectroscopic factors in Eq. 52 correspond to intensities of photoelectron transitions under the approximation that only single-electron detachment contributes to the spectrum. More rigorous simulation of photoelectron intensities require computation of Dyson orbitals with explicit treatment of the wavefunction of injected free electron and will be one of the subjects of our future work.127

5 Results

5.1 Size-Consistency of Energies and Properties

We begin by testing size-consistency of the MR-ADC(2) ionization energies and spectroscopic factors. As for single-reference ADC, the MR-ADC equations are fully connected, which guarantees size-consistency of the MR-ADC energies and transition properties. In practice, however, removing redundancies in the overlap matrix during the solution of the MR-ADC amplitude equations may result in small size-consistency errors.60 As we discussed in Section 3.2, in this work we employ a technique developed by Hanauer and Köhn that restores size-consistency of the MR-ADC results. Table 1 shows deviations from size-consistency of the MR-ADC(2) ionization energies ($\Delta\Omega$) and spectroscopic factors ($\Delta P$) for the \ce(H2O)2 and \ce(HF)2 systems, each composed of two noninteracting monomers with near-equilibrium ($r_{e}$) and stretched geometries ($2\times r_{e}$). The computed size-consistency errors are very small: $\Delta\Omega$ $\sim$ $10^{-5}$ eV and $\Delta P$ $\sim$ $10^{-6}$ on average, with the largest errors of $\Delta\Omega$ = 1.2 $\times$ $10^{-4}$ eV and $\Delta P$ = $-$4.5 $\times$ $10^{-6}$. These remaining errors originate from a finite time step used in the imaginary-time algorithm for solving the semi-internal amplitude equations and become increasingly smaller with a tighter $\Delta_{it}$ parameter (see Section 4 for details). Overall, our numerical results demonstrate size-consistency of the MR-ADC(2) results in the present implementation.

5.2 Small Molecules

In this section, we benchmark the MR-ADC(2) accuracy for predicting ionization energies of small molecules. Table 2 compares vertical ionization energies ($\Omega$) and spectroscopic factors ($P$) of MR-ADC(2) with those obtained by single-reference non-Dyson ADC methods (SR-ADC), equation-of-motion coupled cluster theory with single and double excitations (EOM-CCSD), quasi-degenerate NEVPT2 (QD-NEVPT2), and full configuration interaction (FCI) for a set of eight molecules near their equilibrium geometries (see Section 4 for computational details). In addition to strict second- and third-order SR-ADC (SR-ADC(2) and SR-ADC(3)), Table 2 also presents results of SR-ADC(3) incorporating high-order self-energy corrections, reported in Ref. 94, which we denote as SR-ADC(3+). Out of six approximate methods, the best agreement with FCI is shown by SR-ADC(3+), EOM-CCSD, and QD-NEVPT2. These three methods produce similar mean absolute errors in vertical ionization energies ($\Delta_{\mathrm{MAE}}$ $\sim$ 0.2 eV) with standard deviations from the mean signed error ($\Delta_{\mathrm{STD}}$) ranging from $\sim$ 0.15 to 0.3 eV, as illustrated in Figure 4a. The MR-ADC(2) method shows a similar $\Delta_{\mathrm{STD}}$ error (0.23 eV), but a larger $\Delta_{\mathrm{MAE}}$ error (0.56 eV), which is lower than $\Delta_{\mathrm{MAE}}$ of SR-ADC(2) (0.83 eV), but higher than that of SR-ADC(3) (0.30 eV), indicating that including high-order effects in MR-ADC(2) improves its accuracy relative to SR-ADC(2). For all systems, the MR-ADC(2) ionization energies systematically overestimate energies computed using FCI, showing a good agreement with FCI for energy spacings between electronic states of the ionized systems ($\Delta_{\mathrm{MAE}}$ of 0.11 eV and $\Delta_{\mathrm{STD}}$ of 0.10 eV). The QD-NEVPT2 method shows the best agreement with FCI for energy spacings ($\Delta_{\mathrm{MAE}}$ and $\Delta_{\mathrm{STD}}$ of 0.03 eV), while EOM-CCSD shows larger errors compared to MR-ADC(2) ($\Delta_{\mathrm{MAE}}$ = 0.16 eV, $\Delta_{\mathrm{STD}}$ = 0.27 eV). The MR-ADC(2) spectroscopic factors agree well with those computed using SR-ADC(3) and SR-ADC(3+), with two exceptions observed for the $6\sigma$ state of \ceCS and the $1b_{2}$ state of \ceH2CO. In these cases, the computed spectroscopic factors vary significantly depending on the order of the ADC approximation, suggesting that properties of these photoelectron transitions are significantly affected by electron correlation effects.

To assess performance of MR-ADC(2) when strong correlation is important, we computed its ionization energies and spectroscopic factors for molecules with stretched geometries, where at least one of the bonds is elongated by a factor of two (see Section 4 for details). The MR-ADC(2) results are shown in Table 3, along with those computed using SR-ADC(2), SR-ADC(3), EOM-CCSD, QD-NEVPT2, and FCI. Due to the difficulty of obtaining the FCI energies, we show results only for a few lowest-energy transitions of six molecules. Importance of strong electron correlation for these non-equilibrium geometries is demonstrated by the poor performance of SR-ADC(2) and SR-ADC(3), which show very large deviations from the FCI reference values with $\Delta_{\mathrm{MAE}}$ $>$ 2.5 eV and $\Delta_{\mathrm{STD}}$ $>$ 3 eV. Although SR-ADC(3) shows moderate $\sim$ 0.5 eV errors for single-bond stretching in HF and \ceF2, these errors drastically increase when multiple bonds are elongated, leading to unphysical values of ionization energies that significantly underestimate the FCI results. EOM-CCSD significantly improves prediction of ionization energies over SR-ADC(2) and SR-ADC(3), but still exhibits large errors ($\Delta_{\mathrm{MAE}}$ = 0.56 and $\Delta_{\mathrm{STD}}$ = 0.81 eV, Figure 4b). MR-ADC(2) shows performance similar to that for equilibrium geometries (Table 2), with $\Delta_{\mathrm{MAE}}$ (0.50 eV) and $\Delta_{\mathrm{STD}}$ (0.36 eV) smaller than the corresponding errors for the single-reference methods. The best agreement with FCI is again shown by QD-NEVPT2 with $\Delta_{\mathrm{MAE}}$ = 0.18 eV and $\Delta_{\mathrm{STD}}$ = 0.23 eV. As for the equilibrium geometries, the MR-ADC(2) ionization energies for stretched geometries are systematically overestimated relative to FCI, reproducing energy spacings between ionized states within 0.1 eV for all systems except \ceH2CO, where $\sim$ 0.5 eV errors are observed. QD-NEVPT2 shows a similar performance to MR-ADC(2) for energy spacings with a large error of $\sim$ 0.7 eV for the difference of the $1b_{1}$ and $5a_{1}$ ionization energies of \ceH2CO.

An important advantage of MR-ADC(2) over conventional multi-reference perturbation theories (such as QD-NEVPT2) is that it provides efficient access to spectroscopic properties. We demonstrate this by computing the photoelectron spectrum of \ceC2H4 at equilibrium and stretched geometries in the range between 8.5 and 20 eV, shown in Figure 5. The spectrum at equilibrium geometry exhibits five very intense well-separated peaks corresponding to vertical ionizations in five highest occupied molecular orbitals. All of the computed peaks are systematically shifted by $\sim$ 0.5 eV relative to FCI. The computed spacings between the main peaks are in a good agreement with FCI (Table 2), as well as experimental photoelectron spectrum.128 At the stretched geometry, the MR-ADC(2) spectrum shows four main peaks with significantly decreased intensities, along with several satellite peaks originating from shake-up transitions that involve ionization and simultaneous excitation in the valence orbitals.

5.3 Carbon Dimer

Next, we investigate performance of MR-ADC(2) for simulating photoelectron spectrum of \ceC2, which is a challenging test for ab initio methods, since electronic states of both \ceC2 and \ceC2+ require very accurate description of static and dynamic correlation.129, 130, 131, 132, 133, 14, 134, 135, 136, 137, 138, 139 Table 4 compares results of SR-ADC(3), MR-ADC(2), and QD-NEVPT2 with those from FCI. The MR-ADC(2) photoelectron spectrum, shown in Figure 6, exhibits two very intense peaks for ionizations in the $1\pi_{u}$ and $2\sigma_{u}$ orbitals, corresponding to the $1^{2}\Pi_{u}$ and $1^{2}\Sigma_{u}^{+}$ electronic states of \ceC2+, respectively. For both peaks, MR-ADC(2) is in a good agreement with FCI, showing errors in vertical ionization energies (0.16 and 0.24 eV) within $\Delta_{\mathrm{MAE}}$ and $\Delta_{\mathrm{STD}}$ of small molecules computed in Section 5.2. The MR-ADC(2) results show significant improvement over SR-ADC(3), which underestimates the $1\pi_{u}$ and $2\sigma_{u}$ ionization energies from FCI by 0.65 and 1.14 eV, respectively, indicating that description of multi-reference effects is important for these ionization processes. The best agreement with FCI is demonstrated by QD-NEVPT2, with errors smaller than 0.1 eV.

In addition to the intense peaks, the \ceC2 photoelectron spectrum also exhibits several much weaker (satellite) peaks, which involve ionization in the $1\pi_{u}$ orbital accompanied by single and double $(1\pi_{u})^{3}$ $\rightarrow$ $(3\sigma_{g})^{0}$ excitations (Table 4). Out of four satellite transitions, only two are predicted by SR-ADC(3), with large errors ($>$ 2 eV). For the singly-excited shake-up states of \ceC2+ ($1^{2}\Delta_{g}$, $1^{2}\Sigma_{g}^{-}$, and $1^{2}\Sigma_{g}^{+}$), the largest MR-ADC(2) error is 0.37 eV. However, for the doubly-excited $2^{2}\Pi_{u}$ state, MR-ADC(2) produces a larger 0.66 eV error. The QD-NEVPT2 ionization energies for all four electronic states are within 0.1 eV from the reference FCI values. The large error of MR-ADC(2) for $2^{2}\Pi_{u}$ may be attributed to the importance of differential dynamic correlation effects between this state and the ground state of \ceC2, since in MR-ADC(2) the first-order amplitudes of the effective Hamiltonian are preferentially determined for the latter state (Section 2.2.2), while in QD-NEVPT2 the first-order wavefunction is constructed for each electronic state separately. The description of these differential correlation effects is expected to improve for higher-order MR-ADC approximations and will be a subject of our future research.

5.4 Hydrogen Chains

Finally, we use MR-ADC to study equally-spaced hydrogen chains \ceH10 and \ceH30. Hydrogen chains are one-dimensional models for understanding strong electron correlation in molecules and materials, as well as the hydrogen phase diagram at high pressures.140, 141, 142, 143, 144, 145, 146 An important property of a hydrogen chain is its band gap, which can be calculated as the difference between ionization potential and electron affinity. For equally-spaced chains in the thermodynamic limit, this band gap is believed to be zero at short \ceH-H distances ($r$), corresponding to a metallic phase, and non-zero for long distances, corresponding to an insulator. Recently, Ronca et al. computed local density of states (LDOS) of the Hn chains ($n$ = 10, 30, and 50) at the central hydrogen atom using density matrix renormalization group (DMRG) method with the minimal STO-6G basis set.147 They demonstrated that for near-equilibrium and stretched geometries ($r$ = 1.8 and 3.6 $a_{0}$) LDOS converges to thermodynamic limit already for \ceH50, while for compressed chains ($r$ = 1.4 $a_{0}$) finite size effects are still significant. Although in this study all valence electrons of hydrogen atoms were correlated, importance of dynamic correlation effects beyond those in the minimal one-electron basis was not investigated.

Here, we use MR-ADC to study effect of dynamic correlation and basis set on the density of occupied states in \ceH10 and \ceH30. Figure 7 shows LDOS of \ceH10 for $r$ = 1.4, 1.8, and 3.6 $a_{0}$ computed at the central hydrogen atom using the MR-ADC(0) and MR-ADC(2) methods. We use the full valence (10e, 10o) active space for the CASSCF reference wavefunction and combine MR-ADC with the STO-6G148 and cc-pVTZ basis sets, plotting LDOS for two broadening parameters: 0.05 $E_{\mathrm{h}}$ (Figure 7) and 0.003 $E_{\mathrm{h}}$ (Supporting Information). For the minimal STO-6G basis, results of MR-ADC(0) and MR-ADC(2) are equivalent to FCI. The computed LDOS are in a very good agreement with LDOS obtained by Ronca et al. employing the dynamical DMRG algorithm for all three geometries.147 Next, we consider LDOS computed using MR-ADC(0) with the larger cc-pVTZ basis set. Increasing the basis set shifts LDOS to higher ionization energies, relative to LDOS from FCI/STO-6G. For short bond distances ($r$ = 1.4 and 1.8 $a_{0}$), the largest shifts are observed for the lowest-energy peaks corresponding to the ionization potential of the system ($\sim$ 0.03 and 0.05 $E_{\mathrm{h}}$, respectively). For the stretched chain ($r$ = 3.6 $a_{0}$), increasing the basis set compresses LDOS and shifts the position of its maximum by $\sim$ 0.04 $E_{\mathrm{h}}$. Incorporating dynamic correlation effects from MR-ADC(0) to MR-ADC(2) shifts the computed LDOS further to lower energies. For most of the peaks at short bond distances, the computed shifts are $\leq$ 0.02 $E_{\mathrm{h}}$. For $r$ = 3.6 $a_{0}$, including dynamic correlation does not significantly change position of the first band in the spectrum. Overall, our results suggest that increasing the one-electron basis set and incorporating dynamic correlation effects are similarly important and should be both taken into account in accurate computations of LDOS for hydrogen chains.

An attractive feature of MR-ADC is that it is not limited to describing ionization processes only in active orbitals. We demonstrate this by computing total density of occupied states (DOS) for the \ceH30 chain using MR-ADC(2) with the (10e, 10o) active space. Since for this system we do not include all valence orbitals in the active space, we do not consider the stretched $r$ = 3.6 $a_{0}$ geometry. Figure 8a shows MR-ADC(2) DOS computed using the STO-6G and cc-pVDZ basis sets. For both geometries, DOS computed using MR-ADC(2) with the STO-6G basis closely resembles LDOS of the same system from the DMRG study of Ronca et al.147 Figure 8b plots contributions to MR-ADC(2)/STO-6G DOS from core and active orbitals separately. For the compressed chain ($r$ = 1.4 $a_{0}$), contributions from active orbitals dominate the low-energy part of the spectrum, whereas, for equilibrium geometry ($r$ = 1.8 $a_{0}$), core and active orbitals have similar contributions to DOS already for low ionization energies. Increasing the basis set from STO-6G to cc-pVDZ shifts peaks in DOS to higher energies. As for the \ceH10 chain, the largest shifts are observed for the peak at the first ionization potential.

6 Conclusions

We presented derivation and implementation of second-order multi-reference algebraic diagrammatic construction theory (MR-ADC(2)) for simulating ionization energies and transition properties of strongly correlated systems. In MR-ADC(2), ionization energies and spectral properties are determined from poles and residues of the one-electron Green’s function that is evaluated to second order in multi-reference perturbation theory with respect to a complete active space (CAS) reference wavefunction. In contrast to conventional second-order multi-reference perturbation theories (such as multi-state CASPT2 or NEVPT2), MR-ADC(2) describes ionization in all orbitals (e.g., core and active), does not require using state-averaged wavefunctions to compute higher-energy ionized states, and provides direct access to spectroscopic properties. Although equations of MR-ADC(2) depend on four-particle reduced density matrices, we demonstrated that computation of these large matrices can be completely avoided by constructing efficient intermediates, without introducing any approximations. The resulting MR-ADC(2) implementation has a lower $\mathcal{O}(N_{\mathrm{det}}N^{6}_{\mathrm{act}})$ computational scaling with respect to the number of active orbitals ($N_{\mathrm{act}}$), compared to the $\mathcal{O}(N_{\mathrm{det}}N^{8}_{\mathrm{act}})$ scaling of conventional multi-reference perturbation theories.

We benchmarked accuracy of MR-ADC(2) for predicting ionization energies of eight small molecules, carbon dimer (\ceC2), and hydrogen chains (\ceH10 and \ceH30), against results from full configuration interaction (FCI). For small molecules, MR-ADC(2) shows consistent performance for equilibrium and stretched geometries, with mean absolute errors of $\sim$ 0.5 eV in ionization energies and 0.1 eV in energy separations between ionized states. For \ceC2, MR-ADC(2) predicts energies of the main and singly-excited satellite peaks within 0.4 eV from the FCI reference values, but has a large $\sim$ 0.7 eV error for the doubly-excited satellite transition. The QD-NEVPT2 method shows smaller ($\sim$ 0.1 eV) errors than MR-ADC(2) for all ionized states of \ceC2, providing an improved description of differential dynamic correlation effects, which are important for this system. We expect that these effects will be better described using the higher-order MR-ADC approximations, which will be one of the directions of our future work. Finally, we used MR-ADC(2) to investigate density of occupied states (DOS) in \ceH10 and \ceH30. For \ceH10, our results provide numerical evidence that including dynamic correlation effects beyond those incorporated in a full valence CAS and increasing single-particle basis set have a similar effect on the computed local DOS. Since dynamic correlation is a local phenomenon, we expect that its effect will be similar for longer hydrogen chains as well. For \ceH30, we showed that DOS computed using MR-ADC(2) combined with a small (10e, 10o) active space is in a very good agreement with previously reported results from density matrix renormalization group, incorporating 30 electrons and orbitals in the active space.

Overall, our results suggest that MR-ADC is a promising theoretical approach for computing ionization energies and spectral densities of multi-reference systems and encourage its further development. Future work will be directed towards more efficient implementation of MR-ADC(2) for systems with a large number of electrons and active orbitals, as well as the development of more accurate MR-ADC approximations that will incorporate description of higher-order dynamic correlation effects. We also plan extending our MR-ADC methods to simulations of core-level ionizations in X-ray photoelectron spectroscopy, which has become a widely used tool for experimental investigations of molecules and materials.

\acknowledgement

This work was supported by start-up funds provided by the Ohio State University. Computations were performed at the Ohio Supercomputer Center under projects PAS1317 and PAS1583.149 The authors would like to thank Sandeep Sharma for useful suggestions about carrying out SHCI computations.

\suppinfo

Comparison of MR-ADC(2) results with exact and approximate second-order amplitudes. Equations of the MR-ADC(2) method for the $\mathbf{M}$, $\mathbf{T}$, and $\mathbf{S}$ matrices. Density of states for \ceH10 with a small broadening parameter.