# Excitation energies from relativistic coupled-cluster theory of general excitation rank: Initial implementation and application to the silicon atom and to the molecules

XH (X = As, Sb, Bi)

Mickaël Hubert, Lasse K. Sørensen, Jeppe Olsen, Timo Fleig

► To cite this version:

Mickaël Hubert, Lasse K. Sørensen, Jeppe Olsen, Timo Fleig. Excitation energies from relativistic coupled-cluster theory of general excitation rank: Initial implementation and application to the silicon atom and to the molecules XH (X = As, Sb, Bi). Physical Review A: Atomic, molecular, and optical physics [1990-2015], 2012, 86 (1), pp.012503 1-16. 10.1103/PhysRevA.86.012503. hal-00903408

HAL Id: hal-00903408

https://hal.science/hal-00903408

Submitted on 25 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Excitation energies from relativistic coupled-cluster theory of general excitation rank: Initial implementation and application to the silicon atom and to the molecules XH (X = As, Sb, Bi)

Mickaël Hubert

Laboratoire de Chimie et Physique Quantiques, IRSAMC, Université Paul Sabatier Toulouse III, 118 Route de Narbonne,
F-31062 Toulouse, France

Lasse K. Sørensen and Jeppe Olsen

Theoretical Chemistry, Langelandsgade 140, Aarhus University DK-8000 Århus C, Denmark

Timo Fleig

Laboratoire de Chimie et Physique Quantiques, IRSAMC, Université Paul Sabatier Toulouse III, 118 Route de Narbonne,
F-31062 Toulouse, France

(Received 3 February 2012; published 13 July 2012)

We present an implementation of four-component relativistic coupled-cluster theory for the treatment of electronically excited states of molecules containing heavy elements, allowing for a consistent and accurate treatment of relativistic effects such as the spin-orbit interaction and electron correlations as well as their intertwining. Our approach uses general excitation ranks in the cluster operator and, moreover, allows for the definition of active-space selected excitations of variable excitation rank. Initial applications concern the silicon atom and the heavier pnictogen monohydride molecules, where we focus on the first vertical excitation energy to the Ω = 1 electronic state. We discuss the problem of adequately choosing a reference state (Fermi vacuum) and addressing electron correlation in the presence of effects of special relativity of increasing importance. For the heaviest homolog, BiH, where dynamic electron correlation is of major importance, we obtain vertical excitation energies with a deviation of less than 1% from the experimental value.

DOI: 10.1103/PhysRevA.86.012503

PACS number(s): 31.15.bw, 31.15.am, 31.15.aj

## I. INTRODUCTION

Electronically excited states of small molecules containing heavy atoms play an important role in many research areas of modern physics. In the (ultra)cold molecular sciences [1] there is an increasing interest in experimentally generating molecules in their electronic and rovibrational ground state by photoassociation via an electronically excited state [2]. In astrophysics of stars [3], the understanding of collision processes in stellar atmospheres [4] involves the knowledge of molecular excited states, including both main group and transition-metal atoms. As an example from fundamental physics, various extensions to the standard model of elementary-particle physics postulate electric dipole moments (EDM) of leptons [5]. Modern experiments search for the electron EDM in an electronically excited state of diatomic molecules and molecular ions containing a heavy atom [6]. The accurate determination of the electronic structure in excited states of the relevant molecules is of crucial importance in all of these and other research fields.

At present, the most accurate electronic-structure approach to the calculation of electronically excited states in atoms and molecules is the coupled-cluster (CC) method. Recent progress, including developments for excited states [7], has been documented in a monograph [8] covering this highly active field of many-body theory. When turning to the treatment of heavy elements where relativistic generalizations of these methods are required, the general challenge of implementing such methodology becomes manifest in their scarcity (see [9], and references therein). To date, the only relativistic CC methods for the treatment of molecular excited states are the intermediate Hamiltonian Fock-space CC method (IH FSCC) [10,11] by Visscher, Eliav, and co-workers and higher-order

correlation methods [12] by Hirata and co-workers using the equation-of-motion (EOM) CC formalism [13,14]. IH FSCC is limited in that it is not generally applicable and the treatment of excitation ranks higher than doubles in the wave operator is currently not possible. The method of Hirata et al. is restricted to the use of two-component valence pseudospinors based on a relativistic effective core potential (RECP) including spin-orbit interaction [15]. Such an approach lacks both the rigor and the flexibility of all-electron four-component methods which use a frozen-core approximation for the electrons of atomic cores.

Our developments aim at a rigorous assessment of the electronically excited states of small molecules including heavy elements, a general challenge in the relativistic electronic many-body problem until today [9]. Central elements of our methodology are (1) a rigorous treatment of special relativity using four-component all-electron Dirac Hamiltonians at all stages of the calculation; (2) methods of general excitation rank in the wave operator; and (3) methods based on developments of the wave function in a basis of strings of particle creation operators in second quantization, so-called string-based methods [16].

In this paper we present a relativistic coupled-cluster implementation based on linear-response theory and four-component relativistic Hamiltonian operators for the calculation of molecular excited states. In the following section on general theory (Sec. II) we review the description of electronically excited states in CC theory (Sec. II A) and our previous relativistic CC approach for electronic ground states (Sec. II B). Section III describes our implementation, in particular, the algorithm for calculating the relativistic CC Jacobian matrix. Here we also present an analysis of the computational scaling of our approach. Section IV is
concerned with initial applications of the method. We have chosen an atomic case (Si) featuring excited states of two different kinds: excited states due to the first-order spin-orbit splitting within a spectroscopic term, and states corresponding to a different spectroscopic term. As a second and molecular example we apply our approach to the second-order spin-orbit splitting of the $^3\Sigma$ ground state of heavier pnictogen hydrides, a notoriously difficult problem [17,18] requiring the treatment of static and dynamic electron correlation as well as spin-dependent magnetic interactions accurately. In the final section (Sec. V) we summarize and draw conclusions from our findings.

## II. THEORY

### A. Excited states in coupled-cluster theory

Response theory comprises a general and powerful framework for the calculation of atomic and molecular properties [19] as well as excitation energies based, e.g., on CC wave functions [20]. Here, the simple poles of the linear response function correspond to the excitation energies and occur at the eigenvalues of the CC Jacobian matrix.

An alternative way of deriving the CC Jacobian proceeds by an analogy to configuration interaction (CI) theory. Using CC language the CI Schrödinger equation with subtracted ground-state energy $E_0$ can be rewritten as

$$ (\hat{H} - E_0)(t_0 \hat{1} + \hat{T})|\Phi\rangle = 0, \quad (1) $$

where $|\Phi\rangle$ is the reference (or Fermi vacuum) state, $\hat{T} = \sum_{\mu} t_{\mu} \hat{\tau}_{\mu}$ is the cluster excitation operator with $\hat{\tau}_{\mu} \in \{\hat{\tau}_{i}^{a}, \hat{\tau}_{ij}^{ab}, \dots\}$, $\hat{\tau}_{i}^{a} = \hat{a}_{i}^{\dagger} \hat{a}_{i}$ a single-replacement operator in equal-time second-quantization representation, and $t_{\mu}$ the corresponding expansion coefficient.

Projection with the CI excitation manifold $\langle\psi_{\mu}| = \langle\Phi|\hat{\tau}_{\mu}^{\dagger}$ onto Eq. (1) yields a set of CI coefficient equations for the CI vector function

$$ \Omega_{\mu}^{\text{CI}} = \langle \psi_{\mu} | (\hat{H} - E_0)(t_0 \hat{1} + \hat{T}) | \Phi \rangle = 0. \quad (2) $$

Taking the derivative with respect to all expansion parameters defines an Hermitian CI Jacobian, the matrix elements of which become

$$ \begin{aligned} A_{\mu\nu}^{\text{CI}} &= \frac{\partial}{\partial t_{\nu}} \Omega_{\mu}^{\text{CI}} = \langle \psi_{\mu} | (\hat{H} - E_0) \hat{\tau}_{\nu} | \Phi \rangle \\ &= \langle \psi_{\mu} | \hat{H} | \psi_{\nu} \rangle - E_0 \delta_{\mu\nu}. \end{aligned} \quad (3) $$

Obviously, diagonalization of the matrix $A^{\text{CI}}$ yields excitation energies from CI theory. It is straightforward to construct the analogy in CC theory. Here, the amplitude equations corresponding to Eq. (2) are cast (in linked form) as

$$ \Omega_{\mu}^{\text{CC}} = \langle \psi_{\mu} | e^{-\hat{T}} \hat{H} e^{\hat{T}} | \Phi \rangle = 0, \quad (4) $$

with the same excitation manifold $\langle\Phi|\hat{\tau}_{\mu}^{\dagger}$, and the derivative matrix is obtained as

$$ A_{\mu\nu}^{\text{CC}} = \frac{\partial}{\partial t_{\nu}} \Omega_{\mu}^{\text{CC}} = \langle \psi_{\mu} | e^{-\hat{T}} [\hat{H}, \hat{\tau}_{\nu}] e^{\hat{T}} | \Phi \rangle. \quad (5) $$

Consequently, diagonalization of the matrix $A^{\text{CC}}$ yields excitation energies from CC theory. The difference between Eqs. (3) and (5) can be reduced to the difference in parametrization of

the wave function, linear in CI theory and exponential in CC theory, respectively. The excitation energies $\omega_A$ obtained from the eigenvalue equations

$$ A^{\text{CC}} |\psi_f\rangle = \omega_{Af} |\psi_f\rangle \quad (6) $$

are equivalent to those from the EOM CC theory [13,21].

For reasons of computational efficiency, Eq. (6) is solved iteratively by algorithms similar to direct CI techniques, but in the present case for a non-Hermitian matrix $A^{\text{CC}}$. It has been shown earlier [22] how such linear transformations with the CC Jacobian can be evaluated for CC theory with general excitation levels of the cluster operator. This becomes possible by performing subsequent CI expansions using a general CI program.

### B. Four-component relativistic approach

We have in the present work generalized the nonrelativistic implementation of Ref. [22] to a relativistic formalism where four-component or two-component relativistic Hamiltonian operators may be used from the outset, and our implementation will be described in Sec. III. Our approach to treating special relativity is identical to the one presented in Refs. [23–25]. In summary, the cluster operators $\hat{T} = \sum_{m} \hat{T}_{m}$ are generalized to include the possibility of flipping the Kramers projection of the underlying spinors along with the excitation, e.g., for singles replacements:

$$ \hat{T}_1 = \sum_{ia} \{t_i^a \hat{\tau}_i^a + t_i^a \hat{\tau}_i^a + t_i^{\bar{a}} \hat{\tau}_i^{\bar{a}} + t_i^{\bar{a}} \hat{\tau}_i^{\bar{a}}\}. \quad (7) $$

The same generalization of excitation operators also applies to the operators $\hat{\tau}_v$ in the CC Jacobian matrix, Eq. (5). The approach is therefore Kramers restricted, in the sense that the underlying four-component spinors $\{\phi_i, \phi_{\bar{i}}\}$ form time-reversal partners (Kramers pairs)

$$ \hat{K} \phi_i = \phi_{\bar{i}}, \quad \hat{K} \phi_{\bar{i}} = -\phi_i, \quad (8) $$

and that this symmetry is exploited for reducing the number of unique Hamiltonian one- and two-particle integrals [16]. An arbitrary number of spinor spaces with arbitrary occupation restraints may be used [generalized active spaces (GAS)] [25,26], which allows for the description of the multireference character of electronic states via active-space selected higher excitations. Double point group symmetry has been implemented for the real-valued [27,28] matrix groups $D_{2h}$, $D_{2v}$, and $C_{2v}$. This ensures for these cases a completely real-valued formalism, also when spin-orbit interaction is included. Our implementation is interfaced to a local version of the DIRAC relativistic electronic-structure package [29]. Currently, this local version limits the present method to the use of the four-component Dirac-Coulomb Hamiltonian (in Born-Oppenheimer approximation)

$$ \begin{aligned} \hat{H}^{\text{DC}} = & \sum_{A} \sum_{i} [c(\vec{\alpha} \cdot \vec{p})_i + \beta_i m_0 c^2 + V_{iA}] + \sum_{i,j>i} \frac{1}{r_{ij}} \mathbb{I}_4 \\ & + \sum_{A,B>A} V_{AB}, \end{aligned} \quad (9) $$

