S pectroscopic techniques exploiting X-ray radiation have a long history. Two of the most popular ones, X-ray absorption and X-ray photoemission (also known as X-ray photoelectron or electron spectroscopy for chemical analysis), enable investigation of the local electronic structure in molecules and materials. Today's light sources, which range from synchrotron and X-ray free-electron lasers to table-top Xray instruments based on high harmonic generation, facilitate exciting new experiments, which were merely hypothetical just a few years back. These advances have triggered an explosion of interest of the molecular and material sciences community in X-ray based techniques. [1][2][3] These advances in the experimental tools have been accompanied by a burst of activity in the development of theoretical methods for simulating and interpreting experimental spectra. 4 Conceptually, X-ray spectroscopy is similar to UV-vis spectroscopy, the main difference being the energy scale and, consequently, the type of electronic transitions that are probed. UV-vis radiation induces transitions involving the outer-shell valence electrons, whereas X-ray radiation induces transitions involving inner-shell core electrons. Despite this similarity, the theoretical methods developed for valence spectroscopy are not directly applicable to core-level spectroscopies. 4 Similarly to their valence counterparts, core-level states often have openshell character, but they also exhibit strong orbital relaxation. Thus, their description requires sufficiently flexible basis sets 5,6 and many-body ansaẗze that are capable of tackling static and dynamic correlation as well as orbital relaxation.

The equation-of-motion (EOM) coupled-cluster (CC) framework [7][8][9][10][11][12] is a versatile platform for treating excited and ionized states with open-shell character. Even at the lowest level of the correlation treatment, when only single and double excitations are included in the ansatz, the method has the ability to tackle both dynamic correlation and orbital relaxation quite well. Originally, EOM-CC was developed to study valence states; thus, typical implementations seek the solutions corresponding to the lowest states, which is obviously not suitable for high-energy core-level states. Another complication arises from the fact that the core-level states are embedded in a continuum of valence excited and ionized states, which leads to poor convergence and erratic results. 13,14 The core-valence separation (CVS) scheme, 15 which decouples the valence and core sectors of the Fock space on the basis of the large energy gap between the core and valence orbitals, provides a simple yet effective recipe for extending valence-state methods into the core-level domain and overcoming the convergence issues. CVS has been implemented within various electronic structure methods, [16][17][18][19][20][21] including the EOM-CC family. [22][23][24] The resulting methods have been successfully applied to model a variety of X-ray spectroscopic experiments such as absorption (XAS), 22,23,25,26 photoelectron (XPS), 24,27,28 X-ray emission (XES), and resonant inelastic scattering (RIXS). 14,[29][30][31] At the CC level, the applications have been so far limited to the study of 1s electrons of light elements (that is, the K-edge), due to yet another obstacle toward the quantitative simulation of Xray spectrathe need to include relativistic effects.

Relativistic effects, which become more pronounced at higher energies, can be classified into two categories: scalar and spin-orbit. 32 The first type is not critically important in the context of K-edge spectroscopy, because it results in a constant energy shift of the entire spectrum. 5,27,33 Because K-edge states are mostly affected by scalar relativistic effects, nonrelativistic calculations yield qualitatively correct K-edge spectra. The second type, spin-orbit coupling (SOC), which arises from the coupling between the magnetic moment associated with the spin of the electron and the magnetic field created by the relative motion of charged particles (electrons and nuclei), 34 has a greater impact on the spectra at lower edges. SOC mixes states with different multiplicity, which do not interact in a nonrelativistic framework. It also causes energy splittings of orbitals with nonzero orbital angular momentum (l > 0). Figure 1 shows the SO splitting of the atomic 2p orbitals, illustrating that SOC affects the L-edge spectra by splitting them into two edges: L 2 (or L II ) and L 3 (or L III ).

