We present a formulation and implementation of anisotropic and isotropic electronic circular dichroism (ECD) using the full semi-classical light–matter interaction operator within a four-component relativistic framework. Our treatment uniquely accounts for both beyond-first-order light–matter interactions and relativistic effects, enabling us to investigate the ECD response across the electromagnetic spectrum from optical to x-ray wavelengths where relativistic selection rules and spatial field variations gain increasing importance. We consider the isotropic and oriented ECD across the valence transition and sulfur L- and K-edge transitions in the simplest disulfides, H2S2 and (CH3S)2, and evaluate the influence of the full interaction by comparing to a traditional truncated formulation in the Coulomb gauge (velocity representation). Additionally, we demonstrate that in the relativistic formalism, it is possible to work in the velocity representation, hence keeping order-by-order gauge-origin invariance, contrary to the multipolar gauge, yet being able to distinguish electric and magnetic multipole contributions. Going beyond a first-order treatment in the wave vector is mandatory in the higher-energy end of the soft x-ray region and beyond where the consequent intensity redistribution becomes significant. While the sulfur K-edge absorption spectrum is essentially unaffected by this redistribution, the signed differential counterpart is not: At least third-order contributions are required to describe the differential absorption profile that is otherwise overestimated by a factor of about two. The first-order description deteriorates at higher transition energies (beyond ∼1000 eV) where it may even fail to predict the sign of individual differential oscillator strengths.

Chirality is fundamental to life on Earth, and the origin of homochirality of natural amino acids and sugars remains an unresolved question.1 Enantiomers of a chiral molecule interact differently with chiral objects, such as chiral receptors or left and right circularly polarized light. The former is central in chiral molecular recognition, e.g., drug-receptor binding, while the latter enables spectroscopic detection of chirality. Electronic circular dichroism (ECD) is the lowest-order natural chiroptical response: It measures the differential linear absorption between left and right circularly polarized light in the spectral region of electronic transitions.2 Many biomolecules are optically active either intrinsically due to the presence of a chiral center or induced through structural and electronic perturbations by a chiral environment such as through electronic coupling within a chiral arrangement of achiral chromophores.3,4 ECD in the valence region is therefore extensively used to provide stereochemical, conformational, and binding information of molecular systems.5 Anisotropic circular dichroism allows for additional information, which would otherwise be lost in the isotropic signal. As an example, anisotropic ECD was recently proposed as a way to distinguish different transitions in light-harvesting complex II oriented by its embedding in a membrane.6,7

While the discovery of the natural optical activity in the visible region dates back more than two hundred years to the work of Arago8 and Biot,9 its existence in the x-ray range was only established in 1998: The detection of x-ray natural circular dichroism (XNCD) in non-centrosymmetric crystals10–13 was enabled by the advent of intense third-generation synchrotron x-ray sources with full polarization control. Like other resonant core-level spectroscopies, XNCD exploits the site- and element-specificity of x-ray probes to provide information about the local symmetry around the absorbing atom(s). In addition to chiral crystals, XNCD has also been measured for small organic molecules including several amino acids (see, e.g., Refs. 14–17).

From a multipole expansion of the light–matter interactions in orders of the wave vector, the leading-order contributions to ECD (i.e., first order) are governed by the electric-dipole/magnetic-dipole (E1–M1) and electric-dipole/electric-quadrupole (E1–E2) interference terms. The latter is traceless and hence vanishes for isotropic samples. While the electric-dipole/magnetic-dipole interference dominates in the optical region (in the center-of-mass frame), it is difficult to observe at K(L1)-edges18,19 because the magnetic transition-dipole moment involving the atomic-like 1(2)s-orbitals becomes forbidden within a zeroth-order approximation in the non-relativistic domain.13,20 The electric-dipole/electric-quadrupole interference term, on the other hand, gains intensity in the x-ray region due to its dependency on the spatial variations of the electromagnetic field.20,21 For these reasons, substantial XNCD signals are expected only for oriented samples, although weak electric-dipole/magnetic-dipole contributions have been reported in the absence of orientational order (e.g., in solution or powdered samples19,22).

Theoretical simulations play a key role in the assignment of experimental ECD spectra. Traditionally, these rely on the aforementioned multipole expansion of the ECD signal through first order in the wave vector. However, caution has to be exercised to ensure gauge-origin independent results, e.g., by performing calculations in the velocity representation23 or by employing gauge-including atomic orbitals.24,25 An additional aspect arising in the x-ray regime, where the wavelength of the electromagnetic field approaches molecular dimensions, is the importance of higher-order terms in the expansion. To address these issues in the context of linear absorption of linearly polarized light, we recently proposed using the full semi-classical light–matter interaction operator in both the non-relativistic26,27 and relativistic regimes.28 In addition to providing compatible transformation properties for both light and matter, the relativistic formulation is particularly simple because the light–matter interaction is linear in both the scalar and vector potentials.

In this work, we present a four-component relativistic formulation and implementation of isotropic and anisotropic ECD using the full semi-classical light–matter interaction. This allows us to investigate the ECD response across the electromagnetic spectrum, from optical to x-ray regimes. So far, previous theoretical studies of the isotropic XNCD of molecules have only considered the first-order truncated interaction in a non-relativistic framework.18,22,29–35 In these studies, the only source of non-vanishing response at the L1- and K-edges is the polarization away from the atomic symmetry of the core orbitals. Our implementation accounts for two additional possible contributions: (i) effects of beyond first-order light–matter interactions and (ii) inclusion of relativistic effects, notably spin–orbit coupling that modifies selection rules (in particular, the magnetic transition-dipole selection rule36,37). Consequently, this allows, for the first time, to realistically examine the ECD response of molecules across the valence and x-ray regimes.

As test systems, we consider the simplest disulfide chromophore models, dihydrogen disulfide H2S2 and dimethyl disulfide (CH3S)2. Because of the low disulfide torsional barriers (∼6–11 kcal/mol38–41), the two enantiomeric forms (P- and M-helix) cannot be resolved experimentally. However, the disulfide bridge is an important structural element in proteins,42 where it preferentially occurs in non-planar, chiral conformations (C2 symmetry) and hence displays structurally induced axial chirality. An interesting perspective for complex systems (e.g., proteins) is the potential use of XNCD as a local probe of chirality.30,35 This could potentially complement the delocalized conformational information encoded in valence ECD. Because of its computational tractability, H2S2 has been widely used to benchmark electronic structure methods for the calculation of chiroptical properties.43–49 For the same reason, Goulon et al. also used it to estimate the relative magnitudes of XNCD responses within the first-order truncated interaction and non-relativistic framework, reporting values below the experimental detection limits.18 Here, we revisit the ECD of the disulfide chromophore across the valence and L- and K-edges, going beyond these approximations.

We begin with a description of our ECD formulation using the full semi-classical light–matter interaction, which has been implemented in a development version of the Dirac electronic structure program.50,51 We consider a circularly polarized monochromatic (angular frequency ω) electromagnetic field with electric and magnetic components given by

E(M)(r,t)=12Eωϵ1sin(krωt)+(1)Mϵ2cos(krωt),B(M)(r,t)=12Eωcϵ2sin(krωt)(1)Mϵ1cos(krωt),
(1)

where appears the polarization vectors ϵ1 and ϵ2 that, together with the wave vector k = kek(k=ωc), form a right-handed coordinate system; the field amplitude Eω; and the integer M that specifies the handedness of the field, i.e., left (even) and right (odd) circularly polarized light, respectively, following the International Union of Pure and Applied Chemistry (IUPAC) recommendation.52 Such an electromagnetic wave is conventionally represented in the Coulomb gauge by the scalar and vector potentials,

ϕ(M)(r,t)=0,A(M)(r,t)=12Eωωϵ1coskrωt(1)Mϵ2sinkrωt.
(2)

In the relativistic domain, this leads to an interaction operator for left/right-handed (L/R) circularly polarized light of the form

V̂L/R(t)=ecαAL/R=12Eω±ωT̂L/R(ω)eiωt,T̂L/R(ω)=eω(cαϵL/R)eikr,
(3)

where T̂L/R(ω) is the effective interaction operator (cf. Ref. 28), α are Dirac matrices, and

ϵL/R=12ϵ1±iϵ2=ϵR/L*.
(4)

For notational simplicity, we tacitly assume the summation over electrons in the effective interaction operator in Eq. (3). Considerable simplifications arise from the fact that the relativistic interaction operator is linear in the vector potential and there is no separation into spin and spatial parts.

The linear absorption cross section between initial and final states, |i⟩ and |f⟩, respectively, is then given by

σL/Rω=πωε0cϵL/RT2fω,ωfi,γfi,T=f|eωcαeikr|i.
(5)