where $V_{iA}$ is the potential-energy operator for electron $i$ in the electric field of nucleus $A$, and $V_{AB}$ represents the potential
energy due to the internuclear electrostatic repulsion of the clamped nuclei.

### III. IMPLEMENTATION: ALGORITHM FOR THE RELATIVISTIC COUPLED-CLUSTER JACOBIAN

We now proceed to an outline of our implementation of the eigenvalue equation (6). It may be regarded as a combination of the algorithms described in Refs. [22,23]. Based on the techniques developed for general relativistic CI expansions [30] we evaluate the linear transformation of a coefficient trial vector x with the CC Jacobian

$$J_{\mu}^{\text{CC}} = \sum_{\nu} A_{\mu\nu}^{\text{CC}} x_{\nu} = \sum_{\nu} \langle \psi_{\mu} | e^{-\hat{T}} [\hat{H}, \hat{\tau}_{\nu}] e^{\hat{T}} | \Phi \rangle x_{\nu}. \quad (10)$$

Following Refs. [22,23] Eq. (10) is solved in four steps:

(1) $|a\rangle = e^{\hat{T}}|\Phi\rangle = (\sum_{n=0}^{\infty} \frac{1}{n!} \hat{T}^n)|\Phi\rangle$. The individual terms in the Taylor expansion comprise repeated transformations of the form $\hat{T}|\psi\rangle$. We here employ the modified relativistic GAS CI implementation of Refs. [23,30] and the final ground-state cluster amplitudes in $\hat{T}$. The Taylor expansion truncates naturally upon exhausting the possible excitations on a given reference vector.

(2) $|b\rangle = [\hat{H}, \hat{\tau}_v]|a\rangle = (\hat{H}\hat{\tau}_v - \hat{\tau}_v\hat{H})|a\rangle$. Here, $\hat{H}\hat{\tau}_v|a\rangle$ corresponds to the calculation of a sigma vector [30] from the reference vector $\hat{\tau}_v|a\rangle$. In the second term $-\hat{\tau}_v$ is applied to the sigma vector $\hat{H}|a\rangle$, and the resulting vectors from the two terms are added yielding the commutator.

(3) $|c\rangle = e^{-\hat{T}}|b\rangle = (\sum_{n=0}^{\infty} \frac{(-1)^n}{n!} \hat{T}^n)|b\rangle$. These transformations are evaluated in the same manner as those in step 1.

(4) $J_{\mu}^{\text{CC}} = \langle \psi_{\mu} | c \rangle = \langle \Phi | \hat{\tau}_{\mu}^{\dagger} | c \rangle$. This final step corresponds to the evaluation of a general transition density, which is also possible employing the modified relativistic GAS CI implementation in Refs. [23,30].

Therefore, since the underlying relativistic CI program [30] can treat general excitation levels, we are here immediately able to compute a relativistic CC Jacobian at general excitation rank, both with respect to the cluster operators and the excitation operators.

However, as has been discussed in Refs. [22,23], the present algorithm suffers from an increased operation count compared to conventional (and nonrelativistic) CC implementations for excited states [31]. The increased operation count of CI-based CC has been analyzed earlier [23,32] for ground-state calculations and amounts to a computational scaling of the method as $O^{n+2}V^{n+2}$, where O is the number of occupied orbitals, V is the number of virtual orbitals, and n is the highest excitation rank of the cluster operators. In order to elucidate the scaling of the present algorithm for excited-state calculations, we rewrite the right-hand side of Eq. (5) as

$$A_{\mu\nu}^{\text{CC}} = \langle \psi_{\mu} | e^{-\hat{T}} \hat{H} \hat{\tau}_{\nu} e^{\hat{T}} | \Phi \rangle - \langle \psi_{\mu} | e^{-\hat{T}} \hat{\tau}_{\nu} \hat{H} e^{\hat{T}} | \Phi \rangle. \quad (11)$$

Starting with the second term on the right-hand side of Eq. (11), we reexpress the term

$$e^{-\hat{T}} \hat{\tau}_{\nu} = \hat{\tau}_{\nu} - \sum_{\mu} t_{\mu} \hat{\tau}_{\mu} \hat{\tau}_{\nu} + \frac{1}{2} \left( \sum_{\mu} t_{\mu} \hat{\tau}_{\mu} \right)^2 \hat{\tau}_{\nu} - \dots, \quad (12)$$

which is seen to be a pure deexcitation operator acting on the bra vector $\langle \psi_{\mu} |$. Therefore, the highest excitation rank n of the

excitation manifold $\langle \psi_{\mu} |$ is reduced to $n-m$, where m is the excitation rank of an individual term in Eq. (12). Since $m \ge 1$, the highest excitation rank in $\langle \psi_{\mu} | e^{-\hat{T}} \hat{\tau}_{\nu} | \Phi \rangle$ is $n-1$. This means that in order for $\hat{H}e^{\hat{T}}|\Phi\rangle$ to be connected to this modified excitation manifold, $e^{\hat{T}}|\Phi\rangle$ has to contain excitations up to rank $n+1$ since the Hamiltonian has a maximum down rank of 2. The second term on the right-hand side of Eq. (11) therefore exhibits a computational scaling of $O^{n+1}V^{n+2}$, since in general the highest excitation rank k present in the excitation manifold entails a scaling with the number of occupied orbitals as $O^{k+2}$.

In contrast to this, the first term on the right-hand side of Eq. (11) has no additional cluster operator to the left of the Hamiltonian. This means that for $\hat{H}\hat{\tau}_{\nu}e^{\hat{T}}|\Phi\rangle$ to be connected to the original excitation manifold with rank n, $\hat{\tau}_{\nu}e^{\hat{T}}|\Phi\rangle$ has to contain excitations up to rank $n+2$. Thus, this term is the highest-scaling term of the algorithm and the total algorithm for the CI-based Jacobian scales as $O^{n+2}V^{n+2}$, exactly as does the CI-based algorithm for the ground-state vector function [23,32]. In typical applications, the dimension of the extended space defined by $\hat{H}\hat{\tau}_{\nu}e^{\hat{T}}|\Phi\rangle$ is one to three orders of magnitude larger than the dimension of the excitation manifold. The fact that the cluster operators $\hat{T}$, $\hat{\tau}_{\nu}$, and the Hamiltonian are now relativistic operators has no bearing for the present discussion concerning the orders for the computational scaling. There are, however, increased scaling prefactors in a relativistic algorithm which has been analyzed for commutator-based CC in Ref. [25].

### IV. APPLICATION AND ANALYSIS

In this section we present applications to the silicon atom and to heavier pnictogen hydrides. The silicon atom has been chosen as an initial test case to verify the applicability of our method. We focus on all Russell-Saunders terms originating from the atomic configuration $3p^2$, i.e., $^3P_{2,1,0}$, $^1D_2$, and $^1S_0$. Therefore, the problem comprises excited states arising from the same term due to first-order spin-orbit splitting ($^3P_2$, $^3P_1$, $^3P_0$) and excited states from different terms corresponding to the same electronic configuration. We define the Fermi vacuum determinant as the one where the energetically lowest Kramers pairs are all doubly occupied, i.e., the configuration $1s^22s^22p^63s^23p^2_{1/2}$. The purpose here is to show in a simple way the coupled-cluster description of several different excited states by taking into account spin-orbit interaction.

Turning to the molecules, the pnictogen hydrides are characterized by two valence electrons occupying the $(\pi_{1/2}, \pi_{-1/2})$ and $(\pi_{3/2}, \pi_{-3/2})$ Kramers pairs which are here denoted as the spin-orbit split $\pi$ orbitals assigning $\lambda_{\omega}$ quantum numbers. $\lambda$ is an approximate quantum number as spin-orbit interaction mixes orbitals of different angular momentum projection $m_{\ell}$, e.g., $\sigma$ character into the $\pi$ orbitals, $\sigma_{1/2} - \pi_{1/2}$. Their occupation and character therefore differ depending on the pnictogen atom. Since we describe the systems in a spinor basis $\{\varphi_{\omega}\}$ our natural choice of Fermi vacuum is the closed-shell valence occupation $\pi_{1/2}^1, \pi_{-1/2}^1$. Such a reference state is a good approximation to the wave function in the case of BiH where spin-orbit interaction is strong. However, the multireference (MR) character is expected to become more
and more important toward the lighter homologs where the $\pi_{1/2}$ and $\pi_{3/2}$ spinors become quasidegenerate. The first excited state, with $\Omega = 1$, is predominantly described by $\pi_{3/2}^{1}, \pi_{-1/2}^{1}$ and the excited state with $\Omega = -1$ is predominantly $\pi_{1/2}^{1}, \pi_{-3/2}^{1}$.

Dynamic correlation and relativistic effects are treated on the same footing. The multireference character of states is taken into account via active-space selected higher excitations. Since the GAS-CC method is not strictly invariant to the choice of Fermi vacuum state, we expect a bias of the CC wave function depending on the chosen reference determinant.

### A. Computational details

All calculations were performed with the DIRAC relativistic electronic-structure program package, using the latest version [29] for the Hartree-Fock calculations and integral transformations, and a local development version for the CC calculations.

For our initial test calculations on the silicon atom we have tested different basis sets and resorted to using the atomic-natural-orbital (ANO) Relativistic and Core-Correlating (RCC) basis set [33] and to include only the four valence electrons in the correlation treatment. We employed Dyall's triple-$\zeta$ and quadruple-$\zeta$ basis sets in uncontracted form [34,35] for Bi, Sb, and As. The listed valence and core-correlating functions for the Bi 5d, 6s, and 6p shells, for the Sb 4d, 5s, and 5p shells, and for the As 3d, 4s, and 4p shells have all been included. For H we used Dunning's cc-pVTZ-DK and cc-pVQZ-DK basis sets in uncontracted form [36]. The internuclear distances for AsH, SbH, and BiH are the experimental ones [37,38].

We employed the four-component Dirac-Coulomb Hamiltonian, Eq. (9), throughout. Thus, our models describe one- and two-electron spin-own-orbit coupling and spin-own-orbit-correlation coupling rigorously. We currently do not include spin-spin coupling and spin-other-orbit interactions due to limitations in the implemented Hamiltonian operators. Kramers-paired spinors for the subsequent GAS-CC calculations were obtained from all-electron closed-shell Dirac-Coulomb-Hartree-Fock calculations. In addition, we performed for comparison exemplifying CI and CC calculations based on open-shell average-of-configuration Dirac-Coulomb-Hartree-Fock (DCHF) wave functions. In the open-shell DCHF calculations fractional occupation numbers are introduced in the Fock operator using minimal spaces of Kramers pairs (two electrons in the three 3p Kramers pairs in the case of Si, and two electrons in the two $\pi$ valence Kramers pairs in the molecular cases).

Electron correlations are described in various fashions. On the one hand, we apply the standard CC hierarchy [39] up to full iterative quintuple excitations [CCSDTQP; S = single excitations with respect to the reference state $|\Phi\rangle$, D = double excitations, T = triple, Q = quadruple, and P = pentuple (fivefold) excitations] in the case of SbH. The construction of active spaces is done in an efficient manner by exploiting the GAS concept [26]. In the present case, important subsets of the possible model spaces are denoted as CC($n_m$) models [40] (see Fig. 1).

The silicon atom is treated with the standard CC series CCSD to CCSDTQ (the latter of which in this case corresponds to full valence CC), a series of models excluding the correlation

<figure><img src="image_2.png" /><figcaption>FIG. 1. CC($n_m$) ($n > m$) as a subset of GAS excitation manifolds. $n$ is the maximum number of holes in the occupied subspace; $m$ is the maximum number of particles in the virtual subspace.</figcaption></figure>

of the 3s electrons among each other (S2CC), and different CC($n_m$) schemes. The detailed definitions of these models are to be found in the Appendix.

For all molecular systems we choose to correlate the six valence electrons. The correlation of 3d electrons of As, 4d electrons of Sb, and 5d electrons of Bi does play a role in assessing the ground-state spin-orbit splitting, as studied for the case of BiH by Knecht et al. [41], but the effect is only on the order of $+75 \text{ cm}^{-1}$ for this latter molecule. Also here, we use various electron correlation models, the details for the specification of which are to be found in the Appendix.