These effects can be fully accounted for within a fully relativistic treatment with a four-component Hamiltonian, [35][36][37] but such treatments come with a substantial increase of the computational cost. Various flavors of two-component methods, such as the zeroth-order regular approximation (ZORA) [38][39][40][41] and its infinite-order variant (IORA), 42 the Douglas-Kroll-Hess method, [43][44][45] or the exact two-component (X2C) approach, 46,47 are less demanding; however, these calculations still increase the computational cost by an order of magnitude relative to the nonrelativistic calculation. 48 Furthermore, a variational treatment of the SOC (i.e., when the spin-orbit operator is included at the wave function optimization step) may lead to an imbalanced description of electronic states with different spin projections, resulting in the violation of Kramers' theorem and momentum contamination. 49 Fortunately, in molecules composed of atoms from the first few rows of the periodic table, so-called perturbative treatment of SOC, which entails calculation of the matrix elements of the Breit-Pauli Hamiltonian using nonrelativistic wave functions, is sufficiently accurate while being computationally affordable. 50 This strategy, which has been successfully used within EOM-CC framework, 48,[51][52][53][54][55] has not yet been extended to core-level spectroscopy, which is the focus of this communication. Building upon our previous work, 23,24,[54][55][56] we implemented calculations of SOC within the frozen core (fc) CVS-EOM-CCSD method, thus extending the EOM-CC framework to modeling X-ray absorption and photoelectron spectroscopy at the L and higher (i.e., lower in energy domain) edges.

In this approach, the final states and their properties (energies and transition strengths) are obtained in a two-step procedure. In the first step, the nonrelativistic states are computed using the fc-CVS-EOM-CCSD ansatz:

where ̂ and ̂ are the cluster and the EOM excitation operators, and the subindices v and c refer to the valence and core orbital spaces, respectively. The exact form of the ̂ operator depends on the target-state manifold, giving rise to different flavors of EOM-CC methods. 9 Here, we use EOM-IP (EOM for ionization potentials) to compute ionized states, EOM-EE (EOM for excitation energies) to compute singlet excited states, and EOM-SF (EOM spin-flip) to compute triplet excited states (with spin projection m s = -1). For closed-shell systems, such as those studied here, these calculations employ closed-shell reference states.

In the second step, SOCs are calculated as the matrix elements of the spin-orbit part of the Breit-Pauli Hamiltonian, given in atomic units by 57

where h SO (i) and h SO (i, j) are the one-and two-electron

r i , p i , and l i are the coordinates and the (linear and angular) momenta, respectively, of electron i, and R K and Z K are the coordinates and charge of nucleus K. In second quantization, eq 2 assumes the following form 54 c I a a J a a a a 1 2 1 2 pq pq p q pqrs pqrs p q s r

where I and J refer to the one-and two-electron spin-orbit integrals:

Matrix elements of this operator can be computed by contracting the spin-orbit integrals with the corresponding one-and two-electron transition density matrices (TDM):

where s m a a s m ( , ) ( , )

s m a a a a s m ( , ) ( , )

The calculation of the SOCs, as defined by eq 2, entails the computationally demanding two-electron part of

Fortunately, the two-electron contributions can be effectively evaluated in a mean-field manner. Within this approximation, called spin-orbit mean-field (SOMF), 59 the calculation of SOC requires only the one-electron TDM

with the effective one-electron operator SOMF ̂ with integrals: 51,52,54

where ρ is the state density matrix of the reference determinant. The SOMF approximation simply entails neglecting the nonseparable part of the two-electron transition density matrix. 54 By converting these equations into the atomic orbital basis, one can evaluate the SOMF integrals using efficient algorithms, such as those used for the Fock-matrix builds. The SOMF Hamiltonian can be used as a starting point for more drastic approximations, such as atomic mean-field approximation in which only the diagonal blocks of ρ (in the AO basis) are retained, and introducing one-center approximations for the one and two-electron integrals. 60 In our implementation, we used the SOMF scheme without further approximations.

To obtain the SO-split states, we then construct the SOMF Hamiltonian matrix in the basis of zeroth-order states whose eigenvalues are the energies of the SO-split states. To compute oscillator strengths for the transitions involving the SO-split states, as needed for the simulation of L-edge NEXAFS, we first construct the (non-Hermitian) electric dipole matrices between the ground state Ψ 0 and the zerothorder target states Ψ(s, m s ):

where α denotes the Cartesian components x, y, and z. These transition matrices are then transformed into the new basis of the SO-split states by applying the transformation obtained from the diagonalization of H SOMF :

where the matrix U contains the eigenvectors of H SOMF . (Due to non-Hermiticity of the EOM theory, the geometric averaging of the estimates from eq 15 does not necessarily yield a real positive number; however, for the systems considered here, only very few nonpositive values were observed and their magnitude was negligible.) Finally, the oscillator strengths for the transitions between the ground state and the target SO-coupled state f ̃are computed as

In a similar fashion, the XPS intensities can be estimated from the squared norm of the Dyson orbitals, defined as overlap of the initial N-electron wave function of the neutral and the final (N -1)-electron wave function of the cation. For the XPS of ground-state species, the initial wave function is Ψ 0 and the final (SO-split) wave function is Ψ f . Since in CC theory the bra and ket states are not Hermitian conjugates of each other, we use an additional index L or R to mark left and right wave functions and the respective Dyson orbitals 24,[61][62][63]