The use of Fermi’s golden rule for a transition between discrete states is justified by the inclusion of a (Lorentzian) line shape function f(ω, ωfi, γfi) expressing the finite lifetime of excited states (see the discussion in Ref. 28). Moreover, the line shape function, having the dimension of time, has to be included in order for the absorption cross section to have the correct dimension of area. On the other hand, the dimensionless differential oscillator strength is obtained from the absorption cross section by the pre-factor substitution πωϵ0c2meωe2 and removal of the line shape function. In the interest of compact expressions, we shall in the following carry out our derivations using oscillator strengths rather than absorption cross sections.

Evaluating the differential oscillator strengths between left and right circularly polarized light starting from Eq. (5) provides the expression for the anisotropic differential oscillator strength using the full semi-classical interaction operator,

Δf=fLfR=2meωe2ϵLTϵRT*ϵLT*ϵRT=2meωe2ϵL×ϵRT×T*=i2meωe2ekT×T*,
(6)

where we have used Eq. (4) and the vector relation

a×bc×d=acbdadbc.
(7)

This final, compact expression represents the relativistic extension of the expression previously reported by Hansen and Avery.53 

It is instructive to compare the differential oscillator strength in Eq. (6) (ECD) with the oscillator strength associated with linearly polarized light. The latter is obtained by replacing the complex polarization vector ϵL/R in Eq. (5) by the polarization vector ϵ giving the direction of the electric field vector in the plane perpendicular to the wave vector. One immediately notes that the differential oscillator strength contains no reference to polarization vectors at all, which can be understood as a manifestation of axial symmetry. This simplifies the rotational averaging appropriate for samples such as liquids and gases, with freely tumbling molecules. Specifically, the rotational average involves three angles (θ, ϕ, χ): The first two are associated with an average over propagation direction, whereas the polarization angle χ now becomes redundant.

1. Multipoles in the velocity representation

In our previous work on the absorption of linearly polarized light, we considered two schemes involving truncated light–matter interaction.28 In the first scheme, referred to as the generalized velocity representation, the truncated interaction is obtained via a Taylor expansion of the oscillator strength in orders of the wave vector. This is a relativistic extension of the theory developed by Bernadotte et al.54 Applying this to Eq. (6), we obtain the following expressions for the differential oscillator strength to even and odd order in the wave vector:

Δf2n=2meωe2ekm=0n1m2δm0×ImTn+m×Tnm*=0,
(8)
Δf2n+1=2meωe2ekm=0n1m2ReTn+m+1×Tnm*,
(9)

where

TnϵL/R=f|T̂L/Rn|i,T̂L/Rn=eω1n!icαϵL/Rkrn.
(10)

Note that only odd-order terms are sensitive to the handedness and hence contribute to the expansion. This is a consequence of the time-reversal symmetry of the truncated interaction operators, making the corresponding transition moments real.28 It should be stressed that the (differential) oscillator strength is gauge-origin invariant in each order, whereas the transition moments are intrinsically gauge-origin dependent.28,54

In the second scheme, referred to as generalized length representation, we transform to multipolar gauge (mg).28,55–59 Using the Einstein summation convention, the full interaction operator is then expressed as

V̂mgω=n=0V̂mgnω,V̂mg0ω=Q̂p1Ep[0],V̂mgn0ω=1n+1!Q̂j1jnpn+1Ep;j1jn[n]1n!m̂j1jn1;rnBr;j1jn1[n1],
(11)

where appears electric and magnetic multipole operators,

Q̂j1jnpn+1=erj1rj2rjnrp,
(12)
m̂j1jn1;rn=nn+1rj1rj2rjn1(r×ĵ)r,ĵ=ecα,
(13)

which couple to the electric and magnetic fields and their derivatives,

Fp;j1jn[n]=nFprj1rj2rjnr=0,F=E,B,
(14)

at a selected expansion point, here set to the origin. The corresponding nth-order effective interaction operator is given by

T̂mg;L/R[n]ω=e1n+1!rϵL/Rikrniωnn+1!×ik×ϵL/Rr×cαikrn1.
(15)

In the component form, the effective interaction operator is conveniently expressed as

T̂mg;L/R[n]ω=inϵL/Rpkj1kj2kjnX̂j1jn;p[n]ω,
(16)

in terms of a multipole operator,28 

X̂j1jn;p[n]ω=1n+1!Q̂j1jn,pn+1iω1n!m̂j1jn1;rnεrjnp,
(17)

(εrjnp is the Levi–Cività symbol). However, in contrast to the generalized velocity representation, we now loose the origin independence of oscillator strengths.28 

It is possible, though, to combine the best of both worlds. In our previous work,28 we have shown that it is possible within the generalized velocity representation to distinguish electric and magnetic multipole components, while maintaining origin independence of the oscillator strengths, although its individual components still exhibit this dependency. The operator manipulations required for this separation do not involve commutators with the Hamiltonian and are therefore also valid when using finite basis sets. We here review and expand this result. Using the vector relation in Eq. (7), the effective interaction operator of Eq. (10) can be rewritten as

T̂L/Rn=eω1n+1!icαϵL/Rkrn+nkicα×ϵL/Rrkrn1+nn+1!k×ϵL/Rr×icαkrn1.
(18)

As evident upon comparison to Eq. (15), we have effectively recovered the magnetic multipole part of the operator, and, in fact, writing the corresponding nth-order interaction operator V̂L/Rn out in the component form, we recover the expression of Eq. (11), with the important modification that the electric multipole operator in the length representation Eq. (12) is replaced by its counterpart in the velocity representation

Q̂j1jnpn+1=ieωrj1rjn1cαprjn+ncαjnrp.
(19)

This somewhat curious form is perhaps best understood by noting that the operator can equally well be expressed as

Qj1jnpn+1=ieωcαj1rj2rjnrp++rj1rj2cαjnrp+rj1rj2rjncαpn+1 terms,
(20)

being perfectly symmetric in all indices j1jnp, just as Eq. (12). The more compact expression in Eq. (19) is possible because the relativistic velocity operator cα commutes with the coordinates, contrary to the corresponding non-relativistic operator p/m. Finally, using the commutator relation with the Hamiltonian

irϵL/Rkrn,ĥ=cαϵL/Rkrn+nrϵL/R×kcαkrn1,
(21)

one can show that transition moments of the electric multipole operator in the length and velocity representation will be identical.28 However, very importantly, this generally only holds in a complete basis.

The odd-order oscillator strength of Eq. (9) may now be expressed as

Δf2n+1=2meωe2k2n+1ek;j1ek;j2n+1ekm=0n1m2Xj1jn+m+1n+m+1×Xjn+m+2j2n+1[nm],
(22)

where the components of the real vector Xj1jn[n] are given as the transition moments of the multipole operator in Eq. (17) in the velocity representation, that is, using Eq. (19) for the electric multipole operator.

2. Angular dependence of truncated ECD

In the anisotropic case, both the full and truncated differential oscillator strengths [Eqs. (6) and (22), respectively] are functions of the orientation of the incident radiation. In the case of truncated interaction, we make this angular dependence explicit by writing the unit wave vector as

ek=sinθcosϕ,sinθsinϕ,cosθ.
(23)

The first-order differential oscillator strength can then be expanded as

Δf[1]θ,ϕ=ek;iRij[1]ek;j=Riso[1]+12Rzz[1]Riso[1]3cos2θ1+12Rxz[1]+Rzx[1]sin2θcosϕ+12Ryz[1]+Rzy[1]sin2θsinϕ+12Rxx[1]Ryy[1]sin2θcos2ϕ+12Rxy[1]+Ryx[1]sin2θsin2ϕ,
(24)

where we have introduced the rotational strength tensor Rij[1] (cf. Ref. 60) and its isotropic part Riso[1], known as the rotational strength,

Rij[1]=4meω2ce2Xj1×X[0]i,Riso[1]=13Rxx[1]+Ryy[1]+Rzz[1].
(25)

From Eq. (24), we see that the anisotropic part of Δf[1](θ, ϕ) is a linear combination of d-orbitals weighted by the relevant elements of the rotatory strength tensor. This can be made more explicit by rewriting the first-order differential oscillator strength as

Δf[1](θ,ϕ)=15Riso[1]s+3Rzz[1]Riso[1]dz2+Rxx[1]Ryy[1]dx2y2+Rxz[1]+Rzx[1]dxz+Ryz[1]+Rzy[1]dyz+Rxy[1]+Ryx[1]dxy,
(26)

where we have used a common normalization for all solid harmonics.

From Eq. (22), we see that rotational strength tensors can be generalized to arbitrary odd orders. It will also be seen that the differential oscillator strength Δf2n+1 contains products of 2(n + 1) components of the unit wave vector and is therefore spanned by spherical harmonics of even = 0, 2, …, 2(n + 1).

In the supplementary material, the (first-order) rotational strength tensor is analyzed in terms of symmetry for the C2 point group of the disulfide configurations considered, assuming the rotation axis to coincide with the z-axis. We also show that the rotational strength tensor in the generalized velocity representation is invariant under a shift of origin, also in a finite basis, which generally does not hold true in the generalized length representation. As pointed out above, the separation into electric and magnetic multipole contributions generally depends on the choice of gauge origin and, thus, is not unique. However, the presence of symmetry may introduce additional gauge-origin invariances.61 The E1–M1 and E1–E2 contributions to the Rz*[1] elements are separately invariant under origin shifts along the rotation axes (Sec. S3); in A symmetry, this holds for all elements of the full rotational strength tensor.