Comparative four-component generalized-active-space configuration interaction (GAS-CI) calculations were performed with the KR-CI module [18,42] of the DIRAC program package [29]. This approach makes use of the GAS concept in the same way as our presented GAS-CC method [23]. Closed-shell CI calculations have been performed using the newly-implemented linear symmetry double groups [43,44].

### B. Results and discussion

We present and discuss in this section our results for Si, AsH, SbH, and BiH. The bulk of the GAS-CC and GAS-CI calculations was carried out with a virtual spinor space size for which the total energy has been converged (see Sec. 3 of the Appendix). The principal purpose here is to show for the case of a few representative model systems the performance of standard CC and CC($n_m$) models, and to compare these with a genuine-but linearly parametrized—MR approach, MRCI.

#### 1. Silicon atom

First of all, we performed benchmark CI and MRCI calculations (which are subsets of GAS-CI; see the Appendix) to guide our CC study and to assess the leading effects on excited-state energies. These calculations are presented in Table I. Whereas the os-MRCI calculations consistently yield results close to the FCI values, the truncated closed-shell approaches fail in describing the excited states correctly. Given a balanced starting point for the different electronic states by using average-of-configuration DCHF these are already well
<table><caption>TABLE I. Excitation energies *T* in cm⁻¹ for the ³P₁, ³P₂, ¹D₂, and ¹S₀ excited states of the Si atom, with different relativistic CC models (defined in Tables V and VI), the Dirac-Coulomb Hamiltonian, and a closed-shell Dirac-Hartree-Fock reference state. MRCI calculations based on an open-shell multireference state (os) and a closed-shell single-reference state (cs). We used a complete active space of two electrons in the three 3p orbitals. MRCI models are defined in Table VII. All calculations were performed including single and double excitations of the 3s electrons, except where marked otherwise (S2, only singles from the 3s shell). The basis sets are of ANO-RCC quality; the cutoff for the virtual spinors is set to 10 a.u. (see text).</caption><thead><tr><th colspan="10">Method</th></tr><tr><th colspan="4">os-MR</th><th colspan="4">cs</th></tr><tr><th>State</th><th>CISD</th><th>CISDT</th><th>FCI</th><th>CISD</th><th>CISDT</th><th>FCI</th><th>Expt. [52]</th></tr></thead><tbody><tr><td><sup>3</sup>P<sub>1</sub></td><td>79</td><td>79</td><td>79</td><td>5275</td><td>1768</td><td>89</td><td>77.1</td></tr><tr><td><sup>3</sup>P<sub>2</sub></td><td>229</td><td>228</td><td>228</td><td>7368</td><td>1519</td><td>258</td><td>223.2</td></tr><tr><td><sup>1</sup>D<sub>2</sub></td><td>6475</td><td>6441</td><td>6413</td><td>13447</td><td>7595</td><td>6435</td><td>6299.8</td></tr><tr><td><sup>1</sup>S<sub>0</sub></td><td>15606</td><td>15691</td><td>15551</td><td>20917</td><td>18821</td><td>15574</td><td>15394.4</td></tr></tbody><thead><tr><th colspan="10">Method (all cs)</th></tr><tr><th>State</th><th>S2CCSD</th><th>S2CC(3<sub>2</sub>)</th><th>S2CCSDT</th><th>CCSD</th><th>CC(4<sub>2</sub>)</th><th>CCSDT</th><th>CC(4<sub>3</sub>)</th><th>FCC<sup>a</sup></th><th>Expt. [52]</th></tr></thead><tbody><tr><td><sup>3</sup>P<sub>1</sub></td><td>-1225</td><td>-86</td><td>180</td><td>-784</td><td>-128</td><td>-438</td><td>81</td><td>89</td><td>77.1</td></tr><tr><td><sup>3</sup>P<sub>2</sub></td><td>1318</td><td>687</td><td>277</td><td>1750</td><td>666</td><td>-222</td><td>276</td><td>258</td><td>223.2</td></tr><tr><td><sup>1</sup>D<sub>2</sub></td><td>10422</td><td>6648</td><td>6610</td><td>10726</td><td>6504</td><td>5942</td><td>6446</td><td>6435</td><td>6298.8</td></tr><tr><td><sup>1</sup>S<sub>0</sub></td><td>22392</td><td>20108</td><td>19745</td><td>20121</td><td>15926</td><td>17188</td><td>15591</td><td>15574</td><td>15394.4</td></tr></tbody></table>

<sup>a</sup>Values obtained by an equivalent full CI calculation due to an unresolved instability of our FCC calculation at a 10 a.u. cutoff value. We have obtained identical FCC and FCI values for smaller dimensions of the virtual space confirming the proper functionality of the new CC code.

described by MR-CISD, and higher excitations hardly play a role for energy differences.

In the case of a cs reference state we observe that the results at low excitation levels are largely off the mark, and higher excitations gradually lead to improvements, with only the Full (F) CI and Full (F) CC and CC(4<sub>3</sub>) models yielding accurate results. We rationalize and explain this behavior first by analyzing the Fermi-vacuum determinant of our reference state. The problem is simplified by considering only the two valence *p* electrons. Since we use four-component spinors throughout, we write the determinant in terms of good quantum numbers for the two particles |j(i), m<sub>j</sub>(i)⟩ and further express this determinant in terms of Russell-Saunders coupled states <ML<sub>J</sub>, giving

$$ \left| \frac{1}{2}, \frac{1}{2}; \frac{1}{2}, -\frac{1}{2} \right| = -\sqrt{\frac{2}{3}} {}^3P_0 - \frac{1}{\sqrt{3}} {}^1S_0, \quad (13) $$

the details of which are to be found in the Appendix, Sec. 4. The reference determinant is therefore biased toward the *J* = 0 state of the ³P₀ term and contains a significant admixture from the ¹S₀ term. In order to interpret the excitation energies we expand the singly excited determinants in the same manner, e.g.,

$$ \begin{aligned} \left| \frac{3}{2}, \frac{1}{2}; \frac{1}{2}, \frac{1}{2} \right| &= -\frac{1}{2} {}^3P_1 - \frac{1}{2} {}^3P_2 + \frac{1}{\sqrt{2}} {}^1D_2, \\ \left| \frac{3}{2}, \frac{3}{2}; \frac{1}{2}, \frac{1}{2} \right| &= \sqrt{\frac{2}{3}} {}^1D_2 + \frac{1}{\sqrt{3}} {}^3P_2. \end{aligned} \quad (14) $$

Based on a closed-shell model we ensuing expect a CISD calculation to strongly overestimate the excitation energies of the ³P₁, ³P₂, and ¹D₂ states, since a single excitation is required for their description, leaving them uncorrelated, in contrast to the ground state. The data clearly confirms this. Furthermore, the overestimation for the ¹S₀ state is smaller than, e.g., the

one for the ¹D₀ state, since the former is partially represented in the reference state.

Continuing the argument, a CISDT calculation introduces triple excitations, in addition to the already present single and double excitations. Now, Eqs. (14) show that, for example, the ¹D₂ state is largely represented by singly excited determinants relative to our reference determinant. Since some of the triple excitations in the CISDT model are double excitations combined with single excitations required to qualitatively describe the ¹D₂ excited state, dynamic electron correlation effects are taken into account for the ¹D₂ state, in contrast to the CISD model. We therefore expect the ¹D₂ excitation energy to be much closer to the experimental value in the CISDT model, which the data confirms.

A similar reasoning applies to explain the obtained results for the ¹S₀ state and the ³P₁ and ³P₂ components of the ground-state term. For example, the CISDT model shows the smallest correction for the ¹S₀ state which is already partially correlated in the CISD model. A final, but smaller, correction is obtained by adding quadruple excitations, this correction now being largest for the ¹S₀ state which appears in doubly excited determinants

$$ \begin{aligned} \left| \frac{3}{2}, \frac{3}{2}; \frac{3}{2}, -\frac{3}{2} \right| &= \frac{1}{\sqrt{3}} {}^3P_2 + \frac{1}{\sqrt{6}} {}^1D_2 - \frac{1}{\sqrt{6}} {}^3P_0 + \frac{1}{\sqrt{3}} {}^1S_0, \\ \left| \frac{3}{2}, \frac{1}{2}; \frac{3}{2}, -\frac{1}{2} \right| &= \frac{1}{\sqrt{3}} {}^3P_2 + \frac{1}{\sqrt{6}} {}^1D_2 + \frac{1}{\sqrt{6}} {}^3P_0 - \frac{1}{\sqrt{3}} {}^1S_0, \end{aligned} \quad (15) $$

relative to the reference determinant.

Turning to the CC results in the light of these findings, the CCSD model, containing higher excitations than doubles in disconnected terms (which contribute to the CC energy indirectly due to the coupling of excitations in the CC amplitude equations), yields results of accuracy between CISD and CISDT for the ¹D₀ and ¹S₀ states, but the correction
overshoots for the $^3P_1$ and $^3P_2$ components of the ground-state term. CC(4$_2$) includes the important higher excitations to give qualitatively correct results except for the first-order spin-orbit splitting, where there is a residual error of several hundred cm$^{-1}$. Notably, CCSDT does not improve upon the CC(4$_2$) results. CC(4$_3$) brings about a significant correction for the two first excited states making all of them qualitatively correct, although the excitation level in the active space is the same as in CC(4$_2$). FCC (CCSDTQ) yields only a minor correction to this last model.

The slightly different values obtained with closed-shell and open-shell FCC models reflect the differently polarized core and virtual spinors of the atom depending on the DCHF model used.

We conclude that due to the specific choice of coupling picture, here $j-j$, a Fermi vacuum determinant represented in this coupling picture may not comprise a good description of the electronic ground state. As a consequence, CC models based on this vacuum state and truncated at low excitation ranks may yield large errors in calculated excitation energies. Active-space selected higher excitations largely correct for the ensuing errors in atomic excitation energies. High accuracy in the spin-orbit splitting of the Si atomic ground state (excited states $^3P_{1,2}$) is only achieved if in addition dynamic electron correlations are accounted for through at least triple excitations into the virtual spinor space. CCSDT alone, however, is not accurate enough in describing these states which according to Table IV in Ref. [22] means that for the given case of a poor Fermi vacuum state the excitation energy has to be described at least through fifth order in the fluctuation potential. This is approximately the case for the CC(4$_3$) model which includes important quadruple excitations.

## 2. The individual pnictogen hydrides

(a) Arsenic monohydride—AsH. This system exhibits a small spin-orbit splitting of the $\pi$ Kramers pairs and, therefore, the $\Omega = 0$ ground state is likely to have strong MR character. Indeed, closed-shell CCSD gives a qualitatively wrong energy estimation and even a state inversion which can been seen in Table II. We rationalize this failure of CCSD by again closely analyzing the Fermi vacuum determinant in the molecular case (see the Appendix, Sec. 5). Based on the analysis, our reference state can be qualitatively described as

$$|(m_j)_1; (m_j)_2|_{\Omega} = \left|\left(\frac{1}{2}\right); \left(-\frac{1}{2}\right)\right|_0 \approx c_3|^3\Sigma_0\rangle + c_1|^1\Sigma_0\rangle, \quad (16)$$

where $c_3 \approx -c_1 = \frac{1}{\sqrt{2}}$ in the case of AsH. This in turn means that the true ground state $|^3\Sigma_0\rangle$ is best represented by a linear combination of the determinants $|(\frac{1}{2})\rangle; (-\frac{1}{2})|_0$ and $|(\frac{3}{2})\rangle; (-\frac{3}{2})|_0$. It therefore requires higher CC excitations in our single-reference approach to describe the presence of such strongly contributing determinants in the ground state.

A true MR method such as truncated MRCI, as we apply it here for comparison, gives qualitatively correct results but, being a method with a linear wave function parametrization, is limited by its intrinsic properties. CC(4$_2$) with open-shell spinors describes the MR character qualitatively, but still produces an error of about 80 cm$^{-1}$ compared with experiment.

