Probing chirality across the electromagnetic spectrum with the full semi-classical light–matter interaction

Chirality is fundamental to life on Earth, and the origin of homochirality of natural amino acids and sugars remains an unresolved question.¹ 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.² 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.⁵ 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 Arago⁸ and Biot,⁹ its existence in the x-ray range was only established in 1998: The detection of x-ray natural circular dichroism (XNCD) in non-centrosymmetric crystals^10–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(L₁)-edges^18,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 samples^19,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 representation²³ 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-relativistic^26,27 and relativistic regimes.²⁸ 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 L₁- 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 rule^36,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 H₂S₂ and dimethyl disulfide (CH₃S)₂. Because of the low disulfide torsional barriers (∼6–11 kcal/mol^38–41), the two enantiomeric forms (P- and M-helix) cannot be resolved experimentally. However, the disulfide bridge is an important structural element in proteins,⁴² where it preferentially occurs in non-planar, chiral conformations (C₂ 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, H₂S₂ 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.¹⁸ 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

where appears the polarization vectors ϵ₁ and ϵ₂ that, together with the wave vector k = ke_k$(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.⁵² Such an electromagnetic wave is conventionally represented in the Coulomb gauge by the scalar and vector potentials,

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

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

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

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 $πωϵ0ℏc→2meωℏ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,

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

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

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.

In our previous work on the absorption of linearly polarized light, we considered two schemes involving truncated light–matter interaction.²⁸ 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.⁵⁴ Applying this to Eq. (6), we obtain the following expressions for the differential oscillator strength to even and odd order in the wave vector:

where

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.²⁸ 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

where appears electric and magnetic multipole operators,

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

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

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

in terms of a multipole operator,²⁸ 

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

It is possible, though, to combine the best of both worlds. In our previous work,²⁸ 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

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

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

being perfectly symmetric in all indices $j1…jnp$⁠, 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

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

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

where the components of the real vector $Xj1…jn[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.

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

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

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

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

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 C₂ 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.⁶¹ 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

in terms of spherical harmonics $Yℓm$ 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

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

In the lowest order, the isotropic differential oscillator strength becomes

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 H₂S₂ and (CH₃S)₂ were obtained using the B3LYP^63–66 exchange–correlation functional and the cc-pVTZ^67,68 basis set. Geometry optimizations were performed in Gaussian 16.⁶⁹ To mimic the χ₃ disulfide angle typical for protein structures,⁷⁰ we performed constrained geometry optimizations for χ₃ = −87°, corresponding to M-helical chirality.⁷¹ The restricted excitation window approach^72,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,⁷⁶ and an 86-point Lebedev grid (L_max = 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 PBE0^77,78 exchange–correlation functional and the uncontracted aug-pcX-3⁷⁹ and aug-pc-3^80–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.⁸³ 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 a₀) along the C₂ 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 (CH₃S)₂ 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 (CH₃S)₂ for which experimental gas-phase absorption spectra are available.^84–88 Expectedly, and as shown in Fig. S1, the spectral profiles for H₂S₂ 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 (CH₃S)₂, 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.

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 H₂S₂.⁸⁹ Since the valence orbitals are found to have well-defined spin, we adopt a non-relativistic state notation for the valence states, i.e., 1¹B and 2¹A, respectively. They are separated by ∼0.15 eV. The second valence band originates from transitions into the $σCS*$ orbital (3¹A and 2¹B).

Turning to the x-ray region, the first two bands at the sulfur L_2,3-edge [Fig. 1(b)] are dominated by transitions from the symmetric and antisymmetric combinations of the S 2p_3/2 and 2p_1/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 L₃- and L₂-branches is comparable to the experimental splitting reported for dihydrogen sulfide (∼1.2 eV^90,91). The L₃/L₂ 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 j − j 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 L₂ 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,⁸⁶ originates from excitations to higher-lying orbitals of mixed $σCS*$/Rydberg character. Not surprisingly, the L₁- 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., C₂) 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 Michl⁸⁹ 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 C₂-rule for general chromophores of effective C₂ symmetry.⁹⁹ Woody extended the theoretical analysis to also include the absolute sign of the lowest ECD band within each quadrant,¹⁰⁰ 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.¹⁰⁰ The intensity asymmetry of the lowest-energy valence couplet is attributed to different intrafragment (i.e., CH₃S–) 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).

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 TiCl₄.²⁸ 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 H₂S₂. 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 C₂-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.

From symmetry considerations detailed in the supplementary material, we find that contributions from d_xz and d_yz 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 L₃-edge excitations of A symmetry, this shape is modulated by the $dx2−y2$-contribution, giving the shape of a biconcave disk elongated along the y- and x-axis, respectively, depending on its relative sign. For the L₁- and K-edge excitations, the d_xy-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 d_xy 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 L₃ 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 $dx2−y2$-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, L₃-edge, and L₁-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 d_xy-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,²⁸ 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.

To assess whether the truncated formalism is applicable at higher transition energies, we have performed an additional series of calculations of the K- and L₁-edge isotropic differential oscillator strength for the heavier H₂X₂ analogs (X = Se and Te). As shown in Sec. S2, we find that only the selenide L₁-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 L₁- 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, H₂S₂ and (CH₃S)₂, 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 H₂S₂, The analysis of the L₁- and K-edges of the H₂X₂, 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.