In the case of the full interaction, the factorization of radial and angular parts is possible using the plane wave expansion

eikr=4π=0m=ij(kr)Ym(ek)Ym*(er)
(27)

in terms of spherical harmonics Ym and spherical Bessel functions j. However, it is more cumbersome since it involves infinite sums over angular momentum and will not be pursued further here.

As pointed out above, rotational averaging is simplified for the differential oscillator strength because of the absence of reference to the polarization vectors. It is therefore limited to the solid angle indicating the propagation direction and is defined as

g(r)θ,ϕ=14π02π0πg(r)sinθdθdϕ.
(28)

In the case of the full light–matter interaction, the rotational average is handled numerically using a Lebedev grid.27,28 In the case of the truncated interaction, the rotational average is expressed as

Δf2n+1θ,ϕ=Δfiso2n+1=2meωe2k2n+1ek;j1ek;j2n+1ekθ,ϕm=0n1m2Xj1jn+m+1n+m+1×Xjn+m+2j2n[nm].
(29)

In the lowest order, the isotropic differential oscillator strength becomes

Δf1θ,ϕ=4meω23ce2jXj1×X[0]j=8meω23ce2jf|Q̂j1|if|iωm̂j1|i=Riso[1],
(30)

where we have used that ek;jekθ,ϕ=13ej and that the E1–E2 contribution vanishes since its elements form a traceless matrix.

In practice, the interaction operator is separated into Hermitian and anti-Hermitian parts for compatibility with the quaternion symmetry scheme of Dirac (see Refs. 28 and 62). The implementation of Eq. (6) requires the same integrals as the linear oscillator strength, and thus, our previous integral implementation could be readily extended to ECD.26,27 The integral evaluation over Lebedev points has been parallelized using MPI (Message Passing Interface).

The geometries of H2S2 and (CH3S)2 were obtained using the B3LYP63–66 exchange–correlation functional and the cc-pVTZ67,68 basis set. Geometry optimizations were performed in Gaussian 16.69 To mimic the χ3 disulfide angle typical for protein structures,70 we performed constrained geometry optimizations for χ3 = −87°, corresponding to M-helical chirality.71 The restricted excitation window approach72,73 was used to selectively target the sulfur L- and K-edges. This also eliminates the issue of artificial transitions to quasi-continuum orbitals caused by finite basis set effects that otherwise often interfere with simulations at the L-edge.74,75 A Gaussian model was employed for the nuclear charge distribution,76 and an 86-point Lebedev grid (Lmax = 12) was used for the isotropic averaging of the differential linear absorption based on the full interaction operator. The gauge origin was placed in the center-of-mass (COM), and spatial symmetry was invoked in all cases except for the gauge-origin dependence calculations.

Excitation energies, linear and differential absorption cross sections for the full interaction operator, and the multipole expansions within the generalized velocity gauge were computed using the PBE077,78 exchange–correlation functional and the uncontracted aug-pcX-379 and aug-pc-380–82 basis sets for sulfur and hydrogen, respectively. The pcX-n basis set series was developed for describing core-excitation processes using the ΔSCF (Self-Consistent Field) approach at both the nonrelativistic and relativistic levels. The small component basis sets were generated within the condition of restricted kinetic balance. The relativistic calculations were performed using a Dirac–Coulomb Hamiltonian in which the (SS|SS) integrals are replaced by an interatomic SS energy correction.83 The gauge-origin invariance of our implementation of the full semi-classical formulation of the isotropic and anisotropic rotatory strengths [Eq. (6)] and its first-order truncated counterpart was confirmed numerically by shifting the gauge-origin (from 0 to 100 a0) along the C2 axis. This leads to a redistribution of the E1–E2 and E1–M1 contributions to Rxx[1] and Ryy[1] for transitions of B symmetry. As expected, the results remained unchanged for both the full and truncated formulations (data not shown). Simulated spectra were obtained by convolving the stick spectrum with Gaussian line shape functions with a full width at half maximum (FWHM) of 0.4 eV, and those for (CH3S)2 were shifted by different offsets for each absorption edge to match their experimental counterparts.

Before considering the ECD response of the disulfides, we assign the linear absorption features across the valence and x-ray regions. We initially focus on (CH3S)2 for which experimental gas-phase absorption spectra are available.84–88 Expectedly, and as shown in Fig. S1, the spectral profiles for H2S2 are similar, and because of its greater computational tractability, we consider this minimal disulfide in subsequent analyses.

Figure 1 displays the rotationally averaged linear and differential absorption spectra for valence transition and sulfur L- and K-edge transitions of (CH3S)2, computed using the full interaction operator (green shading) with corresponding oscillator strengths indicated as green sticks. Hereafter, we explicitly indicate the results of the full interaction with the superscript “full.” For comparison, we also provide the lowest non-vanishing terms in the truncated generalized velocity representation (orange lines) for linear and differential absorption, i.e., zeroth- and first-order in the magnitude of the wave vector, respectively. The black sticks at the top indicate the location of the underlying electronic transitions, and black lines (solid and dashed) indicate the experimental absorption spectra.84–88 Apart from uniform shifts necessary to align the lowest-energy band to the respective experimental spectrum, the theoretical spectra capture well both relative intensities and peak splittings.

FIG. 1.

Isotropic linear (top) and differential (bottom) absorption spectra of (CH3S)2: (a) valence, (b) L2,3-edge, (c) L1-edge, and (d) K-edge spectra using the full interaction operator in Eq. (3) (green shadings) or the lowest non-vanishing generalized velocity representation (orange lines). The left axes correspond to (differential) absorption cross sections, whereas (differential) oscillator strengths (sticks) are shown on the right axes. The black sticks indicate the location of all computed transitions, whereas the black lines are experimental spectra.84–88 The stick spectra were convolved with a Gaussian line shape with a FWHM of 0.4 eV. The theoretical absorption spectra have been uniformly shifted to align with the experimental counterparts (shift values indicated in the top panels). The same shifts were applied to the ECD spectra.

FIG. 1.

Isotropic linear (top) and differential (bottom) absorption spectra of (CH3S)2: (a) valence, (b) L2,3-edge, (c) L1-edge, and (d) K-edge spectra using the full interaction operator in Eq. (3) (green shadings) or the lowest non-vanishing generalized velocity representation (orange lines). The left axes correspond to (differential) absorption cross sections, whereas (differential) oscillator strengths (sticks) are shown on the right axes. The black sticks indicate the location of all computed transitions, whereas the black lines are experimental spectra.84–88 The stick spectra were convolved with a Gaussian line shape with a FWHM of 0.4 eV. The theoretical absorption spectra have been uniformly shifted to align with the experimental counterparts (shift values indicated in the top panels). The same shifts were applied to the ECD spectra.

Close modal

The first valence band [Fig. 1(a)] is dominated by the two excitations from the symmetric and antisymmetric combinations of the non-bonding 3p orbital on each sulfur to the lowest unoccupied σSS* orbital (b symmetry). This assignment is consistent with the analysis by Linderberg and Michl on H2S2.89 Since the valence orbitals are found to have well-defined spin, we adopt a non-relativistic state notation for the valence states, i.e., 11B and 21A, respectively. They are separated by ∼0.15 eV. The second valence band originates from transitions into the σCS* orbital (31A and 21B).

Turning to the x-ray region, the first two bands at the sulfur L2,3-edge [Fig. 1(b)] are dominated by transitions from the symmetric and antisymmetric combinations of the S 2p3/2 and 2p1/2 core orbitals to the σSS* orbital. We note that the A/B pair of associated transitions are essentially degenerate (less than 0.1 meV splitting) because of the limited overlap between the core orbitals on each of the sulfur atoms. The calculated spin–orbit splitting of ∼1.35 eV between the L3- and L2-branches is comparable to the experimental splitting reported for dihydrogen sulfide (∼1.2 eV90,91). The L3/L2 branching ratio of ∼1.4:1 (obtained by summing the underlying oscillator strengths) deviates significantly from its statistical value, which is obtained only in the limit of jj coupling.92,93 The second peak in each branch is dominated by excitations to the σCS* orbital and is separated from the first peak by ∼1.4–1.5 eV. Consequently, the second band in the spectrum contains contributions from both branches, whereas the third band is associated with the second peak in the L2 branch. These assignments agree with previous studies.86,94 The energy range considered as well as the basis set used in our calculations does not cover the fourth band in the experimental spectrum that, according to previous work,86 originates from excitations to higher-lying orbitals of mixed σCS*/Rydberg character. Not surprisingly, the L1- and K-edge spectra bear strong resemblance [Figs. 1(c) and 1(d)]: They display two pre-edge features, separated by ∼1.5 eV, which originate from pairs of near-degenerate excitations from the bonding and antibonding combinations of sulfur s-orbitals into the σSS* and σCS* orbitals, respectively.87,95,96