<table><caption>TABLE II. Vertical excitation energies ($T_v$) in cm$^{-1}$ for the $\Omega = 1$ state of the AsH molecule at the experimental bond length of $R_e^{(0+)} \approx R_e^{(1)} \approx 1.5349$ Å [37] with different relativistic CC models (defined in Tables VIII and IX), the Dirac-Coulomb Hamiltonian, and a closed-shell Dirac-Hartree-Fock reference state. “os” refers to an average-of-configuration reference state with an averaging for two electrons in the two $\pi$ orbitals. MRCI calculations are based on an open-shell multireference state except if indicated otherwise (cs). We used a formal core space of four electrons in the $4s,\sigma$ orbitals and a complete active space of two electrons in the two $\pi$ orbitals. The model SD\_CISD corresponds to $a = 2, b = 4$ in Table X. The basis sets are of TZ quality (see text), except where marked otherwise. Converged cutoff for virtual spinors is 10 a.u. except if indicated otherwise.</caption><thead><tr><th>Method</th><th>Δ<sub>SO</sub> (cm<sup>-1</sup>)</th><th>No. CC amplitudes/CI det.</th></tr></thead><tbody><tr><td>csCISD</td><td>11231</td><td>36.015</td></tr><tr><td>csCISDT</td><td>-2899</td><td>1.556.976</td></tr><tr><td>csCISDTQ</td><td>170</td><td>28.650.840</td></tr><tr><td>SD_CISD</td><td>104</td><td>123.472</td></tr><tr><td>SD_CISDT</td><td>103</td><td>2.123.792</td></tr><tr><td>SDT_CISDT</td><td>102</td><td>3.875.024</td></tr><tr><td>SDTQ_CISDTQ</td><td>106</td><td>50.672.784</td></tr><tr><td>CCSD</td><td>-2661</td><td>36.016</td></tr><tr><td>CCSD-10000 a.u.</td><td>-2666</td><td>162.240</td></tr><tr><td>CCSD-26 a.u. (QZ)</td><td>-2707</td><td>147.015</td></tr><tr><td>CCSD-59 a.u. (QZ)</td><td>-2707</td><td>191.535</td></tr><tr><td>CC(4<sub>2</sub>)</td><td>-139</td><td>123.471</td></tr><tr><td>osCC(4<sub>2</sub>)</td><td>39</td><td>123.471</td></tr><tr><td>CC(4<sub>2</sub>)-26 a.u. (QZ)</td><td>-258</td><td>511.771</td></tr><tr><td>osCC(4<sub>2</sub>)-26 a.u. (QZ)</td><td>39</td><td>511.771</td></tr><tr><td>CCSDT</td><td>-71</td><td>1.556.975</td></tr><tr><td>osCCSDT</td><td>-79</td><td>1.556.975</td></tr><tr><td>CC(4<sub>3</sub>)</td><td>68</td><td>2.123.791</td></tr><tr><td>osCC(4<sub>3</sub>)</td><td>71</td><td>2.123.791</td></tr><tr><td>CCSDTQ</td><td>116</td><td>28.650.839</td></tr><tr><td>osCCSDTQ</td><td>107</td><td>28.650.839</td></tr><tr><td>Expt. [37] (T<sub>e</sub> = 2Λ<sub>0</sub>)</td><td>117.7</td><td></td></tr></tbody></table>

(We deduce $T_e$ from the spin-splitting constant $\Lambda_0$ given in Ref. [37]: $T_e = 2\Lambda_0$). CCSDT is still not qualitatively correct, but including higher internal excitations, CC(4$_3$), we observe a positive and qualitatively correct first AsH excitation energy. Comparing osCCSDT, osCC(4$_2$), and osCC(4$_3$) we conclude that the quadruple excitations including the double excitations $\pi_{1/2}^2 \rightarrow \pi_{3/2}^2$ combined with double excitations into the virtual spinors are essential to describe the relative energies of the $\Omega = 0$ and $\Omega = 1$ states. This is confirmed by examining the double-excitation cluster amplitude $T_{\pi_{3/2}\pi_{3/2}\pi_{1/2}\pi_{1/2}}^{\pi_{3/2}\pi_{3/2}}$, which is $\approx 0.45$ for csCCSD and $\approx 0.79$ for csCCSDTQ in the ground state. With higher excitation ranks, open-shell and closed-shell LRCC results become quite the same, as expected. Since the excited state $\Omega = 1$ is essentially obtained by a single excitation $\pi_{1/2}^1 \rightarrow \pi_{3/2}^1$ from our closed-shell Fermi vacuum state, the excitation energy is again correct to the 5th perturbation order for our higher correlated calculation CCSDTQ, according to Table IV in Ref. [22]. CCSDTQ describes both the MR character and dynamic electron correlation accurately. SDTQ-CISDTQ and osCCSDTQ results are almost identical,
as expected, the difference with csCCSDTQ being due to the difference in the molecular orbital basis (see also Sec. IV B3 c). Corresponding closed-shell CI calculations exhibit the same deterioration at lower excitation levels as in the case of the two lowest excited states of the Si atom.

Turning to errors from the employed Dirac-Coulomb Hamiltonian, we estimate the effect of the Gaunt term to be around $-4 \text{ cm}^{-1}$ and correlation effects from the As 3d atomic shell to be $+2 \text{ cm}^{-1}$, according to recent exact two-component MRCISD calculations [45]. Our results indicate that basis-set errors are very small which, however, could be nonzero at CC(4$_3$) or CCSDTQ levels. The experimental excitation energy is well reproduced in our best calculation (csCCSDTQ) with a deviation of less than 1.5%.

(b) *Stibylene*—SbH. SbH is an intermediate case between AsH and BiH in the sense that the spin-orbit splitting of

<table><caption>TABLE III. Vertical excitation energies ($T_v$) in cm$^{-1}$ for the $\Omega = 1$ state of the SbH molecule at the experimental bond length of $R_e^{(0+)} \approx R_e^{(1)} \approx 1.7226$ Å [37] with different relativistic CC models (defined in Tables VIII and IX), the Dirac-Coulomb Hamiltonian, and a closed-shell Dirac-Hartree-Fock reference state. “os” refers to an average-of-configuration reference state with an averaging for two electrons in the two $\pi$ orbitals. MRCI calculations are based on an open-shell multireference state except if indicated otherwise (cs). We used a formal core space of four electrons in the $5s, \sigma$ orbitals and a complete active space of two electrons in the two $\pi$ orbitals. The model SD\_CISD corresponds to $a = 2, b = 4$ in Table X. The basis sets are of TZ quality (see text), except where marked otherwise. Converged cutoff for virtual spinors is 4 a.u. except if indicated otherwise.</caption><thead><tr><th>Method</th><th>Δ<sub>SO</sub> (cm<sup>-1</sup>)</th><th>No. CC amplitudes/CI det.</th></tr></thead><tbody><tr><td>csCISD</td><td>11474</td><td>30.376</td></tr><tr><td>csCISDT</td><td>-1347</td><td>1.205.176</td></tr><tr><td>csCISDTQ</td><td>705</td><td>20.370.586</td></tr><tr><td>SD_CISD</td><td>577</td><td>103.856</td></tr><tr><td>SD_CISDT</td><td>572</td><td>1.640.160</td></tr><tr><td>SDT_CISDT</td><td>563</td><td>2.987.792</td></tr><tr><td>SDTQ_CISDTQ</td><td>582</td><td>35.857.552</td></tr><tr><td>SDTQ_CISDTQP</td><td>582</td><td>192.560.560</td></tr><tr><td>CCSD</td><td>-1070</td><td>30.375</td></tr><tr><td>CCSD-100 a.u.</td><td>-1071</td><td>69.360</td></tr><tr><td>CCSD-6 a.u. (QZ)</td><td>-1073</td><td>82.140</td></tr><tr><td>CCSD-116 a.u. (QZ)</td><td>-1074</td><td>226.935</td></tr><tr><td>CC(4<sub>2</sub>)</td><td>476</td><td>103.855</td></tr><tr><td>osCC(4<sub>2</sub>)</td><td>555</td><td>103.855</td></tr><tr><td>CC(4<sub>2</sub>)-6 a.u. (QZ)</td><td>434</td><td>284.496</td></tr><tr><td>CCSDT</td><td>485</td><td>1.205.175</td></tr><tr><td>osCCSDT</td><td>436</td><td>1.205.175</td></tr><tr><td>CCSDT-6 a.u. (QZ)</td><td>482</td><td>5.376.100</td></tr><tr><td>CC(4<sub>3</sub>)</td><td>599</td><td>1.640.159</td></tr><tr><td>osCC(4<sub>3</sub>)</td><td>575</td><td>1.640.159</td></tr><tr><td>CC(4<sub>3</sub>)-6 a.u. (QZ)</td><td>627</td><td>7.397.616</td></tr><tr><td>CCSDTQ</td><td>641</td><td>20.370.585</td></tr><tr><td>osCCSDTQ</td><td>584</td><td>20.370.585</td></tr><tr><td>CCSDTQP</td><td>645</td><td>152.218.389</td></tr><tr><td>osCCSDTQP</td><td>612</td><td>152.218.389</td></tr><tr><td>Expt. (T<sub>e</sub>) [37]</td><td>654.97</td><td></td></tr></tbody></table>

the $\pi$ spinors becomes appreciable. However, based on our analysis in the Appendix, Sec. 5, we still expect the MR character of the ground state to be large, which is confirmed by the results. These results for stibylene are compiled in Table III. CCSD fails in much the same manner as for AsH, but the error has become smaller. We interpret this behavior by the coefficient $c_3$ now becoming larger than $c_1$ [in Eq. (16)] due to increased spin-orbit effects. Ensuingly, CC(4$_2$) gives a drastic amelioration with a qualitatively correct value and reveals an important contribution of the quadruple excitation including $\pi_{1/2} \rightarrow \pi_{3/2}$. In contrast to AsH, CCSDT no longer improves upon CC(4$_2$) in this case. A value of acceptable accuracy with an error of less than 10% is already achieved with CC(4$_3$), which was not the case in AsH. Nevertheless, quadruple excitations still play a significant role and largely correct for this residual error. Our most accurate value at the CCSDTQP level of 645 cm$^{-1}$ deviates by less than 2% from the experimental result (654.97 cm$^{-1}$ [37]). Comparing with the CI results in Table III the SDTQ-CISDTQP value is already leveled by CC(4$_3$), which is computationally significantly cheaper to perform. Thus, we here encounter a turning point where the advantage of the MRCI method of being a genuine MR approach is surpassed by the CC method due to its superior efficiency in treating higher excitations.

We now turn to residual errors from sources other than the correlation expansion. Considering errors from the truncated Hamiltonian operator, we estimate the effect of the Gaunt term to be $-13 \text{ cm}^{-1}$ and correlation contributions from the Sb atomic 4d shell to be $+5 \text{ cm}^{-1}$ according to recent exact-two-component (X2C)-MRCISD calculations [45]. From our most accurate CC model using different basis sets, CC(4$_3$), we infer a TZ-QZ basis-set error of $+28 \text{ cm}^{-1}$. Adding these estimated residual errors to our single most accurate result, CCSDTQP, the splitting amounts to 665 cm$^{-1}$, comprising a deviation of roughly 1.5% from the experimental value.