where

are the expansion coefficients (or amplitudes) of the left (L) and right (R) Dyson orbitals corresponding to the SO-mixed states on the molecular orbital basis {ϕ p }. The target SO-split state Ψ f ̃is the linear combination of the nonrelativistic cationic states

so that the left and right Dyson orbitals for the SO-split target state

where we introduced the amplitudes of the Dyson orbitals of the original nonrelativistic EOM states 24,61 a a

In the last equality of eqs 20 and 21, the Dyson orbitals are expressed on the atomic orbital basis {χ ν }. The Dyson orbitals of the SO-split states are complex-valued. Complex orbitals can be visualized in different ways, from simply representing their real and imaginary components separately to more sophisticated representations that take their phases into account. 64,65 To visualize the complex Dyson orbitals, we use here the QSimulate-QM program, 66 which implements the algorithm proposed in ref 65.

Figure 2 shows zero-order and SO-coupled right Dyson orbitals for H 2 S, illustrating the effect of the SOC on the ionized states. SOC mixes the nonrelativistic states and changes orbital shapes (i.e., 2p 1/2 and 2p 3/2 orbitals are rotated relative to the original p x , p y , and p z ) and also scrambles spin and space degrees of freedom. The SO-mixed orbitals transform by a different symmetry group (double point group), because the relativistic treatment necessitates using different symmetry groups, as described, for example, in refs 67 and 68. For a C 2v molecule (such as H 2 S), the relativistic states belong to the C̅ 2v double group, which is a non-Abelian group with four one-dimensional irreps of the bosonic type and one two-dimensional irrep of the fermionic type. Wave functions with odd and even numbers of electrons transform according to the fermionic and bosonic representations, respectively. The Dyson orbital spinors (Figure 2) have the representation E × A 1 = E. This is a two-dimensional irrep, which corresponds to a Kramers doublet in the full symmetry group (with timereversal operation). Thus, the apparent shape of the orbital depends on how the basis is selected in this irrep. In the Supporting Information, we present an alternative, symmetrized, rendering of the SO-mixed Dyson orbitals, together with the transformation matrix used to generate it and the unconstrained transformation matrix U.

These SO-mixed Dyson orbitals illustrate the effect of the SOC on the ionized states, in the same fashion as the analysis of the SO-mixed transition density matrices from ref 69 illustrates the effect of SOC on the excited states. We note that this analysis is based on the SO-mixed adiabatic states, in contrast to the SOC NTO analysis presented in ref 56, which is formulated in terms of the nonrelativistic (diabatic) states.

The mixing also affects transition strengths, in the same fashion as it affects oscillator strengths. For example, the relative XPS intensities for the transitions involving the SOsplit states are approximated as

The SOC-CVS-EOM-CCSD approach has been implemented in the Q-Chem electronic structure package. 70,71 The implementation included the extension of the fc-CVS-EOM-CCSD framework to the SF states, in addition to the previously implemented 23,24 IP and EE variants. For the XPS calculations, we used the 6-311+G(3df) basis set with uncontracted core functions, denoted as uC-6-311+G(3df), following the recommendation of a recent benchmark study. 6 Because we focus on L-edges, we only uncontracted the core orbitals, roughly corresponding to n = 2, i.e., the second contracted s-function (leaving the "6" contracted core function untouched) and the two most contracted p-functions. For the XAS calculations, we used uC-6-311(2+,+)G(p,d) augmented with additional Rydberg-type functions whose exponents were generated according to the prescription of Kaufmann et al., 72 and quantum number n = 2.5, ..., 5. Uncontracted bases were used for the active edge only, whereas all other atoms were described by the standard variants of these basis sets. All basis sets are given in the Supporting Information. The number of states included in the SOMF Hamiltonian varies depending on the system, the exact number of states for each system can be found in Supporting Information.

To illustrate the capabilities of the method, we considered several systems for which experimental data are available.