A non-vanishing ECD response in these minimal disulfides results from axial chirality caused by trapping the disulfide bridge in a non-planar (i.e., C2) conformation. As described above, the two lowest-energy transitions in each spectral domain are dominated by an excitation from the bonding or antibonding combinations of the relevant atomic orbitals on the sulfurs into the σSS* orbital. This pairing of transitions manifests as bisignate features in the low-energy region of the ECD spectra. On the basis of the simple Bergson model for the low-energy transitions in the disulfide chromophore,97,98 Linderberg and Michl89 formulated a quadrant rule for the optical activity of the two low-energy valence transitions (dominated by excitations from the symmetric and antisymmetric combinations of non-bonding 3p orbitals on sulfurs to the σSS* orbital) in organic disulfides. This rule relates the sign of the long-wavelength Cotton effect across the four dihedral quadrants and is a specific case of the C2-rule for general chromophores of effective C2 symmetry.99 Woody extended the theoretical analysis to also include the absolute sign of the lowest ECD band within each quadrant,100 providing predictions in agreement with experimental results across different dihedral angles.101–104 

For the M-helical form considered here, we find a negative-first Cotton effect, consistent with the quadrant rule.100 The intensity asymmetry of the lowest-energy valence couplet is attributed to different intrafragment (i.e., CH3S–) contributions to the ECD signal of each transition. At higher energies, the electronic coupling between the (core) orbitals on different fragments decreases, reducing both the energetic splitting between symmetric and antisymmetric combinations as well as the intrafragment ECD contribution itself, since the environment is increasingly achiral on each fragment. As a consequence, the paired core transitions become near-degenerate (energy splitting of a few meV or less) with weak rotational strengths of almost equal magnitudes but opposite signs (see Table I). After additionally accounting for sources of broadening, including finite core-hole lifetimes (∼0.1 and ∼0.5 eV at the sulfur L- and K-edges,105,106 respectively), the bisignate feature may coalesce into a single signal with the absolute sign given by the most intense of the pair of transitions (this will be seen in Fig. S2). In sum, the effective signals are reduced by orders of magnitude in the x-ray region (see the accumulated values in Table I).

TABLE I.

Beyond first-order effects in the isotropically averaged and anisotropic rotatory strengths for the lowest intense valence transition and L- and K-edge transitions of (CH3S)2 for the full semi-classical light–matter interaction operator and in first-order within the Coulomb gauge, computed at the 4c-TD-PBE0 level of theory and the uncontracted aug-pcx-3/aug-pc3 basis set. For the oriented case, θ = ϕ = π/2 for the valence and L2,3-edge transitions and θ = π/2, ϕ = π/4 for the L1- and K-edge transitions [cf. Eq. (23)]. For each spectral region, the accumulated (differential) oscillator strength (ac) across the two pairs of transitions (A/B symmetry) is given in the third line. The numbers in parentheses are exponents of 10.

TransitionΔE (eV)fisofullfiso[0]ΔfisofullΔfiso[1]ΔfkfullΔfk[1]
Valence 
11B(nSσSS*) 4.695 14  7.7977(−03) 7.7974(−03) −1.0856(−04) −1.0856(−04) −6.1154(−05) −6.1147(−05) 
21A(nSσSS*) 4.770 69  2.0797(−03) 2.0793(−03) 8.1247(−05) 8.1244(−05) 1.1759(−04) 1.1758(−04) 
 ac 9.8773(−03) 9.8767(−03) −2.7316(−05) −2.7317(−05) 5.6436(−05) 5.6435(−05) 
L2,3-edge 
6B(2p3/2σSS*) 158.087 21  3.1232(−03) 3.1230(−03) 5.0231(−04) 5.0326(−04) −1.1661(−03) −1.169 2(−03) 
6A(2p3/2σSS*) 158.087 23  2.6587(−03) 2.6566(−03) −4.9798(−04) −4.9892(−04) 1.1726(−03) 1.1757(−03) 
 ac 5.7818(−03) 5.7796(−03) 4.3272(−06) 4.3363(−06) 6.4960(−06) 6.5232(−06) 
10A(2p1/2σSS*) 159.293 75  1.2239(−03) 1.2227(−03) −2.5299(−04) −2.5345(−04) −5.7854(−04) −5.7988(−04) 
10B(2p1/2σSS*) 159.293 87  1.3957(−03) 1.3953(−03) 2.5605(−04) 2.5651(−04) 5.8340(−04) 5.8476(−04) 
 ac 2.6196(−03) 2.6180(−03) 3.0568(−06) 3.0622(−06) 4.8616(−06) 4.8783(−06) 
L1-edge 
4B(2s1/2σSS*) 216.393 99  3.6733(−02) 3.6872(−02) −3.1538(−05) −3.1716(−05) −1.0027(−03) −1.008 5(-03) 
4A(2s1/2σSS*) 216.395 99  2.9891(−04) 1.7881(−04) 3.1469(−05) 3.1646(−05) 1.0279(−03) 1.033 8(-03) 
 ac 3.7032(−02) 3.7051(−02) −6.9310(−08) −6.9637(−08) 2.5196(−05) 2.521 2(-05) 
K-edge 
4B(1s1/2σSS*) 2427.833 34  7.7555(−03) 1.0438(−02) −2.2670(−05) −4.6773(−05) −1.0771(−03) −1.984 2(-03) 
4A(1s1/2σSS*) 2427.833 49  2.6746(−03) 2.4290(−05) 2.2665(−05) 4.6761(−05) 1.0754(−03) 1.982 3(-03) 
 ac 1.0430(−02) 1.0462(−02) −5.6000(−09) −1.1871(−08) −1.6290(−06) −1.821 6(-06) 
TransitionΔE (eV)fisofullfiso[0]ΔfisofullΔfiso[1]ΔfkfullΔfk[1]
Valence 
11B(nSσSS*) 4.695 14  7.7977(−03) 7.7974(−03) −1.0856(−04) −1.0856(−04) −6.1154(−05) −6.1147(−05) 
21A(nSσSS*) 4.770 69  2.0797(−03) 2.0793(−03) 8.1247(−05) 8.1244(−05) 1.1759(−04) 1.1758(−04) 
 ac 9.8773(−03) 9.8767(−03) −2.7316(−05) −2.7317(−05) 5.6436(−05) 5.6435(−05) 
L2,3-edge 
6B(2p3/2σSS*) 158.087 21  3.1232(−03) 3.1230(−03) 5.0231(−04) 5.0326(−04) −1.1661(−03) −1.169 2(−03) 
6A(2p3/2σSS*) 158.087 23  2.6587(−03) 2.6566(−03) −4.9798(−04) −4.9892(−04) 1.1726(−03) 1.1757(−03) 
 ac 5.7818(−03) 5.7796(−03) 4.3272(−06) 4.3363(−06) 6.4960(−06) 6.5232(−06) 
10A(2p1/2σSS*) 159.293 75  1.2239(−03) 1.2227(−03) −2.5299(−04) −2.5345(−04) −5.7854(−04) −5.7988(−04) 
10B(2p1/2σSS*) 159.293 87  1.3957(−03) 1.3953(−03) 2.5605(−04) 2.5651(−04) 5.8340(−04) 5.8476(−04) 
 ac 2.6196(−03) 2.6180(−03) 3.0568(−06) 3.0622(−06) 4.8616(−06) 4.8783(−06) 
L1-edge 
4B(2s1/2σSS*) 216.393 99  3.6733(−02) 3.6872(−02) −3.1538(−05) −3.1716(−05) −1.0027(−03) −1.008 5(-03) 
4A(2s1/2σSS*) 216.395 99  2.9891(−04) 1.7881(−04) 3.1469(−05) 3.1646(−05) 1.0279(−03) 1.033 8(-03) 
 ac 3.7032(−02) 3.7051(−02) −6.9310(−08) −6.9637(−08) 2.5196(−05) 2.521 2(-05) 
K-edge 
4B(1s1/2σSS*) 2427.833 34  7.7555(−03) 1.0438(−02) −2.2670(−05) −4.6773(−05) −1.0771(−03) −1.984 2(-03) 
4A(1s1/2σSS*) 2427.833 49  2.6746(−03) 2.4290(−05) 2.2665(−05) 4.6761(−05) 1.0754(−03) 1.982 3(-03) 
 ac 1.0430(−02) 1.0462(−02) −5.6000(−09) −1.1871(−08) −1.6290(−06) −1.821 6(-06) 