(c) *Bismuth monohydride*—BiH. The heaviest pnictogen homolog, BiH, is a quasi-single-reference case. The spin-orbit splitting of the $\pi$ spinors is significantly larger compared to AsH and SbH, as relativistic effects are sizable in this system. The weak MR character of the ground state renders the closed-shell Fermi vacuum a much better starting point in this case. The results for bismuth monohydride are compiled in Table IV. Despite the fact that BiH is a quasi-single-reference case in the present description, closed-shell CISD is insufficient, because a single excitation is required for describing the excited state, leaving it uncorrelated relative to the ground state. This interpretation is corroborated by the cs-CISDT model where the triple excitations correlating the excited state are added (doubles on top of singles). As expected, and in accord with the cs-CI results CCSD gives qualitatively good results and CC(4$_2$) brings about a large amelioration resulting in an error of less than 1%. It should be noted, however, that CCSD still displays an absolute error of nearly 700 cm$^{-1}$, indicating the remaining bias in our chosen Fermi vacuum state. Again, the importance of the quadruple excitation including $\pi_{1/2} \rightarrow \pi_{3/2}$ contribution is observed. However, a quite accurate value is obtained at this level, in contrast to SbH. Also here, CCSDT does not improve on CC(4$_2$). Our most accurate value at the CC(4$_3$) level is 4931 cm$^{-1}$ with a deviation of less than 0.3% from experiment (4917.1 cm$^{-1}$ [37]). In contrast to the CI
<table><caption>TABLE IV. Vertical excitation energies ($T_v$) in cm$^{-1}$ for the $\Omega = 1$ state of the BiH molecule at an internuclear distance of 1.80 Å (the experimental bond lengths are $R_e^0 = 1.805$ and $R_e^1 = 1.7912$ Å) [37] with different relativistic CC models (defined in Tables VIII and IX), the Dirac-Coulomb Hamiltonian, and a closed-shell Dirac-Hartree-Fock reference state. “os” refers to an average-of-configuration reference state with an averaging for two electrons in the two $\pi$ orbitals. MRCI calculations are based on an open-shell multireference state except if indicated otherwise (cs). We used a formal core space of four electrons in the $6s,\sigma$ orbitals and a complete active space of two electrons in the two $\pi$ orbitals. The model SD\_CISD corresponds to $a = 2, b = 4$ in Table X. The basis sets are of TZ quality (see text), except where marked otherwise. Converged cutoff for virtual spinors is 10 a.u. except if indicated otherwise.</caption><thead><tr><th>Method</th><th>Δ<sub>SO</sub> (cm<sup>-1</sup>)</th><th>No. CC amplitudes/CI det.</th></tr></thead><tbody><tr><td>csCISD</td><td>16710</td><td>61.441</td></tr><tr><td>csCISDT</td><td>4163</td><td>3.475.201</td></tr><tr><td>csCISDTQ</td><td>4854</td><td>83.488.225</td></tr><tr><td>SD_CISD</td><td>4595</td><td>212.137</td></tr><tr><td>SD_CISDT</td><td>4569</td><td>4.769.137</td></tr><tr><td>SDT_CISDT</td><td>4515</td><td>8.738.389</td></tr><tr><td>SDTQ_CISDTQ</td><td>4596</td><td>149.525.014</td></tr><tr><td>CCSD</td><td>4239</td><td>61.440</td></tr><tr><td>osCCSD</td><td>4103</td><td>61.440</td></tr><tr><td>CCSD-104 a.u.</td><td>4243</td><td>168.540</td></tr><tr><td>CCSD-11 a.u. (QZ)</td><td>4300</td><td>144.060</td></tr><tr><td>CCSD-86 a.u. (QZ)</td><td>4302</td><td>277.440</td></tr><tr><td>CC(4<sub>2</sub>)</td><td>4867</td><td>212.136</td></tr><tr><td>osCC(4<sub>2</sub>)</td><td>4763</td><td>212.136</td></tr><tr><td>CC(4<sub>2</sub>)-11 a.u. (QZ)</td><td>4892</td><td>501.408</td></tr><tr><td>CCSDT</td><td>4855</td><td>3.475.200</td></tr><tr><td>CC(4<sub>3</sub>)</td><td>4931</td><td>4.769.136</td></tr><tr><td>Expt. [37] (T<sub>c</sub>)</td><td>4917.1</td><td></td></tr></tbody></table>

results given in Table IV, the SDTQ-CISDTQ value does not reach the quality of the computationally significantly cheaper CC(4$_2$). A highly accurate description of both effects due to special relativity and electronic correlation is given here with the CC(4$_3$) model.

There are still some residual errors from different sources. A TZ-QZ basis-set correction of +25 cm$^{-1}$ is obtained by comparing CC(4$_2$) results. Considering errors from the truncated Hamiltonian operator, we estimate the effect of the Gaunt term to be -60 cm$^{-1}$ and correlation contributions from the Bi atomic 5d shell to be +75 cm$^{-1}$ according to exact-two-component (X2C) -MRCISD calculations performed by Knecht et al. [41]. The remaining deviation may be attributed to the fact that the potential curves for the 0$^+$ and 1 states are no longer parallel in the case of BiH, thus our vertical excitation energies slightly overshoot the experimental values for $T_c$. Therefore, residual errors are expected to largely compensate each other, which confirms the accuracy of our highest-level results.

### 3. Discussion of theoretical aspects across the series

In this section we draw a comparison between the three molecules focusing on selected theoretical issues.

(a) CC($n_m$) models and excitation rank of $\hat{T}$. In Fig. 2 we show the convergence evolution of the various CC models

for the vertical excitation energy of the $\Omega = 1$ state of the three molecules. CC(4$_2$) is set between CCSD and CCSDT and brings some essential quadruple excitations. It roughly equals CCSDT in quality, but at a lower cost. CC(4$_3$) is a good compromise between CCSDT and CCSDTQ; it provides a high accuracy and avoids full quadruple excitations into the virtual space. For BiH, high accuracy is reached using this model (deviation of 0.3%). These CC models depend on the excitation rank of the operator $\hat{T}$, which for the standard models CCSD to CCSDTQP is 2 to 5, respectively. For the CC($n_m$) model, the rank of $\hat{T}$ depends on the active-space structure. As we elucidate in Table IX, for GAS I (core spinors) we use a maximum rank of 2 to perform double core excitations toward the virtual spinors. In GAS II the maximum rank is 4 ($n = 4$), but those quadruple excitations are restricted to the $\pi_{3/2}$ spinors. Finally, in GAS III the rank is 2 or 3 ($m = 2$ or 3) in order to perform double or triple excitations toward virtual spinors. The CC($n_m$) approach enables a flexible adaptation for ground-state and excited-state calculations by taking into account important classes of excitation in the $\hat{T}$ operator for a system-tailored description.

(b) Multireference problem—Comparison of GAS-CC and MRCI. The multireference character on these three systems decreases toward the heaviest homolog BiH. This character is linked to the energy difference between the $\pi_{1/2}$ and the $\pi_{3/2}$ spinors. With GAS-CC, which is not a true MR CC approach, we have to impose a single-reference Fermi vacuum $\pi_{1/2}^1\pi_{-1/2}^1$. For lighter systems than BiH, we ensuinglly introduce a certain bias into CC wave function. However, we can compensate for this flawed point of departure with higher excitation ranks. In Fig. 3 we show a comparison between comparable CI and CC models in terms of deviation from experiment for the three molecules. In AsH where MR effects are strong, MR-CI remains superior to GAS-CC up to the level of triple excitations into the virtual spinor space. For SbH, multireference effects are still significant but we obtain a slightly better description at CC(4$_3$) level surpassing MR-CISDT (see Table III). The single-reference dominated system BiH is significantly better described with CC(4$_2$) already, improving on MR-CISD by more than 250 cm$^{-1}$ (see Table IV) with the same number of wave function parameters.

(c) Spinors from closed-shell or open-shell optimization. In closed-shell optimizations on a system with a near degeneracy of states, the energy gap between the occupied and the unoccupied valence spinors is largely overestimated. This gap becomes much more realistic in the open-shell models. Open-shell spinors could be used in cases of strong near degeneracy and where the excitation level must be kept low. For AsH, it gives a significant amelioration for CC($n_m$) models (see Table II). However, as the near degeneracies decrease, closed-shell approaches become the better choice in our molecular series (see Tables III and IV). A systematic difference between cs and os approaches remains even at very high excitation ranks, due to the fact that the different valence models lead to different polarization of the core and virtual spinors. Since the spinor basis is truncated both in the occupied and virtual space, the two models do not yield identical results in the FCI/FCC limit. In the cs case the spinors are optimized for the reference determinant used in the correlated approach. In contrast to this, os spinors
<figure><img src="image_3.png" /><figcaption>FIG. 2. Convergence of various closed-shell CC models for the three molecules. *T* values are *T*<sub>v</sub> taken from Tables II, III, and IV; experimental values are *T*<sub>e</sub><sup>expt</sup> from [37].</figcaption></figure>

comprise an averaging over several states and therefore do
not correspond to the reference state used in the correlation
approach. Due to this inconsistency, we consider os-CC as a
pragmatic approach in certain cases, but csCC results at high
excitation ranks our qualitatively best values.

V. CONCLUSION

An implementation of a general excitation rank relativistic
coupled cluster is presented with which electronically excited
states can be calculated at high accuracy using linear response
theory. It has been demonstrated that the relativistic GAS-CC
approach is applicable to atomic and molecular electronically
excited states, for which we have chosen showcase systems
exhibiting strong effects of both relativistic and electron

correlation origin. We regard these findings largely as proof of
principle for our method.

We conclude from the present study that within the GAS-CC approach both the multireference character and the importance of dynamic electron correlation on relative energies can be addressed efficiently. The former is achieved by adding active-space selected higher excitations to the standard CC expansion. For BiH (and to some degree also SbH) where the ground state is dominated by a single Slater determinant in the relativistic picture the quality of the GAS-CC results surpasses that of a linear wave-function expansion such as relativistic CI theory, even if the latter is applied as a genuine multireference approach. In cases where our chosen Fermi vacuum determinant is no longer the dominant contributor to the electronic ground state (Si atom, AsH, to some degree SbH) we find that higher CC excitations, at least up to full triples,
<figure><img src="image_4.png" /><figcaption>FIG. 3. Deviation from experiment [37] in percent for MRCI and closed-shell CC models for the three molecules. Values are taken from Tables II, III, and IV. $\Delta T$ values are calculated from $T_v$ taken from Tables II, III, and IV and experimental values $T_e^{\text{expt.}}$ from [37]. MRCI models are built according to Table X with $a = 2, b = 4$ for SD\_CISD, $a = 1, b = 3$ for SDT\_CISDT, $a = 0, b = 2$ for SDTQ\_SDTQ, and $a = 0, b = 1$ for SDTQ\_SDTQP.</figcaption></figure>