The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter These systems are listed in the Supporting Information, along with their structural parameters. In most cases, we used structures optimized at the CCSD(T)/cc-pCVQZ level of theory, taken from ref 73. For thiophene, we used an MP2/cc-pVTZ optimized structure. All Cartesian coordinates are given in the Supporting Information. Table 1 shows the first three core-ionization energies (IEs) for H 2 S, OCS, SO 2 , CS 2 , and C 4 H4S 2 computed with SOC-CVS-EOMIP-CCSD/uC-6-311+G(3df) and compares them with the experimental values. The table also shows the energy difference with respect to the first IE (ΔE). The corresponding zeroth-order energies, calculated at the nonrelativistic fc-CVS-EOMIP-CCSD/uC-6-311+G(3df) level of theory, are given in Table S3 in the Supporting Information. The energies are assigned to the ionization of the electrons from the 2p orbitals of the atom marked in bold. The results show how spin-orbit coupling splits the nearly degenerate 2p orbitals into two sets: one 2p 1/2 orbital and two near-degenerate 2p 3/2 orbitals, as explained in Figure 1 and illustrated in Figure 2. The gap between these two sets is due to the SOC, whereas the small energy difference between the 2p 3/2 orbitals arises from nonspherically symmetric molecular environment (this splitting is called molecular field splitting). Depending on the system, the absolute deviations of the computed IEs relative to the experimental values are on the order of 0.1-1.0 eV, which corresponds to relative errors of the order of 0.05-0.5%. Figure 3 shows the computed XPS spectra of the thiophene molecule, illustrating the spectroscopic signatures of the SOC. At the nonrelativistic level, the 2p orbitals, although already slightly nondegenerate due to the molecular field, are still close enough so that the spectrum has only one peak. The inclusion of the SOC splits this peak into two, with the intensity ratio 2:1, corresponding to the ionization of two 2p 3/2 and the one 2p 1/2 . After a shift of +0.5 eV the SO-corrected spectrum agrees well with the experiment, in terms of both the intensity ratio and the energy splitting.

Figure 4 compares the computed L-edge NEXAFS for SiH 4 , SO 2 , C 4 H 4 S, and Ar with the experimental spectra. The raw theoretical data (energies and oscillator strengths) are provided in the Supporting Information.

Figure 4a shows the L-edge spectra of silane with and without SOCs. As in the XPS example above, the first peak is split into two upon inclusion of the SOC and agrees well with the experiment. The shifted IEs are also in good agreement with the experiment. However, our convoluted spectrum does not reproduce the highly structured set of bands above 104.4 eV observed in the experiment, possibly due to using the same empirical Gaussian broadening function for all computed states, regardless of their actual lifetime.

Figure 4b shows the spectra for the argon atom. Theory and experiment agree well, after a small shift of +0.7 eV is applied. The first band at around 244.5 eV is due to the 2p 3/2 → 4s transition, whereas the second band at around 246.5 eV corresponds to the 2p 1/2 → 4s transition. The third and fourth bands contain contributions from the 2p 3/2 → 5s,3d and 2p 1/2 → 5s,3d transitions, respectively. In the case of sulfur dioxide (see Figure 4c) the agreement between theory and experiment is also quite good for the first peaks (zoomed-in region), after a shift of -0.6 eV is applied. The agreement deteriorates slightly at higher energies.

Finally, Figure 4d illustrates the spectra for a larger molecule, thiophene. In this case, the theoretical spectrum is in excellent agreement with the experimental one in the entire energy range and without an energy shift.

In conclusion, we have presented the first implementation of the perturbative inclusion of spin-orbit effects within coupledcluster theory to describe core-level spectroscopy. This has been achieved by utilizing a general framework for calculating the SOCs from spinless one-particle density matrices computed for the fc-CVS-EOM-CCSD wave functions. This methodological advance enables the calculation of SOcorrected ionization and excitation energies by a simple twostep procedure. In the first step, the nonrelativistic states are computed using appropriate variants of the EOM-CC family of methods; the choice of the method is determined by the target states, i.e., EOM-IP for ionized states and EOM-EE/SF for excited states. In the second step, these zeroth-order states are mixed by the perturbation due to the SO part of the Breit-Pauli Hamiltonian, giving rise to the SO-corrected energies and intensities. The examples illustrate the capabilities of the new method to accurately and efficiently simulate L-edge XAS and XPS.

* sı Supporting Information State University for his insightful comments on relativistic symmetry. S.C. thanks Professor Lan Cheng from the Johns Hopkins University for valuable discussions during the initial phase of this project. M.L.V. and S.C. acknowledge financial support from DTU Chemistry (Ph.D. grant to M.L.V) and from the Independent Research Fund DenmarkNatural Sciences (Research-Project-2 grant No. 7014-00258B to S.C.). The work at USC was supported by the U.S. National Science Foundation (No. CHE-1856342 to A.I.K.).