The linear absorption profiles with the full interaction and the electric-dipole approximation essentially coincide across the four spectral regions. However, as shown by the underlying oscillator strengths in Table I, non-dipolar effects at the K-edge lead to intensity redistribution among the underlying near-degenerate transitions (i.e., unrelated to the arbitrary mixing allowed for degenerate states). Nonetheless, the overall spectral profiles within and beyond the electric-dipole approximation are essentially identical because of the nearly overlapping transitions. This is consistent with our previous findings for Cl K-edge absorption in TiCl4.28 In contrast, the beyond-first-order effects become evident in the differential K-edge absorption profile because of the signed nature of the underlying quantities. This leads to a factor-of-two overestimation of the ECD within the conventional first-order treatment. Introducing third-order contributions largely corrects this discrepancy at the sulfur K-edge, but going to higher orders in the expansion is not a general remedy, as will be discussed below.

To better understand the nature of the chiral response across spectral regions, we computed the underlying anisotropic differential oscillator strength distributions, considering now the smaller H2S2. Figure 2 shows the full ECD distributions (points), compared with the first-order truncated counterparts (surfaces). The solid angle represents the propagation direction, the distance from the origin (COM) indicates the magnitude of the associated signal, and the color indicates its sign. Note that different scaling factors have been applied across the transitions (see the upper right corner of each inset). The C2-rotation axis coincides with the z-axis, whereas the disulfide bond is along the x-axis. The shapes of the anisotropic distributions can be understood by decomposing the first-order signals into isotropic and d-orbital contributions [Eq. (26)]. The resulting orbital weights of the excitations plotted in Fig. 2 are reported in Table II.

FIG. 2.

Comparison of full and truncated differential oscillator strength Δf(θ, ϕ) across the spectral regions (valence, L3-edge, L1-edge, and K-edge). Transitions of [(a)–(d)] A and [(e)–(h)] B symmetry in H2S2. The black arrow points along the direction of the wave vector for the anisotropic ECD intensity given in Table I. The truncated ECD is represented by the smooth surface, whereas the full ECD is shown as individual points generated with a 5810-point Lebedev grid (Lmax = 131). Blue: negative and red: positive ECD signal. Note that different scaling factors (upper right corner) have been applied. The corresponding isotropic differential oscillator strengths are indicated in each inset.

FIG. 2.

Comparison of full and truncated differential oscillator strength Δf(θ, ϕ) across the spectral regions (valence, L3-edge, L1-edge, and K-edge). Transitions of [(a)–(d)] A and [(e)–(h)] B symmetry in H2S2. The black arrow points along the direction of the wave vector for the anisotropic ECD intensity given in Table I. The truncated ECD is represented by the smooth surface, whereas the full ECD is shown as individual points generated with a 5810-point Lebedev grid (Lmax = 131). Blue: negative and red: positive ECD signal. Note that different scaling factors (upper right corner) have been applied. The corresponding isotropic differential oscillator strengths are indicated in each inset.

Close modal
TABLE II.

Weights of contributions from solid harmonics [(24)] to the differential cross section Δf(θ, ϕ) of selected transitions in H2S2. The weights have been scaled by the absolute value of the s-contribution. The numbers in parentheses are exponents of 10.

IrrepExcitationΔfiso[1]sdz2dx2y2dxy
A Valence 5.176(−05) 1.000 −0.447 0.300 −0.056 
L3-edge −9.118(−05) −1.000 0.447 0.414 −0.031 
L2-edge −2.824(−04) −1.000 0.447 0.408 −0.024 
L1-edge 7.316(−05) 1.000 −0.447 −0.862 8.817 
K-edge 1.657(−04) 1.000 −0.447 −0.772 9.108 
B Valence −7.624(−05) −1.000 −0.328 −0.090 0.314 
L3-edge 9.186(−05) 1.000 −0.132 −0.592 −0.005 
L2-edge 2.883(−04) 1.000 −0.128 −0.591 −0.008 
L1-edge −7.312(−05) −1.000 2.503 −0.347 −8.262 
K-edge −1.657(−04) −1.000 0.377 0.813 −9.125 
IrrepExcitationΔfiso[1]sdz2dx2y2dxy
A Valence 5.176(−05) 1.000 −0.447 0.300 −0.056 
L3-edge −9.118(−05) −1.000 0.447 0.414 −0.031 
L2-edge −2.824(−04) −1.000 0.447 0.408 −0.024 
L1-edge 7.316(−05) 1.000 −0.447 −0.862 8.817 
K-edge 1.657(−04) 1.000 −0.447 −0.772 9.108 
B Valence −7.624(−05) −1.000 −0.328 −0.090 0.314 
L3-edge 9.186(−05) 1.000 −0.132 −0.592 −0.005 
L2-edge 2.883(−04) 1.000 −0.128 −0.591 −0.008 
L1-edge −7.312(−05) −1.000 2.503 −0.347 −8.262 
K-edge −1.657(−04) −1.000 0.377 0.813 −9.125 

From symmetry considerations detailed in the supplementary material, we find that contributions from dxz and dyz vanish for excitations of both A and B symmetries. In A symmetry, Rzz[1] is also zero by symmetry, such that the s- and dz2-contributions come in a fixed ratio, giving a toroid in the xy-plane. For valence and L3-edge excitations of A symmetry, this shape is modulated by the dx2y2-contribution, giving the shape of a biconcave disk elongated along the y- and x-axis, respectively, depending on its relative sign. For the L1- and K-edge excitations, the dxy-contribution completely dominates. This is also the case for the corresponding excitations in the B symmetry. It may be noted that in the A symmetry, the E1–M1 contribution to the dxy term is zero by symmetry, whereas in the B symmetry, the E1–E2 and E1–M1 contributions are of similar magnitude and the same sign. Continuing to the L3 excitation of the B symmetry, the angular plot resembles that of its counterpart in the A symmetry, albeit with opposite overall sign and rotated π/2 about the molecular axes. The latter can be understood from the relative weights of the dz2- and dx2y2-contributions, as seen in Table II. Finally, the angular plot of the valence excitation in the B symmetry is a biconcave disk elongated along the z-axis, arising from the positive interference between s- and dz2-contributions contrary to the negative interference observed for its counterpart of the A symmetry.

Next, we compare these conventional first-order truncated ECD distributions with their full counterparts. For the valence, L3-edge, and L1-edge transitions, the anisotropic distributions virtually coincide, thereby confirming the validity of the first-order truncated description also for the anisotropic signal in these energy regimes. On the other hand, the full and truncated ECD distribution for the K-edge transitions is seen to have the same overall shape, but markedly different size, such that the factor-of-two overestimation at first order of the isotropic response (Tables I and S1) arises largely from an overall scaling. Closer inspection of the angular distribution of the full ECD distribution reveals that the lobes are not strictly perpendicular as in a dxy-orbital, hence indicating the contributions from solid harmonics of higher even angular momentum. To investigate this further, we provide the order-by-order contributions together with the full anisotropic ECD distribution in Fig. 3. Although Δf[3] ( = 0, 2, 4) resembles Δf[1], the inclusion of higher-order solid harmonics in the former leads to non-orthogonal lobes. Furthermore, the two distributions differ by an overall sign, such that the inclusion of Δf[3] decreases the ECD signal. A possible issue is the rate of convergence of the truncated interaction toward the full one. Our implementation of truncated interaction for linear absorption allows us to go to arbitrary order, a unique functionality of the Dirac code. In our previous work on the absorption of linearly polarized light,28 we demonstrated that the truncated treatment converges to the full interaction upon inclusion of higher-order terms, but for higher-energy transitions (photon energies beyond ∼1000 eV), the convergence behavior was too slow for practical applications. In the present case, we assessed the convergence of the isotropic differential oscillator strength expansion at the sulfur K-edge and found that the relative error is below the threshold of 1% at seventh order (Table III). Indeed, as shown in Fig. 3, the Δf[5] distribution ( = 0, 2, 4, 6) is minute. A similar convergence rate is found for the isotropic linear oscillator strength. From these data, it follows that the rate of convergence is acceptable for the application to the sulfur K-edge.

FIG. 3.

Convergence of the truncated anisotropic ECD distributions, i.e., Δf[2n+1](θ, ϕ), for the K-edge transitions of (a) A and (b) B symmetry in H2S2. The superscript notation Δf(→m) indicates accumulated contributions up through mth order. The black arrows point along the direction of the wave vector for the anisotropic ECD intensity given in Table S1.

FIG. 3.

Convergence of the truncated anisotropic ECD distributions, i.e., Δf[2n+1](θ, ϕ), for the K-edge transitions of (a) A and (b) B symmetry in H2S2. The superscript notation Δf(→m) indicates accumulated contributions up through mth order. The black arrows point along the direction of the wave vector for the anisotropic ECD intensity given in Table S1.

Close modal
TABLE III.

Contributions to the isotropic linear and differential oscillator strength (Δfiso[2n+1] and fiso[2n]) at various orders for H2S2, n = 0, 1, 2, 3, compared to the result of the full interaction for the two 1s1/2σSS* transitions of the A/B symmetry. The superscript notation □(→m) indicates accumulated contributions up through mth order. In particular, the “ac”-labeled row indicates the (differential) oscillator strengths accumulated through 8th (9th) order. The errors upon truncation are defined as %δΔfiso(2n+1)=|(Δfiso(2n+1)Δfisofull)/Δfisofull|×100%. The results were obtained using the 4c-TD-PBE0 level of theory and the uncontracted aug-pcx-3/aug-pc3 basis set. The numbers in parentheses are exponents of 10.