have to be included for achieving high accuracy. In such cases true multireference CC (such as Mukherjee's Mk-CC [46]) where a number of reference determinants is treated on equal footing would seem to be the better choice. It is planned to implement such a genuine MR approach into our relativistic methodology.

In ongoing work we are generalizing a computationally more efficient commutator-based evaluation of the CC Jacobian matrix to the four-component relativistic formalism. This improvement will lead to a code with the optimal computational scaling of conventional CC theory also in the calculation of excited states, and will allow us to increase the number of explicitly correlated electrons. On the technical side this is carried out by merging the relativistic commutator-based GAS-CC [25] with the approach described in this paper and including the new and more efficient code for the relativistic CC Jacobian.

## ACKNOWLEDGMENTS

This work was granted access to the HPC resources of CALMIP under the allocation 2010-p1050 and 2011-p1050. We furthermore acknowledge technical support by Anthony Scemama with the calculations on the local computing cluster. L.K.S. acknowledges the Villum Fonden for financial support. We thank Radovan Bast for some helpful comments on the calculation of expectation values with Dirac.

## APPENDIX: TECHNICAL DETAILS ON ACTIVE SPINOR SPACES

### 1. The Si atom

For the silicon atom four different correlation model hierarchies are defined. We use three GAS for the active spinors (see Table V with 3s spinors in GAS I, $3p_{1/2}$ spinors in GAS II,
<table><caption>TABLE V. General active space models for Si with three GAS for the standard CC hierarchy. “Min. el.” represents the minimum accumulated number and “Max. el.” the maximum accumulated number of electrons after consideration of a given GAS. $a = 0, 1$ for SD2CC and S2CC, respectively. $b \in \{a, \dots, 2\}$ for CCSDTQ (FCC), CCSDT, and CCSD, respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>1</td><td>a</td><td>2</td><td>3s</td></tr><tr><td>II</td><td>2</td><td>b</td><td>4</td><td>3p<sub>1/2</sub></td></tr><tr><td>III</td><td>X</td><td>4</td><td>4</td><td>3p<sub>3/2</sub>, 3p<sub>3/2</sub> +<br>virtual Kr. pairs</td></tr></tbody></table>

and virtual spinors in GAS III). This first specification allows for defining the standard CC hierarchy and another hierarchy where only up to one hole in the space of the 3s spinors is allowed (S2CC). The third specification allows for the definition of various CC($n_m$) correlation models, for which we use four GAS (see Table VI). This particular GAS structure accounts for a selected set of higher excitations, here up to quadruple excitations which decompose, for instance, in the case of CC(4$_2$) into double excitations from 3s to the virtual spinors combined with double excitations from the 3p$_{1/2}$ to the 3p$_{3/2}$ spinors. Finally, for MRCI calculations, which are genuine MR calculations, we use a specific GAS configuration as shown in Table VII, since here there is no need to define a Fermi vacuum determinant.

## 2. The pnictogen monohydrides

In the molecular cases we use three different correlation model hierarchies. The standard CC series (CCSD through CCSDTQP) is defined by two GAS for the active spinors (see GAS Table VIII). For the CC($n_m$) correlation methods, the minimum number of GAS required is four (see GAS Table IX). Similar to the atomic case, the more pronounced the MR character of the state in question, the more important the amplitude $r_{\pi_{1/2}\pi_{3/2}}^{\pi_{3/2}\pi_{3/2}}$ becomes. The CC($n_m$) models partially account for the MR character by introducing such higher excitations which are expected to give large contributions to the states in question. We finally use an extra GAS configuration suited for MRCI calculations (see GAS Table X). These complete-active-space (CAS)-MRCI calculations were performed to provide results from more standard approaches which are compared with CC models.

<table><caption>TABLE VI. Si general active space models with four GAS for the CC($n_m$) hierarchy (see Fig. 1). $a = 1$ for S2CC(3$_2$). if $a = 0$: $b = 1, 2$ for SD2CC(4$_3$) and SD2CC(4$_2$), respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>1</td><td>a</td><td>2</td><td>3s</td></tr><tr><td>II</td><td>2</td><td>a</td><td>4</td><td>3p<sub>1/2</sub></td></tr><tr><td>III</td><td>4</td><td>b</td><td>4</td><td>3p<sub>3/2</sub>, 3p<sub>3/2</sub></td></tr><tr><td>IV</td><td>X</td><td>4</td><td>4</td><td>Virtual Kr. pairs</td></tr></tbody></table>

<table><caption>TABLE VII. Si general active space models with three GAS for MRCI hierarchy. All are MRCISD2_CAS2in3 type. $a = 0, 1, 2$ for SDTQ-4 (FCI-4), SDT-4, and SD-4, respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>1</td><td>0</td><td>2</td><td>3s</td></tr><tr><td>II</td><td>4</td><td>a</td><td>4</td><td>3p<sub>1/2</sub>, 3p<sub>3/2</sub>, 3p<sub>3/2</sub></td></tr><tr><td>III</td><td>X</td><td>4</td><td>4</td><td>Virtual Kr. pairs</td></tr></tbody></table>

## 3. Virtual spinor spaces

GAS-CCSD and CAS-MRCI calculations were performed with increasing sizes of virtual spinor spaces. It is a standard procedure in four-component electronic-structure calculations with uncontracted Gaussian basis sets to use a truncation energy value for the virtual spinors (see, e.g., Ref. [47]) and to perform the correlation calculation in the resulting subspace. We have in all cases converged the excitation energies with respect to this subspace dimension using the CCSD model.

## 4. Coupling pictures and determinants

### a. L-S coupling.

For the sake of simplicity we adopt a two-particle approximation, i.e., we restrict ourselves to the electronic configuration $np^2$. All states will be written as $(2S+1)L_J(m_J)$ in accord with the Russell-Saunders convention, and determinants as $|L_{m_L}m_S L_{m_L}m_S|$. In addition, we will use the shorthand notation $|\alpha\rangle = |S = \frac{1}{2}, m_S = \frac{1}{2}\rangle$ and $|\beta\rangle = |S = \frac{1}{2}, m_S = -\frac{1}{2}\rangle$ for spin states, and $|P_+\rangle = |L = 1, m_{L=1}\rangle$, $|P_0\rangle = |L = 1, m_{L=0}\rangle$, and $|P_-\rangle = |L = 1, m_{L=-1}\rangle$.

In order to find the expansion of states $(2S+1)L_J(m_J)$ in terms of determinants $|L_{m_L}m_S L_{m_L}m_S|$ we start out from the state with $\max(m_J)$ for a given $(2S+1)L_J$ and apply shift operators $\hat{J}_-$ to construct the states up to $\min(m_J)$. All $m_J$ states are individually normalized. We obtain

$$^1S_0: \qquad ^1S_0(0) = \frac{1}{\sqrt{3}}(|P_-\alpha P_+\beta| - |P_0\alpha P_0\beta| + |P_+\alpha P_-\beta|). \quad (A1)$$

$$^1D_2: \qquad ^1D_2(2) = |P_+\alpha P_+\beta|, \qquad (A2)$$

$$^1D_2(1) = \frac{1}{\sqrt{2}}(|P_0\alpha P_+\beta| + |P_+\alpha P_0\beta|), \qquad (A3)$$

$$^1D_2(0) = \frac{1}{\sqrt{6}}(|P_-\alpha P_+\beta| + 2|P_0\alpha P_0\beta| + |P_+\alpha P_-\beta|), \quad (A4)$$

<table><caption>TABLE VIII. AsH, SbH, and BiH general active space models with two GAS for the standard CC hierarchy. $a = 4, 3, 2, 1$ for CCSD-6, CCSDT-6, CCSDTQ-6, and CCSDTQP-6, respectively. $n = 4, 5, 6$ for AsH, SbH, and BiH, respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>3</td><td>a</td><td>6</td><td>ns, σ<sub>1/2</sub>, π<sub>1/2</sub></td></tr><tr><td>II</td><td>X</td><td>6</td><td>6</td><td>π<sub>3/2</sub> + virtual Kr. pairs</td></tr></tbody></table>
<table><caption>TABLE IX. AsH, SbH, and BiH general active space models with four GAS for the SD4_CC(n<sub>m</sub>) hierarchy (see Fig. 1). a = 3,4 for CC(4<sub>3</sub>) and CC(4<sub>2</sub>), respectively. n = 4,5,6 for AsH, SbH, and BiH, respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>2</td><td>2</td><td>4</td><td>ns, σ<sub>1/2</sub></td></tr><tr><td>II</td><td>3</td><td>2</td><td>6</td><td>π<sub>1/2</sub></td></tr><tr><td>III</td><td>3</td><td>a</td><td>6</td><td>π<sub>3/2</sub></td></tr><tr><td>IV</td><td>X</td><td>6</td><td>6</td><td>Virtual Kr. pairs</td></tr></tbody></table>

$$ ^1D_2(-1) = \frac{1}{\sqrt{2}}(|P_0\alpha P_-\beta| + |P_-\alpha P_0\beta|), \quad (A5) $$

$$ ^3P_0: \qquad ^1D_2(-2) = |P_-\alpha P_-\beta|. \qquad (A6) $$

$$ ^3P_0(0) = \frac{1}{\sqrt{3}} \left[ |P_0\alpha P_-\alpha| + |P_+\beta P_0\beta| - \frac{1}{\sqrt{2}}(|P_+\alpha P_-\beta| + |P_+\beta P_-\alpha|) \right]. \quad (A7) $$

$$ ^3P_1: \qquad ^3P_1(1) = \frac{1}{2}(\sqrt{2}|P_+\alpha P_-\alpha| - |P_+\alpha P_0\beta| - |P_+\beta P_0\alpha|), \quad (A8) $$

$$ ^3P_1(0) = \frac{1}{\sqrt{2}}(|P_0\alpha P_-\alpha| - |P_+\beta P_0\beta|), \quad (A9) $$

$$ ^3P_1(-1) = \frac{1}{2}(-\sqrt{2}|P_+\beta P_-\beta| + |P_0\alpha P_-\beta| + |P_0\beta P_-\alpha|). \quad (A10) $$

$$ ^3P_2: \qquad ^3P_2(2) = |P_+\alpha P_0\alpha|, \quad (A11) $$

$$ ^3P_2(1) = \frac{1}{2}(\sqrt{2}|P_+\alpha P_-\alpha| + |P_+\alpha P_0\beta| + |P_+\beta P_0\alpha|), \quad (A12) $$

$$ ^3P_2(0) = \frac{1}{\sqrt{6}}[|P_0\alpha P_-\alpha| + |P_+\beta P_0\beta| + \sqrt{2}(|P_+\alpha P_-\beta| + |P_+\beta P_-\alpha|)], \quad (A13) $$

$$ ^3P_2(-1) = \frac{1}{2}(\sqrt{2}|P_+\beta P_-\beta| + |P_0\alpha P_-\beta| - |P_0\beta P_-\alpha|), \quad (A14) $$

$$ ^3P_2(-2) = |P_0\beta P_-\beta|. \quad (A15) $$

<table><caption>TABLE X. AsH, SbH, and BiH general active space models with three GAS for MRCI hierarchy. a = 0,1,2 for SDTQ-4, SDT-4, and SD-4, respectively. b ≥ a,b = 1,2,3,4 for CISDTQP-6, CISDTQ-6, CISDT-6, and CISD-6, respectively. n = 4,5,6 for AsH, SbH, and BiH, respectively. X: Number of virtual Kramers pairs.</caption><thead><tr><th>GAS</th><th>Kramers pairs<br>per irrep. E<sub>1/2</sub></th><th>Min. el.</th><th>Max. el.</th><th>Shell types</th></tr></thead><tbody><tr><td>I</td><td>2</td><td>a</td><td>4</td><td>ns, σ<sub>1/2</sub></td></tr><tr><td>II</td><td>4</td><td>b</td><td>6</td><td>π<sub>1/2</sub>, π<sub>3/2</sub></td></tr><tr><td>III</td><td>X</td><td>6</td><td>6</td><td>Virtual Kr. pairs</td></tr></tbody></table>

To find the expansion of a given determinant in terms of $(2S+1)L_J$ states we have to invert the matrix $X$ in

$$ \vec{s} = X\vec{d}, \quad (A16) $$

where $\vec{s}$ is a vector of states and $\vec{d}$ is a vector of determinants. $X$ is orthonormal, so solving

$$ \vec{d} = X^{-1}\vec{s} \quad (A17) $$

is easy since $X^{-1} = X^T$.

We then find that the various determinants for the subspace $m_J = 0$ can be expressed as

$$ |P_0\alpha P_-\alpha| = \frac{1}{\sqrt{3}} {}^3P_0 + \frac{1}{\sqrt{2}} {}^3P_1 + \frac{1}{\sqrt{6}} {}^3P_2, \quad (A18) $$

$$ |P_+\beta P_0\beta| = \frac{1}{\sqrt{3}} {}^3P_0 - \frac{1}{\sqrt{2}} {}^3P_1 + \frac{1}{\sqrt{6}} {}^3P_2, \quad (A19) $$

$$ |P_-\alpha P_+\beta| = \frac{1}{\sqrt{3}} {}^1S_0 + \frac{1}{\sqrt{6}} {}^1D_2 + \frac{1}{\sqrt{6}} {}^3P_0 - \frac{1}{\sqrt{3}} {}^3P_2, \quad (A20) $$

$$ |P_+\alpha P_-\beta| = \frac{1}{\sqrt{3}} {}^1S_0 + \frac{1}{\sqrt{6}} {}^1D_2 - \frac{1}{\sqrt{6}} {}^3P_0 + \frac{1}{\sqrt{3}} {}^3P_2, \quad (A21) $$

$$ |P_0\alpha P_0\beta| = -\frac{1}{\sqrt{3}} {}^1S_0 + \frac{2}{\sqrt{6}} {}^1D_2. \quad (A22) $$

Notice sign changes for determinants such as $|P_-\alpha P_+\beta|$. For the other $m_J$ values we find the following expressions:

$m_J = 1:$

$$ |P_+\alpha P_-\alpha| = \frac{1}{\sqrt{2}} {}^3P_1 + \frac{1}{\sqrt{2}} {}^3P_2, \quad (A23) $$

$$ |P_+\alpha P_0\beta| = \frac{1}{\sqrt{2}} {}^1D_2 - \frac{1}{2} {}^3P_1 + \frac{1}{2} {}^3P_2, \quad (A24) $$

$$ |P_+\beta P_0\alpha| = -\frac{1}{\sqrt{2}} {}^1D_2 - \frac{1}{2} {}^3P_1 + \frac{1}{2} {}^3P_2. \quad (A25) $$

$m_J = -1:$

$$ |P_+\beta P_-\beta| = -\frac{1}{\sqrt{2}} {}^3P_1 + \frac{1}{\sqrt{2}} {}^3P_2, \quad (A26) $$

$$ |P_0\alpha P_-\beta| = \frac{1}{\sqrt{2}} {}^1D_2 + \frac{1}{2} {}^3P_1 + \frac{1}{2} {}^3P_2, \quad (A27) $$

$$ |P_0\beta P_-\alpha| = -\frac{1}{\sqrt{2}} {}^1D_2 + \frac{1}{2} {}^3P_1 + \frac{1}{\sqrt{2}} {}^3P_2. \quad (A28) $$

$m_J = 2:$

$$ |P_+\alpha P_0\alpha| = {}^3P_2, \quad (A29) $$

$$ |P_+\alpha P_+\beta| = {}^1D_2. \quad (A30) $$

$m_J = -2:$

$$ |P_0\beta P_-\beta| = {}^3P_2, \quad (A31) $$

$$ |P_-\alpha P_-\beta| = {}^1D_2. \quad (A32) $$
b. J-J coupling

First we want to find the various states $(J_1, J_2)_J$ that arise
from J-J coupled spinors. These five states are

$$
(3/2,3/2)_2, (3/2,3,2)_0, (3/2,1/2)_2, (3/2,1/2)_1, (1/2,1/2)_0. \tag{A33}
$$

We expand each spinor $(J, m_J)$ in terms of nonrelativistic
spin orbitals using the notation introduced above and the
corresponding Clebsch-Gordan coefficients:

$$
J = 3/2:
$$

$$
(3/2, 1/2) = \frac{1}{\sqrt{3}}(\sqrt{2}P_0\alpha + P_+\beta), \quad (\text{A35})
$$

$$
(3/2, -1/2) = \frac{1}{\sqrt{3}}(\sqrt{2}P_0\beta + P_{-\alpha}), \quad (\text{A36})
$$

$$
(3/2, -3/2) = P_{-\beta}. \tag{A37}
$$

$J = 1/2:$

$$
(1/2, 1/2) = \frac{1}{\sqrt{3}}(-P_0\alpha + \sqrt{2}P_+\beta), \quad (\text{A38})
$$

$$
(1/2, -1/2) = \frac{1}{\sqrt{3}}(P_0\beta - \sqrt{2}P_{-\alpha}). \quad (\text{A39})
$$

This allows us to write the $m_J$ components of the various
$(J_1, J_2)_J$ states, denoted as $(J_1, J_2)_{J,m_J}$, as

$$
(1/2, 1/2)_{0,0} = |(1/2, 1/2), (1/2, -1/2)| = \frac{1}{3}(-|P_0\alpha P_0\beta| + \sqrt{2}|P_+\beta P_0\beta| + \sqrt{2}|P_0\alpha P_{-\alpha}| + 2|P_{-\alpha}P_+\beta|), \quad (\text{A40})
$$

$$
(3/2, 3/2)_{2,2} = |(3/2, 3/2), (3/2, 1/2)| = \frac{1}{\sqrt{3}}(\sqrt{2}|P_+\alpha P_0\alpha| + |P_+\alpha P_+\beta|), \quad (A41)
$$

$$
(3/2, 3/2)_{2,1} = |(3/2, 3/2), (3/2, -1/2)| = \frac{1}{\sqrt{3}}(\sqrt{2}|P_+\alpha P_0\beta| + |P_+\alpha P_{-\alpha}|), \quad (A42)
$$

$$
(3/2, 3/2)_{2,0} = \frac{1}{\sqrt{2}} \left[ \begin{aligned}
& |(3/2, 1/2), (3/2, -1/2)| + |(3/2, 3/2), (3/2, -3/2)| \\
& \qquad = \frac{1}{\sqrt{2}} \left[ \frac{1}{3} \left( \begin{aligned}[t]
& 2|P_0\alpha P_0\beta| + \sqrt{2}|P_0\alpha P_{-\alpha}| \\
& + \sqrt{2}|P_+\beta P_0\beta| + |P_{-\alpha}P_+\beta| \end{aligned} \right) + |P_+\alpha P_{-\alpha}| \right]
\end{aligned} \right], \quad (\text{A43})
$$

$$
(3/2, 3/2)_{2,-1} = |(3/2, 1/2), (3/2, -3/2)| = \frac{1}{\sqrt{3}}(-\sqrt{2}|P_{-\beta}P_0\alpha| - |P_{-\beta}P_+\beta|), \quad (\text{A44})
$$

$$
(3/2, 3/2)_{2,-2} = |(3/2, -1/2), (3/2, -3/2)| = \frac{1}{\sqrt{3}}(-\sqrt{2}|P_{-\beta}P_0\beta| - |P_{-\beta}P_{-\alpha}|), \quad (\text{A45})
$$

$$
\begin{align}
(3/2, 3/2)_{0,0} &= \frac{1}{\sqrt{2}} \left( |(3/2, 1/2), (3/2, -1/2)| - |(3/2, 3/2), (3/2, -3/2)| \right) \nonumber \\
&= \frac{1}{\sqrt{2}} \left[ \frac{1}{3} \left( 2|P_0\alpha P_0\beta| + \sqrt{2}|P_0\alpha P_{-\alpha}| + \sqrt{2}|P_+\beta P_0\beta| + |P_{-\alpha}P_+\beta| \right) - |P_+\alpha P_{-\alpha}| \right], \tag{A46}
\end{align}
$$

$$
(3/2, 1/2)_{2,2} = |(3/2, 3/2), (1/2, 1/2)| = \frac{1}{\sqrt{3}}(-|P_+\alpha P_0\alpha| + \sqrt{2}|P_+\alpha P_+\beta|), \quad (\text{A47})
$$

$$
\begin{align}
(3/2, 1/2)_{2,1} &= \frac{1}{2}[\sqrt{3}|(3/2, 1/2), (1/2, 1/2)| + |(3/2, 3/2), (1/2, -1/2)|] \nonumber \\
&= \frac{1}{2}\left(\sqrt{3}|P_0\alpha P_+\beta| + \frac{1}{\sqrt{3}}|P_+\alpha P_0\beta| - \sqrt{\frac{2}{3}}|P_+\alpha P_{-\alpha}|\right), \tag{A48}
\end{align}
$$

$$
(3/2, 1/2)_{2,0} = \frac{1}{\sqrt{2}} [ |(3/2, 1/2), (1/2, -1/2)| + |(3/2, -1/2), (1/2, 1/2)| ] = \frac{1}{3\sqrt{2}} (2\sqrt{2}|P_0\alpha P_0\beta| - |P_0\alpha P_{-\alpha}| \\
- |P_+\beta P_0\beta| + 2\sqrt{2}|P_{-\alpha}P_+\beta| ), \qquad (\text{A49})
$$

$$
\begin{align}
(3/2, 1/2)_{2,-1} &= \frac{1}{2} [\sqrt{3} |(3/2, -1/2), (1/2, -1/2)| + |(3/2, -3/2), (1/2, 1/2)|] \nonumber \\
&= \frac{1}{2} \left( -\sqrt{3} |P_0\beta P_{-\alpha}| - \frac{1}{\sqrt{3}} |P_{-\beta} P_0\alpha| + \sqrt{\frac{2}{3}} |P_{-\beta} P_+\beta| \right), \tag{A50}
\end{align}
$$

$$
(3/2, 1/2)_{2,-2} = |(3/2, -3/2), (1/2, -1/2)| = \frac{1}{\sqrt{3}}(|P_{-\beta}P_0\beta| - \sqrt{2}|P_{-\beta}P_{-\alpha}|), \quad (\text{A51})
$$

$$
(3/2, 1/2)_{1,1} = \frac{1}{2} [ |(3/2, 1/2), (1/2, 1/2)| - \sqrt{3} |(3/2, 3/2), (1/2, -1/2)| ] = \frac{1}{2} ( |P_0\alpha P_+\beta| - |P_+\alpha P_0\beta| + \sqrt{2} |P_+\alpha P_{-\alpha}| ), \quad (\text{A52})
$$

$$
(3/2, 1/2)_{1,0} = \frac{1}{\sqrt{2}} [ -|(3/2, 1/2), (1/2, -1/2)| + |(3/2, -1/2), (1/2, 1/2)| ] = \frac{1}{\sqrt{2}} ( |P_0\alpha P_{-\alpha}| - |P_+\beta P_0\beta| ), \quad (\text{A53})
$$
$$ (3/2, 1/2)_{1,-1} = \frac{1}{2}[-|(3/2, -1/2), (1/2, -1/2)| + \sqrt{3}|(3/2, -3/2), (1/2, 1/2)|] \\ = \frac{1}{2}(|P_0\beta P_{-\alpha}| - |P_{-\alpha} P_0\beta| + \sqrt{2}|P_{-\beta} P_{+\beta}|). \qquad (A54) $$

### c. J-J coupled states in terms of L-S coupled states and back

Writing the $(1/2, 1/2)_0$ state out in terms of the determinants from Eqs. (A18)–(A32) we find

$$ (1/2, 1/2)_0 = \frac{1}{\sqrt{3}} {}^1S_0 + \sqrt{\frac{2}{3}} {}^3P_0. \qquad (A55) $$

The remaining states can accordingly be expressed as

$$ (3/2, 3/2)_2 = \sqrt{\frac{2}{3}} {}^3P_2 + \frac{1}{\sqrt{3}} {}^1D_2, \qquad (A56) $$

$$ (3/2, 3/2)_0 = \frac{1}{\sqrt{3}} {}^3P_0 - \sqrt{\frac{2}{3}} {}^1S_0, \qquad (A57) $$

$$ (3/2, 1/2)_2 = \sqrt{\frac{2}{3}} {}^1D_2 - \frac{1}{\sqrt{3}} {}^3P_2, \qquad (A58) $$

$$ (3/2, 1/2)_1 = {}^3P_1. \qquad (A59) $$

We close this section with the inverse expansion of Russell-Saunders terms in terms of J-J coupled states:

$$ {}^1S_0 = \frac{1}{\sqrt{3}}(1/2, 1/2)_0 - \sqrt{\frac{2}{3}}(3/2, 3/2)_0, \qquad (A60) $$

$$ {}^3P_0 = \sqrt{\frac{2}{3}}(1/2, 1/2)_0 + \frac{1}{\sqrt{3}}(3/2, 3/2)_0, \qquad (A61) $$

$$ {}^3P_1 = (3/2, 1/2)_1, \qquad (A62) $$

$$ {}^3P_2 = \sqrt{\frac{2}{3}}(3/2, 3/2)_2 - \frac{1}{\sqrt{3}}(3/2, 1/2)_2, \qquad (A63) $$

$$ {}^1D_2 = \frac{1}{\sqrt{3}}(3/2, 3/2)_2 + \sqrt{\frac{2}{3}}(3/2, 1/2)_2. \qquad (A64) $$

The expansion of determinants over relativistic spinors in terms of Russell-Saunders terms—as referred to in the main body of the paper—can be deduced by combining Eqs. (A40)–(A54) with Eqs. (A60)–(A64).

## 5. Molecular determinants and states—Choice of Fermi vacuum

### a. Wave function of molecular states

For the lightest homolog spin-orbit interaction is a perturbation to electrostatic effects. Furthermore, it is known from similar systems with an approximate valence $\pi^2$ configuration that $\sigma-\pi$ mixing due to spin-orbit interaction is negligibly small in $4p$ element molecules [48]. We therefore start out from a molecular two-electron valence wave function for the expected electronic ground state which can be written as

$$ |{}^3\Sigma_0\rangle = c_0|{}^3\Sigma_{M_j=0}\rangle + c'|{}^1\Sigma_{M_j=0}\rangle, \qquad (A65) $$

the $|{}^1\Sigma\rangle$ state corresponding to the $\pi^2$ configuration being the main perturber. We estimate the mixing coefficient $c'$ for a first-order perturbation correction to the wave function with the one-electron spin-orbit Hamiltonian in Pauli approximation:

$$ c' = \frac{\langle {}^3\Sigma | \hat{H}^{SO} | {}^1\Sigma \rangle}{E_1\Sigma - E_3\Sigma} \approx \frac{151}{7050}. \qquad (A66) $$

The value of $151 \text{ cm}^{-1}$ has been obtained by using an effective nuclear charge of 7.44 a.u. for a $4p$ electron in As [49] and an expectation value $\langle \frac{1}{\vec{r}} \rangle = 7.0$ a.u. from Ref. [50] for calculating the spin-orbit matrix element. The energy difference of $7050 \text{ cm}^{-1}$ has been calculated using the LUCITA module of the DIRAC program package [29] in Dyall’s spin-orbit free approximation [51] to the Dirac-Coulomb Hamiltonian.

Normalizing the total wave function thus gives us an estimated contribution of roughly 0.03% of the $|{}^1\Sigma_{M_j=0}\rangle$ to the molecular ground state, which can safely be neglected, even if two-electron spin-orbit contributions were accounted for in addition. This means that the nonrelativistic $|{}^{2S+1}\Lambda_{\pm\Omega}\rangle$ wave functions

$$ |{}^3\Sigma_0\rangle = \frac{1}{2}[\pi_+(1)\pi_-(2) - \pi_-(1)\pi_+(2)][\alpha(1)\beta(2) + \beta(1)\alpha(2)], \qquad (A67) $$

$$ |{}^1\Sigma_0\rangle = \frac{1}{2}[\pi_+(1)\pi_-(2) + \pi_-(1)\pi_+(2)][\alpha(1)\beta(2) - \beta(1)\alpha(2)], \qquad (A68) $$

are a good approximation to the molecular states in AsH, where we use the notation $[symbol]_{m_i}(j)$ denoting $\lambda_{m_i}(\vec{r}_j)$ for the spatial wave function of particle j and the spin part in accordance with the definition in Sec. A 4 a.

We finally expand the results in Eqs. (A67) and (A68) into Cartesian components, according to $\pi_+ = -\frac{1}{\sqrt{2}}(\pi_x + i\pi_y)$ and $\pi_- = \frac{1}{\sqrt{2}}(\pi_x - i\pi_y)$, yielding

$$ |{}^3\Sigma_0\rangle = \frac{i}{2}[\pi_x(1)\pi_y(2) - \pi_y(1)\pi_x(2)][\alpha(1)\beta(2) + \beta(1)\alpha(2)], \qquad (A69) $$

$$ |{}^1\Sigma_0\rangle = -\frac{1}{2}[\pi_x(1)\pi_x(2) + \pi_y(1)\pi_y(2)][\alpha(1)\beta(2) - \beta(1)\alpha(2)]. \qquad (A70) $$

### b. Choice of Fermi vacuum

In order to compare our Fermi vacuum state to the estimated molecular wave functions we have carried out a Mulliken population analysis of the AsH valence spinors as a function of internuclear distance. The results are to be found in Figure 4. The spinors underlying the figure are those energetically highest and still doubly occupied (HOMS). At the equilibrium bond length they are energetically well separated from the bonding spinors by 0.125 a.u. We therefore construct our
<figure><img src="image_5.png" /><figcaption>FIG. 4. Mulliken population analysis of the HOMS spinor $m_j = \frac{1}{2}$ as a function of internuclear distance. Since the closed-shell DCHF model does not lead to physically correct dissociation, we exploit the information from close to the equilibrium bond distance only.</figcaption></figure>

two-electron Fermi vacuum state from the HOMS

$$ (m_j = +\frac{1}{2}) = \frac{1}{\sqrt{2}}(-\pi_x - i\pi_y)\beta, $$

$$ (m_j = -\frac{1}{2}) = \frac{1}{\sqrt{2}}(+\pi_x - i\pi_y)\alpha, $$

where we have represented the molecular spinors by their principal character. The form of the spatial part has been obtained from the MO-AO expansion coefficients of the Dirac-Coulomb Hartree-Fock calculation, and the spin function from computing the expectation value $\langle\varphi_{j,m_j}|\hat{s}_z|\varphi_{j,m_j}\rangle$ for the respective spinors $\varphi_{j,m_j}$. The Kramers partner has been deduced by applying the time-reversal operator to a given spinor.

Using this information we can rewrite our Fermi vacuum state as

$$ \begin{aligned} &|(m_j)_1; (m_j)_2| \\ &= \left|\left(m_j = \frac{1}{2}\right); \left(m_j = -\frac{1}{2}\right)\right| \\ &= \frac{1}{2\sqrt{2}} \left\{ \begin{aligned} &[-\pi_x(1) - i\pi_y(1)]\beta(1)[\pi_x(2) - i\pi_y(2)]\alpha(2) \\ &+ [\pi_x(2) + i\pi_y(2)]\beta(2)[\pi_x(1) - i\pi_y(1)]\alpha(1) \end{aligned} \right\} \\ &= \frac{i}{2\sqrt{2}} [\pi_x(1)\pi_y(2) - \pi_y(1)\pi_x(2)][\alpha(1)\beta(2) + \beta(1)\alpha(2)] \\ &\quad + \frac{1}{2\sqrt{2}} [\pi_x(1)\pi_x(2) + \pi_y(1)\pi_y(2)][\alpha(1)\beta(2) \\ &\quad - \beta(1)\alpha(2)]. \end{aligned} \quad (A71) $$

Comparing Eqs. (A67) and (A68) with Eq. (A71) shows that our Fermi vacuum state from a relativistic calculation represents the true ground state only to roughly 50% and contains an equally large contribution from the excited $|^1\Sigma_0\rangle$ state.

[1] J. Doyle, B. Friedrich, R. V. Krems, and F. Masnou-Seeuws, Eur. Phys. J. D 31, 149 (2004).

[2] W. C. Stwalley, P. L. Gould, and E. E. Eyler, in Cold Molecules, edited by R. V. Krems, W. C. Stwalley, and B. Friedrich (CRC Press, Boca Raton, 2009), Chap. 5.

[3] M. Asplund, N. Grevesse, A. J. Sauval, and P. Scott, Annu. Rev. Astron. Astrophys. 47, 481 (2009).

[4] P. S. Barklem, A. K. Belyaev, M. Guitou, N. Feautrier, F. X. Gadéa, and A. Spielfiedel, Astron. Astrophys. 530, A94 (2011).

[5] E. D. Commins, Adv. At., Mol., Opt. Phys. 40, 1 (1999).

[6] A. E. Leanhardt, J. L. Bohn, H. Loh, P. Maletinsky, E. R. Meyer, L. C. Sinclair, R. P. Stutz, and E. A. Cornell, J. Mol. Spectrosc. 270, 1 (2011).

[7] V. V. Ivanov, D. I. Lyakh, and L. Adamowicz, in Recent Progress in Coupled Cluster Methods, edited by P. Čársky, J. Paldus, and J. Pittner, Vol. 11, (Springer, Heidelberg, Germany, 2010), Chap. 9.

[8] Recent Progress in Coupled Cluster Methods, edited by P. Čársky, J. Paldus, and J. Pittner (Springer, Heidelberg, 2010).

[9] T. Fleig, Chem. Phys. 395, 2 (2012).

[10] L. Visscher, E. Eliav, and U. Kaldor, J. Chem. Phys. 115, 9720 (2001).

[11] E. Eliav, M. J. Vilkas, Y. Ishikawa, and U. Kaldor, J. Chem. Phys. 122, 224113 (2005).

[12] S. Hirata, T. Yanai, R. J. Harrison, M. Kamija, and P.-D. Fang, J. Chem. Phys. 126, 024104 (2007).

[13] J. F. Stanton and R. J. Bartlett, J. Chem. Phys. 98, 7029 (1993).

[14] M. Musial and R. J. Bartlett, J. Chem. Phys. 129, 134105 (2008).

[15] P. A. Christiansen and K. S. Pitzer, J. Chem. Phys. 73, 5160 (1980).

[16] T. Fleig, in Relativistic Methods for Chemists, edited by M. Barysz and Y. Ishikawa, Vol. 10, (Springer, Heidelberg, 2010), Chap. 10.

[17] T. Zeng, D. G. Fedorov, and M. Klobukowski, J. Chem. Phys. 131, 124109 (2009).

[18] S. Knecht, H. J. Aa Jensen, and T. Fleig, J. Chem. Phys. 132, 014108 (2010).

[19] J. Olsen and P. Jørgensen, J. Chem. Phys. 82, 3235 (1984).

[20] H. Koch and P. Jørgensen, J. Chem. Phys. 93, 3333 (1990).

[21] K. Emrich, Nucl. Phys. A 351, 379 (1981).

[22] K. Hald, P. Jørgensen, J. Olsen, and M. Jaszuński, J. Chem. Phys. 115, 671 (2001).

[23] T. Fleig, L. K. Sørensen, and J. Olsen, Theor. Chem. Acc. 118, 347 (2007); 118, 979(E) (2007).

[24] L. K. Sørensen, T. Fleig, and J. Olsen, Z. Phys. Chem. 224, 671 (2010).

[25] L. K. Sørensen, J. Olsen, and T. Fleig, J. Chem. Phys. 134, 214102 (2011).
[26] T. Fleig, J. Olsen, and C. M. Marian, *J. Chem. Phys.* **114**, 4775 (2001).

[27] T. Saue and H. J. Aa Jensen, *J. Chem. Phys.* **111**, 6211 (1999).

[28] L. Visscher and T. Saue, *J. Chem. Phys.* **113**, 3996 (2000).

[29] DIRAC, a relativistic *ab initio* electronic structure program, Release DIRAC10 (2010), written by T. Saue, L. Visscher, and H. J. Aa Jensen, with new contributions from R. Bast, K. G. Dyall, U. Ekstrøm, E. Eliav, T. Enevoldsen, T. Fleig, A. S. P. Gomes, J. Henriksson, M. Iliaš, Ch. R. Jacob, S. Knecht, H. S. Nataraj, P. Norman, J. Olsen, M. Pernpointner, K. Ruud, B. Schimmelpfennig, J. Sikkema, A. Thorvaldsen, J. Thyssen, S. Villaume, and S. Yamamoto.

[30] T. Fleig, J. Olsen, and L. Visscher, *J. Chem. Phys.* **119**, 2963 (2003).

[31] O. Christiansen, H. Koch, P. Jørgensen, and J. Olsen, *Chem. Phys. Lett.* **256**, 185 (1996).

[32] J. Olsen, *J. Chem. Phys.* **113**, 7140 (2000).

[33] B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov, and P. O. Widmark, *J. Phys. Chem. A* **108**, 2851 (2005).

[34] K. G. Dyall, *Theor. Chim. Acta* **109**, 284 (2003).

[35] K. G. Dyall, *Theor. Chim. Acta* **115**, 441 (2006).

[36] Thom. H. Dunning, *J. Chem. Phys.* **90**, 1007 (1989).

[37] K. P. Huber and G. Herzberg, *Molecular Spectra and Molecular Structure* (Van Nostrand Reinhold, New York, 1979).

[38] K. P. Huber and G. Herzberg, *Constants of Diatomic Molecules* (data prepared by J. W. Gallagher and R. D. Johnson, III)

in *NIST Chemistry WebBook*, edited by P. J. Linstrom and W. G. Mallard, NIST Standard Reference Database Number 69, (National Institute of Standards and Technology, Gaithersburg) http://webbook.nist.gov (retrieved December 10, 2008).

[39] R. J. Bartlett and M. Musial, *Rev. Mod. Phys.* **79**, 291 (2007).

[40] A. Köhn and J. Olsen, *J. Chem. Phys.* **125**, 174110 (2006).

[41] S. Knecht, H. J. Aa Jensen, and T. Fleig, *J. Chem. Phys.* **128**, 014108 (2008).

[42] T. Fleig, H. J. Aa Jensen, J. Olsen, and L. Visscher, *J. Chem. Phys.* **124**, 104106 (2006).

[43] Stefan Knecht (Odense) (unpublished).

[44] J. Loras, S. Knecht, and T. Fleig (unpublished).

[45] Jessica Loras, Master thesis, Laboratory of Quantum Chemistry and Quantum Physics (LCPQ), University Paul Sabatier, Toulouse, 2011.

[46] U. S. Mahapatra, B. Datta, and D. Mukherjee, *J. Chem. Phys.* **110**, 6171 (1999).

[47] T. Fleig, *Phys. Rev. A* **72**, 052506 (2005).

[48] J.-B. Rota, S. Knecht, T. Fleig, D. Ganyushin, T. Saue, F. Neese, and H. Bolvin, *J. Chem. Phys.* **135**, 114106 (2011).

[49] E. Clementi and D. L. Raimondi, *J. Chem. Phys.* **38**, 2686 (1963).

[50] J.-P. Desclaux, *At. Data Nucl. Data Tables* **12**, 311 (1973).

[51] K. G. Dyall, *Introduction to Relativistic Quantum Chemistry* (Chemistry Department, University of Odense, Denmark, 1995).

[52] W. C. Martin and R. Zalubas, *J. Phys. Chem. Ref. Data* **12**, 323 (1983).