4A(1s1/2σSS*);ωA=2427.88349 (eV)
nfiso[2n]Δfiso[2n+1]%δfiso(2n)%δΔfiso(2n+1)
5.5576(−05) 1.6569(−04) 97.99 109.10 
3.4276(−03) −1.0818(−04) 26.07 27.42 
−8.1305(−04) 2.4151(−05) 3.36 3.06 
9.9900(−05) −2.2416(−06) 2.57(−1) 2.32(−1) 
−7.5661(−06) 1.2938(−07) 1.69(−2) 3.96(−1) 
ac 2.7625(−03) 7.9240(−05) ⋯ ⋯ 
full 2.7630(−03) 7.9240(−05) ⋯ ⋯ 
n 4B(1s1/2σSS*);ωB=2427.83334 (eV) 
1.0512(−02) −1.6570(−04) 35.24 109.10 
−3.4601(−03) 1.0818(−04) 9.27 27.42 
8.1318(−04) −2.4153(−05) 1.19 3.06 
−9.9897(−05) 2.2413(−06) 9.60(−2) 2.24(−1) 
7.5513(−06) −1.2878(−07) 1.43(−3) 3.98(−1) 
ac 7.7729(−03) −7.9433(−05) ⋯ ⋯ 
full 7.7728(−03) −7.9246(−05) ⋯ ⋯ 
4A(1s1/2σSS*);ωA=2427.88349 (eV)
nfiso[2n]Δfiso[2n+1]%δfiso(2n)%δΔfiso(2n+1)
5.5576(−05) 1.6569(−04) 97.99 109.10 
3.4276(−03) −1.0818(−04) 26.07 27.42 
−8.1305(−04) 2.4151(−05) 3.36 3.06 
9.9900(−05) −2.2416(−06) 2.57(−1) 2.32(−1) 
−7.5661(−06) 1.2938(−07) 1.69(−2) 3.96(−1) 
ac 2.7625(−03) 7.9240(−05) ⋯ ⋯ 
full 2.7630(−03) 7.9240(−05) ⋯ ⋯ 
n 4B(1s1/2σSS*);ωB=2427.83334 (eV) 
1.0512(−02) −1.6570(−04) 35.24 109.10 
−3.4601(−03) 1.0818(−04) 9.27 27.42 
8.1318(−04) −2.4153(−05) 1.19 3.06 
−9.9897(−05) 2.2413(−06) 9.60(−2) 2.24(−1) 
7.5513(−06) −1.2878(−07) 1.43(−3) 3.98(−1) 
ac 7.7729(−03) −7.9433(−05) ⋯ ⋯ 
full 7.7728(−03) −7.9246(−05) ⋯ ⋯ 

To assess whether the truncated formalism is applicable at higher transition energies, we have performed an additional series of calculations of the K- and L1-edge isotropic differential oscillator strength for the heavier H2X2 analogs (X = Se and Te). As shown in Sec. S2, we find that only the selenide L1-edge converges at a sufficient rate, whereas it is too slow to practically converge for all other edges. In both cases, the first-order treatment not only overestimates the rotatory strengths of individual transitions but also most critically incorrectly predicts the signs of (the intense) pairs of near-degenerate transitions. It thereby appears that the breakdown of the first-order description at the sulfur K-edge, and not at the sulfur L-edges, is a consequence of its order-of-magnitude higher transition energy (∼215 and 2427 eV at the L1- and K-edges).

We have reported the first implementation and application of the anisotropic and isotropic ECD signal using the full semi-classical light–matter interaction operator within a four-component relativistic framework. This simultaneous account of beyond-first-order light–matter interactions and relativistic effects provides two additional sources of ECD, which become increasingly important at high photon energies. The linear form of the light–matter interaction operator in the relativistic domain further enabled straightforward extension to a multipole-based scheme in the velocity representation that allows for the traditional (albeit, in general, ambiguous) decomposition into electric and magnetic contributions while retaining order-by-order gauge-origin independence.

The presented approach was used to investigate the ECD response of two prototypical disulfides, H2S2 and (CH3S)2, across the electromagnetic spectrum, from valence to core transitions. To quantify the implications of higher-order effects, we compared the results of the full interaction to those obtained within the traditional lowest-order non-vanishing (i.e., first-order) truncated generalized velocity representation. Going beyond the electric-dipole approximation at the sulfur K-edge leads to non-negligible intensity redistribution among near-degenerate transitions but with no visible implications on the linear absorption profile. On the other hand, the differential absorption profile is affected by such redistribution because of its signed nature. This leads to an overall factor-of-two overestimation. By examining the shapes of the underlying anisotropic ECD distributions, we find this discrepancy to largely originate from an overall scaling that is corrected upon introducing third-order contributions. Critically, the first-order treatment deteriorates at higher transition energies (beyond ∼1000 eV) where it may even fail to predict the sign of individual differential oscillator strengths. At such energies, going to higher orders is not a practical remedy because of the slow convergence of the truncated interaction—the full interaction is a must.

From a practical point of view, the full interaction operator provides a compact and inherently gauge-origin invariant treatment of light–matter interactions, which is both implementation- and computation-wise competitive with traditional multipole-based schemes for oriented and isotropic linear spectroscopies. For this reason, we believe that the use of the full light–matter interaction will become the new standard for theoretical x-ray spectroscopies. The question of how to efficiently handle the rotational averaging for nonlinear light–matter interactions will be explored in future work.

See supplementary material for linear and differential absorption spectra and corresponding tabulated oscillator strengths for H2S2, The analysis of the L1- and K-edges of the H2X2, X = S, Se, and Te series, and the analysis of the structure and origin dependence of the rotational strength tensor.

N.H.L. acknowledges financial support from the Carlsberg Foundation (Grant No. CF16-0290). M.v.H. acknowledges funding from the French Ministry of Higher Education and Research. Computing time from CALMIP (Calcul en Midi-Pyrenées) and SNIC (Swedish National Infrastructure for Computing) at the National Supercomputer Center (NSC) are gratefully acknowledged.

There are no conflicts to declare.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
A.
Dorta-Urra
and
P.
Bargueño
,
Symmetry
11
,
661
(
2019
).
2.
L. D.
Barron
,
Molecular Light Scattering and Optical Activity
(
Cambridge University Press
,
2004
).
3.
P. E.
Schipper
and
A.
Rodger
,
J. Am. Chem. Soc.
105
,
4541
(
1983
).
4.
N. H.
List
,
J.
Knoops
,
J.
Rubio-Magnieto
,
J.
Idé
,
D.
Beljonne
,
P.
Norman
,
M.
Surin
, and
M.
Linares
,
J. Am. Chem. Soc.
139
,
14947
(
2017
).
5.
N.
Berova
,
K.
Nakanishi
, and
R. W.
Woody
,
Circular Dichroism: Principles and Applications
(
John Wiley & Sons
,
2000
).
6.
D.
Lindorfer
and
T.
Renger
,
J. Phys. Chem. B
122
,
2747
(
2018
).
7.
P.
Akhtar
,
D.
Lindorfer
,
M.
Lingvay
,
K.
Pawlak
,
O.
Zsiros
,
G.
Siligardi
,
T.
Jávorfi
,
M.
Dorogi
,
B.
Ughy
,
G.
Garab
,
T.
Renger
, and
P. H.
Lambrev
,
J. Phys. Chem. B
123
,
1090
(
2019
).
8.
M.
Arago
,
Mem. Inst. Fr.
1
,
93
(
1811
); available at https://babel.hathitrust.org/cgi/pt?id=ucm.5326746608&view=image&seq=105.
9.
J.-B.
Biot
,
Mem. Inst. Fr.
1
,
1
(
1812
); available at https://babel.hathitrust.org/cgi/pt?id=ucm.5326746626&view=page&seq=11.
10.
L.
Alagna
,
T.
Prosperi
,
S.
Turchini
,
J.
Goulon
,
A.
Rogalev
,
C.
Goulon-Ginet
,
C. R.
Natoli
,
R. D.
Peacock
, and
B.
Stewart
,
Phys. Rev. Lett.
80
,
4799
(
1998
).
11.
J.
Goulon
,
C.
Goulon-Ginet
,
A.
Rogalev
,
V.
Gotte
,
C.
Malgrange
,
C.
Brouder
, and
C. R.
Natoli
,
J. Chem. Phys.
108
,
6394
(
1998
).
12.
J.
Goulon
,
C.
Goulon-Ginet
,
A.
Rogalev
,
V.
Gotte
,
C.
Malgrange
, and
C.
Brouder
,
J. Synchrotron Radiat.
6
,
673
(
1999
).
13.
B.
Stewart
,
R. D.
Peacock
,
L.
Alagna
,
T.
Prosperi
,
S.
Turchini
,
J.
Goulon
,
A.
Rogalev
, and
C.
Goulon-Ginet
,
J. Am. Chem. Soc.
121
,
10233
(
1999
).
14.
M.
Tanaka
,
K.
Nakagawa
,
A.
Agui
,
K.
Fujii
, and
A.
Yokoya
,
Phys. Scr.
2005
,
873
.
15.
K.
Nakagawa
,
T.
Matsui
,
Y.
Izumi
,
A.
Agui
,
M.
Tanaka
, and
T.
Muro
,
Radiat. Phys. Chem.
78
,
1198
(
2009
).
16.
Y.
Izumi
,
A.
Imazu
,
A.
Mimoto
,
M.
Tanaka
,
K.
Nakagawa
,
M.
Tanaka
,
A.
Agui
, and
T.
Muro
,
J. Phys.: Conf. Ser.
190
,
012209
(
2009
).
17.
Y.
Izumi
,
M.
Tanabe
,
A.
Imazu
,
A.
Mimoto
,
M.
Tanaka
,
A.
Agui
,
T.
Muro
, and
K.
Nakagawa
,
J. Chem. Phys.
138
,
074305
(
2013
).
18.
J.
Goulon
,
F.
Sette
,
C.
Moise
,
A.
Fontaine
,
D.
Perey
,
P.
Rudolf
, and
F.
Baudelet
,
Jpn. J. Appl. Phys., Part 1
32
,
284
(
1993
).
19.
A.
Rogalev
,
J.
Goulon
,
F.
Wilhelm
, and
A.
Bosak
, in
Magnetism and Synchrotron Radiation
, edited by
E.
Beaurepaire
,
H.
Bulou
,
F.
Scheurer
, and
K.
Jean-Paul
(
Springer Berlin Heidelberg
,
Berlin, Heidelberg
,
2010
), pp.
169
190
.
20.
C.
Brouder
,
J. Phys.: Condens. Matter
2
,
701
(
1990
).
21.
J.
Goulon
, in
Rayonnement Synchrotron Polarisé, Electrons Polarisés et Magnétisme
, edited by
E.
Beaurepaire
,
B.
Carrière
, and
J.-P.
Kappler
(
IPCMS
,
Mittelwihr, Haut-Rhin
,
1990
), pp.
333
386
.
22.
S.
Turchini
,
N.
Zema
,
S.
Zennaro
,
L.
Alagna
,
B.
Stewart
,
R. D.
Peacock
, and
T.
Prosperi
,
J. Am. Chem. Soc.
126
,
4532
(
2004
).
23.
A. E.
Hansen
and
T. D.
Bouman
, “
Natural chiroptical spectroscopy: Theory and computations
,” in
Advances in Chemical Physics
(
John Wiley & Sons
,
2007
), pp.
545
644
.
24.
F.
London
,
J. Phys. Radium
8
,
397
(
1937
).
25.
K. L.
Bak
,
A. E.
Hansen
,
K.
Ruud
,
T.
Helgaker
,
J.
Olsen
, and
P.
Jørgensen
,
Theor. Chim. Acta
90
,
441
(
1995
).
26.
N. H.
List
,
J.
Kauczor
,
T.
Saue
,
H. J. A.
Jensen
, and
P.
Norman
,
J. Chem. Phys.
142
,
244111
(
2015
).
27.
N. H.
List
,
T.
Saue
, and
P.
Norman
,
Mol. Phys.
115
,
63
(
2017
).
28.
N. H.
List
,
T. R. L.
Melin
,
M.
van Horn
, and
T.
Saue
,
J. Chem. Phys.
152
,
184110
(
2020
).
29.
S.
Villaume
and
P.
Norman
,
Chirality
21
,
E13
(
2009
).
30.
Y.
Zhang
,
J. R.
Rouxel
,
J.
Autschbach
,
N.
Govind
, and
S.
Mukamel
,
Chem. Sci.
8
,
5969
(
2017
).
31.
O.
Takahashi
,
M.
Kimoto
, and
L. G. M.
Pettersson
,
Chem. Phys.
450-451
,
109
(
2015
).
32.
V.
Kimberg
and
N.
Kosugi
,
J. Chem. Phys.
126
,
245101
(
2007
).
33.
A.
Jiemchooroj
,
U.
Ekström
, and
P.
Norman
,
J. Chem. Phys.
127
,
165104
(
2007
).
34.
V.
Carravetta
,
O.
Plachkevytch
,
O.
Vahtras
, and
H.
Ågren
,
Chem. Phys. Lett.
275
,
70
(
1997
).
35.
O.
Plashkevych
,
V.
Carravetta
,
O.
Vahtras
, and
H.
Ågren
,
Chem. Phys.
232
,
49
(
1998
).
36.
J.
Kaniauskas
,
Liet. Fiz. Rink
14
,
463
(
1974
) (in Russian; English translation available from the present authors); available at http://www.itpa.lt/LFD/profile.html.
37.
I. P.
Grant
,
Relativistic Quantum Theory of Atoms and Molecules: Theory and Computation
(
Springer Science & Business Media
,
2007
), Vol. 40, pp.
673
713
.
38.
S.
Urban
,
E.
Herbst
,
P.
Mittler
,
G.
Winnewisser
,
K. M. T.
Yamada
, and
M.
Winnewisser
,
J. Mol. Spectrosc.
137
,
327
(
1989
).
39.
E.
Herbst
and
G.
Winnewisser
,
Chem. Phys. Lett.
155
,
572
(
1989
).
40.
G. S.
Maciel
,
P. R. P.
Barreto
,
F.
Palazzetti
,
A.
Lombardi
, and
V.
Aquilanti
,
J. Chem. Phys.
129
,
164302
(
2008
).
41.
V.
Aquilanti
,
M.
Ragni
,
A. C. P.
Bitencourt
,
G. S.
Maciel
, and
F. V.
Prudente
,
J. Phys. Chem. A
113
,
3804
(
2009
).
42.
J. M.
Thornton
,
J. Mol. Biol.
151
,
261
(
1981
).
43.
A.
Rauk
,
J. Am. Chem. Soc.
106
,
6517
(
1984
).
44.
T.-K.
Ha
and
W.
Cencek
,
Chem. Phys. Lett.
182
,
519
(
1991
).
45.
M.
Pericou-Cayere
,
M.
Rerat
, and
A.
Dargelos
,
Chem. Phys.
226
,
297
(
1998
).
46.
C.
Diedrich
and
S.
Grimme
,
J. Phys. Chem. A
107
,
2524
(
2003
).
47.
M.
Pecul
and
K.
Ruud
, in
Response Theory and Molecular Properties (A Tribute to Jan Linderberg and Poul Jørgensen)
, Advances in Quantum Chemistry, edited by
h. J. Aa.
Jensen
(
Academic Press
,
2005
), Vol. 50, pp.
185
212
.
48.
L.
Bednárová
,
P.
Bour
, and
P.
Malon
,
Chirality
22
,
514
(
2010
).
49.
M.
Scott
,
D. R.
Rehn
,
S.
Coriani
,
P.
Norman
, and
A.
Dreuw
,
J. Chem. Phys.
154
,
064107
(
2021
).
50.
DIRAC, a relativistic ab initio electronic structure program, Release DIRAC21 (
2021
), written by
R.
Bast
,
A. S. P.
Gomes
,
T.
Saue
,
L.
Visscher
, and
H. J. Aa.
Jensen
, with contributions from
I. A.
Aucar
,
V.
Bakken
,
K. G.
Dyall
,
S.
Dubillard
,
U.
Ekström
,
E.
Eliav
,
T.
Enevoldsen
,
E.
Faßhauer
,
T.
Fleig
,
O.
Fossgaard
,
L.
Halbert
,
E. D.
Hedegård
,
T.
Helgaker
,
B.
Helmich–Paris
,
J.
Henriksson
,
M.
Iliaš
,
C. R.
Jacob
,
S.
Knecht
,
S.
Komorovský
,
O.
Kullie
,
J. K.
Lærdahl
,
C. V.
Larsen
,
Y. S.
Lee
,
N. H.
List
,
H. S.
Nataraj
,
M. K.
Nayak
,
P.
Norman
,
G.
Olejniczak
,
J.
Olsen
,
J. M. H.
Olsen
,
A.
Papadopoulos
,
Y. C.
Park
,
J. K.
Pedersen
,
M.
Pernpointner
,
J. V.
Pototschnig
,
R.
Di Remigio
,
M.
Repisky
,
K.
Ruud
,
P.
Sałek
,
B.
Schimmelpfennig
,
B.
Senjean
,
A.
Shee
,
J.
Sikkema
,
A.
Sunaga
,
A. J.
Thorvaldsen
,
J.
Thyssen
,
J.
van Stralen
,
M. L.
Vidal
,
S.
Villaume
,
O.
Visser
,
T.
Winther
, and
S.
Yamamoto
(available at , see also http://www.diracprogram.org).
51.
T.
Saue
,
R.
Bast
,
A. S. P.
Gomes
,
H. J. Aa.
Jensen
,
L.
Visscher
,
I. A.
Aucar
,
R.
Di Remigio
,
K. G.
Dyall
,
E.
Eliav
,
E.
Fasshauer
,
T.
Fleig
,
L.
Halbert
,
E. D.
Hedegård
,
B.
Helmich-Paris
,
M.
Iliaš
,
C. R.
Jacob
,
S.
Knecht
,
J. K.
Laerdahl
,
M. L.
Vidal
,
M. K.
Nayak
,
M.
Olejniczak
,
J. M. H.
Olsen
,
M.
Pernpointner
,
B.
Senjean
,
A.
Shee
,
A.
Sunaga
, and
J. N. P.
van Stralen
,
J. Chem. Phys.
152
,
204104
(
2020
).
52.
S. E.
Braslavsky
,
Pure Appl. Chem.
79
,
293
(
2007
).
53.
A. E.
Hansen
and
J.
Avery
,
Chem. Phys. Lett.
13
,
396
(
1972
).
54.
S.
Bernadotte
,
A. J.
Atkins
, and
C. R.
Jacob
,
J. Chem. Phys.
137
,
204106
(
2012
).
55.
F.
Bloch
, in
W. Heisenberg und die Physik unserer Zeit
, edited by
F.
Bopp
(
Vieweg & Sohn
,
Braunschweig
,
1961
).
56.
D. H.
Kobe
,
Am. J. Phys.
50
,
128
(
1982
).
57.
A. M.
Stewart
,
J. Phys. A: Math. Gen.
32
,
6091
(
1999
).
58.
T.
Saue
, in
Relativistic Electronic Structure Theory. Part 1. Fundamentals
, edited by
P.
Schwerdtfeger
(
Elsevier
,
Amsterdam
,
2002
), p.
332
.
59.
P.
Norman
,
K.
Ruud
, and
T.
Saue
,
Principles and Practices of Molecular Properties: Theory, Modeling, and Simulations
(
Wiley
,
Hoboken, NJ
,
2018
).
60.
T. B.
Pedersen
and
A. E.
Hansen
,
Chem. Phys. Lett.
246
,
1
(
1995
).
61.
H.-G.
Kuball
,
G.
Sieber
,
S.
Neubrech
,
B.
Schultheis
, and
A.
Schönhofer
,
Chem. Phys.
169
,
335
(
1993
).
62.
T.
Saue
and
H. J. Aa.
Jensen
,
J. Chem. Phys.
111
,
6211
(
1999
).
63.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
,
J. Phys. Chem.
98
,
11623
(
1994
).
64.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
65.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
,
Can. J. Phys.
58
,
1200
(
1980
).
66.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
67.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
68.
D. E.
Woon
and
T. H.
Dunning
,
J. Chem. Phys.
98
,
1358
(
1993
).
69.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, gaussian 16, Revision C.01,
Gaussian, Inc.
,
Wallingford, CT
,
2016
.
70.
D. B.
Craig
and
A. A.
Dombkowski
,
BMC Bioinf.
14
,
346
(
2013
).
71.
G.
Pescitelli
,
L.
Di Bari
, and
N.
Berova
,
Chem. Soc. Rev.
40
,
4603
(
2011
).
72.
M.
Stener
,
G.
Fronzoni
, and
M.
de Simone
,
Chem. Phys. Lett.
373
,
115
(
2003
).
73.
C.
South
,
A.
Shee
,
D.
Mukherjee
,
A. K.
Wilson
, and
T.
Saue
,
Phys. Chem. Chem. Phys.
18
,
21010
(
2016
).
74.
M.
Kadek
,
L.
Konecny
,
B.
Gao
,
M.
Repisky
, and
K.
Ruud
,
Phys. Chem. Chem. Phys.
17
,
22566
(
2015
).
75.
T.
Fransson
,
D.
Burdakova
, and
P.
Norman
,
Phys. Chem. Chem. Phys.
18
,
13591
(
2016
).
76.
L.
Visscher
and
K. G.
Dyall
,
At. Data Nucl. Data Tables
67
,
207
(
1997
).
77.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
78.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
79.
M. A.
Ambroise
and
F.
Jensen
,
J. Chem. Theory Comput.
15
,
325
(
2019
).
80.
F.
Jensen
,
J. Chem. Phys.
115
,
9113
(
2001
).
81.
F.
Jensen
,
J. Chem. Phys.
116
,
7372
(
2002
).
82.
F.
Jensen
,
J. Chem. Phys.
117
,
9234
(
2002
).
83.
L.
Visscher
,
Theor. Chem. Acc.
98
,
68
(
1997
).
84.
I.
Tokue
,
A.
Hiraya
, and
K.
Shobatake
,
Chem. Phys.
130
,
401
(
1989
).
85.
C. H.
Hearn
,
E.
Turcu
, and
J. A.
Joens
,
Atmos. Environ., Part A
24
,
1939
(
1990
).
86.
K.
Schnorr
,
A.
Bhattacherjee
,
K. J.
Oosterbaan
,
M. G.
Delcey
,
Z.
Yang
,
T.
Xue
,
A. R.
Attar
,
A. S.
Chatterley
,
M.
Head-Gordon
,
S. R.
Leone
, and
O.
Gessner
,
J. Chem. Phys. Lett.
10
,
1382
(
2019
).
87.
M.
Ochmann
,
A.
Hussain
,
I.
von Ahnen
,
A. A.
Cordones
,
K.
Hong
,
J. H.
Lee
,
R.
Ma
,
K.
Adamczyk
,
T. K.
Kim
,
R. W.
Schoenlein
,
O.
Vendrell
, and
N.
Huse
,
J. Am. Chem. Soc.
140
,
6554
(
2018
).
88.
A. P.
Hitchcock
,
S.
Bodeur
, and
M.
Tronc
,
Physica B
158
,
257
(
1989
).
89.
J.
Linderberg
and
J.
Michl
,
J. Am. Chem. Soc.
92
,
2619
(
1970
).
90.
S.
Svensson
,
A.
Naves de Brito
,
M. P.
Keane
,
N.
Correia
, and
L.
Karlsson
,
Phys. Rev. A
43
,
6441
(
1991
).
91.
S.
Svensson
,
A.
Ausmees
,
S. J.
Osborne
,
G.
Bray
,
F.
Gel’mukhanov
,
H.
Ågren
,
A.
Naves de Brito
,
O.-P.
Sairanen
,
A.
Kivimäki
,
E.
Nõmmiste
,
H.
Aksela
, and
S.
Aksela
,
Phys. Rev. Lett.
72
,
3021
(
1994
).
92.
B. T.
Thole
and
G.
van der Laan
,
Phys. Rev. A
38
,
1943
(
1988
).
93.
P. S.
Bagus
,
H.
Freund
,
H.
Kuhlenbeck
, and
E. S.
Ilton
,
Chem. Phys. Lett.
455
,
331
(
2008
).
94.
R. B.
Bernini
,
L. B. G.
da Silva
,
F. N.
Rodrigues
,
L. H.
Coutinho
,
A. B.
Rocha
, and
G. G. B.
de Souza
,
J. Chem. Phys.
136
,
144307
(
2012
).
95.
Y.
Baba
,
K.
Yoshii
, and
T. A.
Sasaki
,
J. Chem. Phys.
105
,
8858
(
1996
).
96.
I. J.
Pickering
,
M.
Barney
,
J. J. H.
Cotelesage
,
L.
Vogt
,
M. J.
Pushie
,
A.
Nissan
,
R. C.
Prince
, and
G. N.
George
,
J. Phys. Chem. A
120
,
7279
(
2016
).
97.
G.
Bergson
,
Ark. Kemi
12
,
233
(
1958
).
98.
G.
Bergson
,
Ark. Kemi
18
,
409
(
1961
).
99.
G.
Wagnière
and
W.
Hug
,
Tetrahedron Lett.
11
,
4765
(
1970
).
100.
R. W.
Woody
,
Tetrahedron
29
,
1273
(
1973
).
101.
M.
Carmack
and
L. A.
Neubert
,
J. Am. Chem. Soc.
89
,
7134
(
1967
).
102.
R. M.
Dodson
and
V. C.
Nelson
,
J. Org. Chem.
33
,
3966
(
1968
).
103.
G.
Claeson
,
Acta Chem. Scand.
22
,
2429
(
1968
).
104.
U.
Ludescher
and
R.
Schwyzer
,
Helv. Chim. Acta
54
,
1637
(
1971
).
105.
J. L.
Campbell
and
T.
Papp
,
At. Data Nucl. Data Tables
77
,
1
(
2001
).
106.
C.
Nicolas
and
C.
Miron
,
J. Electron Spectrosc. Relat. Phenom.
185
,
267
(
2012
).

Supplementary Material