X-ray linear and non-linear spectroscopy of the ESCA molecule

Recent developments in multidimensional electronic spectroscopy, giving access to extremely short laser pulses with remarkable phase stability, have provided the tools required to follow the ultrafast dynamics of molecular systems at an unprecedented level of accuracy. The rapid development of X-ray free-electron laser (FEL) facilities^1–3 holds the promise to extend this approach to the rapidly growing field of novel multidimensional nonlinear spectroscopic techniques, ranging from two dimensional (2D) X-ray four wave mixing spectroscopy, X-ray Stimulated Raman,⁴ to time resolved X-ray diffraction.⁵ X-ray FELs offer unprecedented temporal and spatial resolution, making them the ideal choice for nonlinear spectroscopy of core excitations, electronic coherences, quantum control, and nonlinear dynamics. Besides the high temporal resolution, which can be pushed down to the attosecond time scale, X-ray techniques are of interest because of their element- and site-selectivity: different elements are spectroscopically resolved due to the large difference in binding energy of core electrons (∼100 eV difference for C, N, and O). Moreover, it is possible to distinguish the signals of the same element within a molecule, by detecting the small differences in the binding energy caused by the composition of its immediate surrounding. This chemical shift makes X-rays an adequate tool for structural analysis,^6,7 like heteronuclear NMR.

Recently, the use of X-ray techniques for obtaining excited-state molecular dynamics has grown considerably.⁸ In these experiments, a molecular target is excited optically, creating a valence electronic state wave-packet. After some delay, the excited state electronic structure is probed by an X-ray pulse. Several methods are available to track the photoinduced dynamics in the X-ray domain: X-ray photoelectron spectroscopy (XPS), for example, employs X-rays to ionize a core electron from a molecule and measure the energy of the ejected electron. X-ray absorption near edge structure (XANES) spectroscopy is used to elucidate core-to-valence excitations.⁹ X-ray emission spectroscopy (XES) studies the characteristic wavelengths of photons emitted by valence electrons when filling a previously generated core.¹⁰ It is also possible to image the excited state dynamics using off-resonant X-ray sum-frequency generation diffraction.^11,12

Multidimensional nonlinear spectroscopic techniques in the X-ray regime are expected to provide a novel diagnostic tool for tracing electronic and structural dynamics in molecular materials. While these are not yet available in current experimental facilities, the creation of new modeling capabilities and the development of efficient simulation protocols for the description of single-core and multiple-core excited state energetics and dynamics and for the interpretation of their spectroscopies are a timely and growing research field.^13–23 

For the work outlined in this article, we compute and investigate the signal contributions to the X-ray linear and nonlinear spectra of a paradigmatic system, namely, the ESCA molecule, to show how different core electronic excitations are coupled in a molecule and explain their energetic positioning and brightness. The ESCA molecule (ethyl trifluoroacetate, CF₃–CO–O–CH₂–CH₃, shown in Fig. 1), has been widely studied in X-ray spectroscopy^7,24–29 due to its small size and unique X-ray photoelectron core spectra. Since the ESCA molecule has four carbon atoms each in a different electronic environment, the XPS peak for each carbon is separated by 2 eV, such that each carbon can be excited independently within the same energy range. As such, the ESCA molecule offers the ability to study how nearby core electrons couple to valence excitations, without the need to change the probe energy by hundreds or thousands of electronvolt. The carbon K-edge is accessible by the seeded X-ray Free-Electron Lasers (XFEL) at Fermi, where coherent control of multiple pulses at short wavelength has been demonstrated by operating the FEL in the so-called echo-enabled harmonic generation (EEHG).^30,31 This recent achievement allows for a complete manipulation of the energy (from 200 eV up to 400 eV),³⁰ polarization, phase, and intensity of the produced pulses and have the potential to set the stage for the experiments proposed in this article, i.e., coherent four wave mixing X-ray nonlinear optical techniques.

Molecular core transitions absorb in the hundreds of electronvolt spectral range and are thus preceded by a vast number of valence states. In fact, core transitions are immersed in the continuum of valence ionized transitions. A common practical strategy for targeting specifically desired core-excitations is to neglect the coupling to the valence transitions (an approximation known as core-valence separation³²) and to the continuum.³³ Within this approximation, we apply linear response time-dependent density functional theory (LR-TDDFT)^34,35 and restricted active space self consistent field (RASSCF) methods corrected by second order perturbation (RASPT2). Specifically, at the TDDFT level, core-excitations^36–44 are targeted via the restricted excitation window (REW) technique.^37,42,45,46 At the RASSCF level, the goal is achieved by including the core orbitals of the atoms of interest in the active space (AS), i.e., the list of orbitals and electrons of paramount importance for describing the electronic structure of the system, in combination with a projection technique based on a selective removal of all configuration state functions (CSFs) with a certain occupation from a given subspace.¹⁹ 

We present linear and third-order nonlinear 2D X-ray absorption, XANES, and two-dimensional coherent X-ray spectroscopy (2DCXS), respectively, at the carbon K-edge (∼300 eV). For the latter, we simulate 2D four wave mixing spectra using the sum-over-state (SOS) approach.⁴⁷ In theoretical simulations, one may selectively switch on and off given contributions. We had exploited this possibility in the analysis, disentangling signals pertaining different core (pairs) which may overlap in a real experimental spectrum.

Since TDDFT cannot adequately describe double excitations,⁴⁸ only single core-excitations are treated at this level. On the contrary, the projection technique in combination with RASSCF is applicable to single and multiple core-transitions. Here, the computational complexity arises due to the sheer number of double core-excitations that must be simultaneously accounted for. To give an idea, already ten virtual orbitals in the active space give rise to 200 double core-transitions per core pair. In the case of ESCA with four carbon centers (and six different core pairs), this would require to consider more than 1500 double core-transitions coupled to 40 single-core transitions. Instead of a brute force approach, we focus on the simulation and detailed analysis of the 2DCXS signal for selected pairs of carbon cores by including a limited number of representative virtual orbitals (i.e., of π−, σ− and Rydberg character) in the active space. In this way, we categorize the signals that would be observed in the realistic experiment according to their spectral position and intensity. Furthermore, we quantify the spatial extension of the couplings between individual core-transitions and the associated spectral shifts. The findings represent the first step toward establishing 2DCXS as the optical counterpart of 2D-NMR.

Weak excitations such as absorption spectroscopy can be described by perturbative quantum calculations. Here, the applied fields are weaker than the internal molecular fields, and the XANES can be easily obtained. We have employed two methods for the computation of single- and double-core excitations, namely, TDDFT and RASSCF/RASPT2.

Based upon its balance between speed and accuracy, linear-response (LR) TDDFT^34,35 is a popular weak-field method that involves computing roots in the frequency-dependent molecular response, provided only small changes to electron density from the ground state.^49,50 TDDFT does not reproduce double excitaitons, which limits its ability to capture nonlinear processes. With care, however, TDDFT has been widely applied to calculate challenging processes such as double excitations,^51,52 core ionization dynamics,⁵³ ultrafast charge migration,⁵⁴ strong-field ionization,^55–58 charge transfer,^59–62 and X-ray absorption.^36–44 

The Casida approach is often used to simulate core absorption spectra with TDDFT.⁴⁹ In this framework, a set of Davidson iterations are performed, yielding the eigenvalues and eigenvectors of the linear response matrix.⁶⁵ These are then used to calculate the frequencies and oscillator strengths of the absorption peaks.⁴⁹ Since the entire history of the electron density is not known, TDDFT calculations often employ the adiabatic approximation,^34,66,67 thereby assuming that the density does not change much over time. Although weak perturbations effectively describe single excitations, the system cannot be resonantly driven to the excited state with adiabatic TDDFT,^68–70 and doubly excited states are not adequately described and strongly depend on the initial conditions.^71,72 Some novel methods attempt to address these challenges directly such as nonadiabatic DFT^73,74 and added memory effects.^75–79 Others employ two-step calculations with well-defined initial conditions, which have shown notable success in describing the double excitations and excited-state dynamics.^53,80,81

In this work, we use LR-TDDFT to compute the XANES spectra for the four carbon K-edges in the ESCA molecule. In addition, the nonlinear 2D X-ray spectrum is computed using the sum-over-states (SOS) approach⁴⁷ for the carbon cores in the ESCA molecule. Based solely on the single excitations, since TDDFT does not provide the double excitaitons, we simulated ground-state bleaching (GSB) and stimulated emission (SE) contributions (see Sec. II C and the diagrams in the  Appendix, Fig. 6).

The following protocol was employed to compute the singly excited states used to obtain the linear K-edge core absorption and the above-mentioned selected contribution of the 2DCXS signal:

1. the ground state geometry is optimized for an isolated ESCA molecule using the PBE0 functional and a mixed basis of aug-cc-pVTZ for all atoms except for carbon; the carbon atoms used the Sapporo-TZP basis set to provide a better description of the core and associated excitations; the PBE0 functional was used in the TDDFT module of the development version of the NWChem computational chemistry software package;⁸² 

2. for this geometry, mixed basis, and functional, the states contributing to the linear K-edge absorption spectrum were computed using the linear response TDDFT module of NWChem with a restricted excitation window (REW)^37,42,45,46 for the carbon cores; the states associated with each carbon core were then isolated to facilitate the analysis of the linear core absorption peaks.

All calculations were performed in the gas-phase employing a development version of the NWChem computational chemistry software package.⁸² 

The restricted active space self consistent field (RASSCF) approach,⁸³ from the family of multiconfiguration wavefunction based methods, offers a neat way to target desired core excitations, which are indeed placed hundreds of electronvolt above a vast number of other states: single-core transitions lie above all the preceding valence excitations, while double-core excitations appear only above an extremely large number of both valence and single-core excitations.

RASSCF subdivides the full active space (AS) into three subspaces: RAS1, with a fixed upper limit of holes; RAS2, where all possible permutations of electrons within the orbitals are considered; RAS3, with a maximal number of electrons. The possible configurations over which the wavefunction is expanded are built according to these subspaces rules. RASSCF offers a higher flexibility with respect to the more commonly employed complete active space approach (CASSCF) which does not distinguish between subspaces,^84,85 thus effectively removing “dead wood” from the list of configuration state functions (CSFs) without sacrificing accuracy. The missing (dynamical) correlation can be added on top of a RASSCF calculation through second-order perturbation corrections (RASPT2).^86–88 

The strategy to treat single- and double-core excitations within the RASSCF/RASPT2 framework is based on in a novel projection technique recently implemented in OpenMolcas,⁸⁹ based on a selective removal of all CSFs with a certain occupation from a given subspace. The highly excited states (HEXS) method sets to zero the CI coefficients of all CSFs with the maximum occupation from the selected subspace(s), thus effectively projecting out all undesired valence transitions preceding energetically the core-transitions.¹⁹ The HEXS approaches offer intriguing possibilities in combination with the generalized version of the RASSCF technique, known as GASSCF,^90,91 which removes the restriction of working with only three subspaces. Applications of the HEXS technique can be found in Refs. 13–15.

Below we outline the protocol for computing both linear and nonlinear spectra (relying on single core-virtual excitations in the former case, and on both single and double core-virtual excitations in the latter). To speed up the calculations, we take advantage of the C_S symmetry of the molecule.

• A Hartree-Fock calculation is performed;

• The four cores active spaces are constructed as follows (see Fig. 2):

  • •

    RAS1: one core-orbital is placed herein and kept frozen to avoid orbital rotation (using the OpenMolcas SUPSYM keyword); the upper limit of holes is set to 1;

  • •

    RAS2: another core-orbital is placed herein and kept frozen as in RAS1;

  • •

    RAS3: 10 virtual orbitals of either a′ or a″ symmetry are included, together with one occupied valence orbital per symmetry; these occupied orbitals are required as initial tests omitting occupied valence orbitals demonstrated considerable instabilities in the subsequent second-order perturbative correction (RASPT2). The maximum of excitations in this subspace was set to four.

In the following, we designate the active space as RAS(6, 1, 4; 1, 1, 11), where the first three indices denote the number of electrons, the upper limit of holes (in RAS1), and the upper limit of excitations (in RAS3), respectively; while the last three indies denotes the number of orbitals included in the three subspaces: in order, RAS1, RAS2, and RAS3. Two separate computations were performed with this active space setup, one for each of the two symmetries: in the a′ computations, the RAS3 subspace is therefore filled with one occupied valence orbital and 10 virtual orbitals of a′ symmetry, while in the a″ set of computations, the RAS3 orbitals are of a″ symmetry. Note that the XANES protocol does not require per se the inclusion of two core orbitals in the active space. In this work, this is done to assure uniformity with the protocol for simulating nonlinear spectroscopy (see below).

• (c)

  RASSCF calculations were performed separately for each of the four carbons for the following states:

  • •

    Single core-to-virtual excitations: for each of the two symmetries a′ and a″, the HEXS keyword was utilized to force the creation of the core-hole, either from RAS1 or from RAS2, and up to 10 states were included in the state-averaging, i.e., SA-10-core-RAS(6, 1, 4; 1, 1, 11);

  • •

    GS: a state specific calculation was performed for each active space (i.e., for each chosen couple of core that are placed in RAS1 and RAS2).

The choice of number of states in the state-averaging procedure aims at describing at least one core-virtual transition per virtual orbital. However, double core + valence → virtual excitations may occasionally appear among the computed states. Due to being essentially dark from the GS, these states do not affect the outcome of the simulations.

• (d)

  Single state (SS)-RASPT2 was performed on top of each RASSCF calculation to account for the dynamical correlation, setting the ionization-potential electron-affinity shift^92,93 of 0.0 hartree; to reduce problems with intruder states, an imaginary shift⁹⁴ of 0.3 hartree has been applied;

• (e)

  Transition dipole moments between the GS and the single core-virtual manifolds were computed through the RAS state interaction (RASSI) routine.

• (a)

  A Hartree-Fock calculation is performed;

• (b)

  For each of the three selected core orbitals pairs C₁–C₂, C₁–C₃, and C₁–C₄ the active space (see Fig. 2) is constructed as follows:

  • •

    RAS1: the core-orbital of C₁ is placed herein and kept frozen to avoid orbital rotation (using the OpenMolcas SUPSYM keyword); the upper limit of holes is set to 1;

  • •

    RAS2: the core-orbital of the opposite carbon in the pair (C₂, C₃, or C₄) is placed herein and also kept frozen;

  • •

    RAS3: one occupied valence orbital per symmetry is placed herein; furthermore, a limited number of two virtual orbitals of each symmetry, i.e., a′ and a″, are included with a maximum of six excitations; the nature of the orbitals is discussed in Sec. III D.

In the following, we designate the active space as RAS(8, 1, 6; 1, 1, 6) where RAS3 contains three orbitals of a′ and a″ symmetry, respectively.

• (c)

  RASSCF calculations were performed separately for each of the three pairs C₁–C₂, C₁–C₃, and C₁–C₄ for the following sets of states:

  • •

    Double core-to-virtual excitations: between 10 and 15 states of a′ and a″, symmetry was computed in a state-averaging procedure, applying the HEXS keyword to both RAS1 and RAS2 to create a double core-hole;

  • •

    Single core-to-virtual excitations: between two and three states of either a′ and a″, symmetry was computed in a state-averaging procedure, applying the HEXS keyword either to RAS1 or RAS2 to create a single-core-hole;

  • •

    GS: a state specific calculation was performed;

Molecular orbitals were allowed to freely relax in the double core-to-virtual calculations. Subsequently, orbital relaxation was restricted only within the active space when computing single core-to-virtual excitations, thus effectively prohibiting rotations out of the active space. This strategy is necessary to assure active space uniformity and allows for the automation of the computational protocol when working with a limited number of virtual orbitals in the active space. We note that a larger number of virtual orbitals would ensure the completeness of the K-edge spectrum signal and would not require to impose restrictions on the orbital relaxation procedure.

The choice of number of states in the state-averaging procedure aims at describing all single and double core-virtual transitions that can be accessed within the given active space and give rise to various signal contributions in the nonlinear spectrum [i.e., ground-state bleaching, stimulated emission, and excited state absorption (ESA); see Sec. II C].

• (d)

  Multistate (MS)-RASPT2 was performed on top of each RASSCF calculation to account for the dynamical correlation, setting the ionization-potential electron-affinity shift^92,93 of 0.0 hartree; to reduce problems with intruder states, an imaginary shift⁹⁴ of 0.3 hartree has been applied;

• (e)

  Transition dipole moments were computed between the GS and the single-core manifolds, as well as between the single- and double-core manifolds through the RASSI routine.

In both linear- and nonlinear spectroscopy protocols, scalar relativistic effects have been included by using a second-order Douglas-Kroll-Hess Hamiltonian in combination with a relativistic atomic natural orbital basis set (ANO-RCC).⁹⁵ A triple-ζ basis function set augmented with two sets of d-functions was used for carbon, oxygen, and fluorine atoms, whereas an additional set of f-functions was used to improve the description of carbon, thus giving rise to ANO-RCC: C[4s3p2d1f], O[4s3p2d], F[4s3p2d], and H[2s1p]. A density-fitting approximation of the electron repulsion integrals has been used, knows as Cholesky decomposition.⁹⁶ 

All calculations were performed in the gas-phase with the OpenMolcas suite.⁸⁹ 

In the SOS approach, we divided the manifold of states into a (single) ground state (GS—g), a single excitation manifold (containing all the states e that are reached from the GS, i.e., the single-core excitations), and a double excitation manifold (containing all the states f that can be reached from the first excitation manifold, i.e., the double-core excitations). Transitions from the GS to the first excitation manifold have energy gap ω_eg and transition dipole moment μ_eg, while ω_fe and μ_fe represents transitions from first to second excitation manifold. XANES spectra are produced by employing only GS and first excitation manifold, while nonlinear (third order) techniques also need the second excitation manifold.

The quantities that enter in the simulation of spectra are energy differences (ω_eg and ω_fe), transition probabilities (transition dipole moments, μ_eg and μ_fe), and line-shape broadening. The former two quantities were provided by the two different level of theory (TDDFT and RASSCF/RASPT2, respectively), and the stick spectrum was dressed with a Lorenzian broadening. This broadening is the natural choice as the main contribution of the line-shape broadening in X-ray spectra should come from the extremely short lifetime of the core-excited states. The linear absorption A(ω) can be obtained from the collection of ω_eg and μ_eg values as

Here, γ is the line-width parameter governing the width of the broadening, which was set to 0.5 eV in all the simulations employed here. This corresponds to a full-width-at-half-maximum broadening of 1.0 eV.

For the simulation of 2DCXS signals, we consider coherent all X-ray four-wave mixing processes, carried out by subjecting the molecule to a sequence of three pulses, with a controlled time delay between each other. The pulses have wave-vectors k₁, k₂, and k₃ and carrier frequency ω₁, ω₂, and ω₃, respectively. A fourth pulse, k₄ − ω₄, allows for the heterodyne detection of the signal. The coherent nonlinear response generated in the k_I = −k₁ + k₂ + k₃ and in the k_II = k₁ −k₂ + k₃ phase matching directions, called rephasing and nonrephasing contributions, respectively, is recorded as a function of the delays t₁, t₂, and t₃ between the consecutive pulses. In this paper, we computed both rephasing and nonrephasing contributions to the 2DCXS maps at time t₂ = 0, and sum them together to obtain the (quasiabsorptive) spectra shown in Sec. III. 2DCXS is obtained as a two-dimensional frequency correlation plot along excitation and detection frequencies Ω₁ and Ω₃, by a Fourier transform of the signal with respect to the times t₁ and t₃. Expressions for the response functions and Feynman diagrams of the different contributions are given in the  Appendix.

The various contributions to the 2DCXS signal are denoted ground-state bleaching (GSB), stimulated emission (SE), and excited-state absorption (ESA). Their names reflect the different physical processes activated by the interaction with the sequence of pulses.⁹⁷ 

We now present the results obtained, at the two level of theory, by simulating both linear and nonlinear X-ray spectra.

The carbon K-edge linear absorption (XANES) spectra are shown in Fig. 3 for TDDFT (a) and RASSCF/RASPT2 (b). The spectra share the same overall spectral shape, characterized by a number of peaks originating from the different core-specific excited states. Notably, the spectra exhibit the same total signal intensity. The molecular orbital character of the core transitions is consistent. In general, the strength of a given core-to-virtual transition is determined by the extent of the spatial overlap between the starting (core) and the arrival (virtual) orbitals. Accordingly, the lowest bright transition in both calculations (lying at ∼278 eV and 288 eV in TDDFT and RASSCF/RASPT2 spectra, respectively) is of 1s(C2) → π^* origin where the arrival orbital of π^* character has the highest overlap with the core orbital of C2 (orbitals shown in Fig. 2). Further peaks, associated with core-to-σ* transitions, are distributed over a broad range of ∼10 eV depending on the starting core-orbital and the type of arrival σ^* orbital (i.e., C–H, C–C, C–O, or C–F). Rydberg-type orbitals are diffuse and exhibit only a small overlap with the core-orbitals. Consequently, core-to-Rydberg transitions are weak to dark, giving rise to a uniform nonzero background. Due to the various arrival orbitals available to the different cores in the here simulated XANES spectra, each core contributes to several bands(see core-specific spectra in Fig. 3). This is at variance with XPS, where each of the four cores contributes to a different band, well separated from the rest.⁷ A quantitative analysis of the underlying signal is provided in Table I.

Comparing the different cores, the two levels of theories are in excellent agreement for the C1 and C2 carbons. We observe however notable discrepancies in the C3 and C4 edges. These are likely due to well-known Coulomb and self-interaction errors,^98–101 which strongly depend on the choice of DFT functional. The choice of functional also affects the overall energy shift of the core (Δ ∼ 10 eV relative to the RASSCF/RASPT2 results). As shown in the figure, it has small effects on the spacing and relative intensities within the edge, but the overall composition and spectral information remains largely consistent. Using a tuned, range-separated hybrid functional^41,102 and other established methods^103,104 have been shown to improve such errors, provided additional computational resources. We note, however, that the TDDFT XANES results at our reported level of theory are sufficient for the comparison with the high-level RASSCF/RASPT2 calculations.

We should mention here that the higher the energy of a considered core-excited state, the higher the probability of other phenomena to occur, as, e.g., photo-ionization of the excited core electron. This will translate in a reduced intensity and increased broadening (due to reduced lifetime) of the signals lying close/slightly above the core-ionization threshold. These effects, as well as the continuum of photo-ionized states, are not included in the present calculations.

In Fig. 4, we report a sequence of all X-ray four-wave mixing quasiabsorptive 2D maps simulated at t₂ = 0 (see the  Appendix). To reduce the computational cost, instead of simulating the full four-center 2D map, we present 2D maps of three subsystems, each comprising the core-orbital of C1, in combination with C4 [Figs. 4(a)–4(c)], C3 [Figs. 4(d)–4(f)], and C2 [Figs. 4(g)–4(i)]. Furthermore, we consider a restricted number of orbitals and of excited states, as detailed in Sec. II B 2. These simplifications facilitate the assignment of the peak patterns and establish trends which we expect to recover in the total spectrum. In particular, the choice of the three subsystems is aimed at exploring how core-orbital spatial proximity affects the spectral signatures of the various core-excitations. In this framework, C1–C2 represents an adjacent core couple, C1–C4 represents a distant core couple, and C1–C3 represents an intermediate situation.

In general, one expects to see in a 2D spectrum a pattern of diagonal and off-diagonal peaks of different signs, where positive signals pertain GSB and SE contributions, associated with transitions between the GS and singly excited manifold, whereas negative signals refer to ESA contributions, associated with transitions between the singly and doubly excited manifolds (see the  Appendix). In Fig. 4, we show separately the GSB-SE contribution (left column) and the ESA contribution (central column), whereas the total signal is presented in the right column. The XANES spectrum for the given pair of atoms (i.e., a sub-spectrum of the total XANES in Fig. 3) is also reported along both the excitation and detection energy axes, enabling one to identify the core-specific origin of the various peaks and cross-peaks in the maps. A mapping of the core-excited states computed in this restricted setup to those obtained in the single-core excitations with the larger active space is reported in Table I.

Along the diagonal, GSB and SE peaks mirror transitions which are already detected in the XANES spectrum and correspond to pumping and probing the same core-virtual transition. We note that in third-order nonlinear spectroscopy, transition dipole moments enter to forth power (compared to second power in XANES), thus effectively suppressing weak (such as core-to-Rydberg) transitions. The GSB and SE spectra reveal a number of off-diagonal peaks as well, created by addressing different cores with the pump and the probe pulses. In fact, an off-diagonal bleach signal emerges for each couple of core-transitions. Some of these bleach signals are expected to be canceled by corresponding ESA signals when the coupling between the transitions is negligible or null. This is observed in the spectrum of the distant couple C1–C4 [Figs. 4(a)–4(c)]. On the diagonal, we see an intense peak at 295.6 eV belonging to a C1-σ^*(C1–F) excitation (labeled 2) and a weak peak at around 288 eV due to two C4-Ry^* excitations (labeled 1′). The two centers are distant enough (∼4.7 Å) so that C1 and C4 transitions are independent of each other. As a consequence off-diagonal signals coupling, the C1 and C4 peaks are absent in the total map due to the cancellation of the positive GSB-SE and the negative ESA signals (highlighted by the black circles Fig. 4).

The second row of plots [Figs. 4(d)–4(f)] shows the signals associated with the C1–C3 couple. On the diagonal, we see an intense peak at 289.6 eV belonging to a C3–σ^*(C3–O) excitation (labeled 1) and two weaker peaks at 293.5 eV (labeled 2) and 298.5 eV (labeled 3) due to transitions from C1 to π^*(C2) and σ^*(C3–O), respectively. As the two cores lie spatially closer (∼3.7 Å) with respect to the C1–C4 couple, the obtained spectra display clear, albeit weak, off-diagonal features indicative of the intercore coupling. No GSB-SE/ESA cross-peak cancellation is observed for transitions 1 and 3 and the ESA appear red-shifted by ∼3.5 eV with respect to the off-diagonal bleach signal (peaks at Ω₁/Ω₃ = 298.5/286.0 eV and peaks at Ω₁/Ω₃ = 289.6/295.0 eV). Instead, we observe a near complete cancellation of the cross-peaks between transitions 1 and 2 [highlighted by black circles in Figs. 4(d)–4(f)].

The third row [Figs. 4(g)–4(i)] shows the signals associated with the adjacent C1–C2 couple (∼1.55 Å separation). On the diagonal, we observe several intense peaks at 288.1 eV, associated with a C2-π^*(C2), as well as 295.7 eV (two overlapping transition) and 298.3 eV, associated with C1-σ^*(C1–F), C2-σ^*(C1–C2), and C1-σ^*(C1–C2) transitions, labeled 2, 2′, and 3, respectively. This subsystem produces spectra with the largest number of off-diagonal features. Here, all the transitions appear to be coupled, and it will not be straightforward, from an experimental point of view, to disentangle the origin of the different signals observed. For this reason, the identification of core-specific signals and core-core coupling terms in the simulations is a necessity to effectively identify the many—overlapping—spectral features.

The origin and the extent of the shift of the double excitations (giving rise to ESA signals) with respect to the single excitations are detailed in Sec. III D. Here, we simply discuss one further observation, namely, that there is a class of positive (bleach) cross-peaks that is never going to be canceled by the corresponding double excitations. These are cross-peaks coming from different—bright—excited states pertaining the same core orbital (i.e., peaks at Ω₁/Ω₃ = 295.7/288.1 eV and Ω₁/Ω₃ = 288.1/295.7 eV coupling transitions 1 and 2′, as well as peaks at Ω₁/Ω₃ = 295.7/298.3 eV and Ω₁/Ω₃ = 298.3/295.7 eV coupling transitions 2 and 3). The corresponding double excitation, in which a second electron is excited from the same core, does appear dramatically blue-shifted (by tens of electronvolt and hence not shown in the 2D maps in Fig. 4). Indeed, the energy required to extract a second electron from the same core-excited orbital is much higher than the one required to excite the first electron due to the reduced shielding of the positive nuclear charge that the removal of the first electron produced.

Finally, an additional source of signals in the 2D spectrum, not shown in Fig. 4, is worth mentioning. This is due to the possibility of de-exciting a core excited system (i.e., a SE processes) to a valence excited state of nπ^*, σπ^*, or σσ^* character instead of the electronic GS. This de-excitation process is the electronic counterpart of the de-excitation to a vibrationally “hot” GS state observed in optical spectroscopy. Our computations suggest that in ESCA, these signals would appear in a spectral region red-shifted by ∼10 eV with respect to the one documented in this work. We elaborate further on this point in Sec. IV.

The 2D X-ray spectra for the ground state bleaching and stimulated emission components of 2DCXS spectra are shown in Fig. 5, employing the same set of states already presented for the XANES spectra of Fig. 3. The same XANES spectra of Fig. 5 are also reported on the side of the maps to facilitate the identification of the various core-signals. As with the linear absorption spectra, the TDDFT (a) and RASSCF/RASPT2 (b) results show reasonably good agreement. Recall that the TDDFT spectrum has a slightly narrower energy range. As such, the corresponding 2D spectrum spans ∼1 eV less than the RASSCF/RASPT2 result for the same peaks. Similarly to the XANES, the main differences between the two maps come from the peaks corresponding to C3 and C4. As a result of the compression in the TDDFT spectrum, the C3 and C4 peaks are red-shifted by ∼0.4 eV, relative to the RASSCF/RASPT2 result. Despite this shift, however, the two levels of theory show good agreement, as shown by the same overall intensity, energy range, and carbon K-edge contributions in each peak of the two spectra.

The many cross-peaks observed in the GSB-SE maps should be superimposed with ESA signals in the total map, as already shown for the reduced system studied in Fig. 4. In particular, we expect that

• •

  Cross-peaks pertaining noninteracting transitions (e.g., transitions of distant cores) will disappear, canceled by corresponding ESA signals;

• •

  Cross-peaks pertaining interacting transitions will only be partially canceled or not canceled at all, according to the strength of the interaction which is revealed by the shift of the ESA signals with respect to the GSB peak position;

• •

  Cross-peaks pertaining transitions computed from the same core will never be canceled as the corresponding ESA signals are dramatically blue shifted.

As evident in the 2D correlation plots presented in Secs. III B and III C, one notable strength of 2D spectroscopy lies in its ability to resolve intercore couplings as a function of the core proximity. At this point, we should clarify what we mean by coupling between transitions. The coupling we refer to is what in other contexts (as, e.g., in the description of Frenkel excitons) is known as quartic coupling.⁹⁷ This implies that double excitations are not additive, i.e., exciting a given core C_i in the presence of an already excited core C_j is not the same as exciting it when C_j is not excited. Quartic couplings shift the excitation energies and affect the transition dipole moments. Indeed, we observe very large quartic couplings which can reach up to 10 eV red-shift. Instead, we observe only a very weak quadratic coupling between core-transitions, which allows us to split (i) single-core computations in four different core-specific calculations (obtaining the total spectrum as a sum of the four subspectra) and (ii) double-core computations over couple of cores (e.g., C1–C2, C1–C3, and C1–C4). We now describe the physical background of the spectral shifts.

The quartic coupling is associated with orbital relaxations following a core-excitation: the dramatic decrease in the electronic shielding increases the electronegativity of the core-excited center. This lowers the energy of all σ- and π-type orbitals (i.e., non-diffuse orbitals) pertaining this center. Due to the electronic stabilization, we observe decoupling of the aforementioned orbitals from the manifold of generally delocalized σ- and π-type molecular orbitals and formation of localized orbitals around the core-hole. These relaxation events have a three-fold effect. First, due to the electronic stabilization which affects also the core orbital, a subsequent transition from the same core results in a strong blue-shift of the corresponding ESA (∼50 eV in the case of carbon). Second, subsequent transitions from any other core to the stabilized orbitals experience a significant red-shift with respect to core-excitations from the ground state. The shift magnitude is determined by the separation of the first core-hole/electron pair. In particular, the more diffuse is the orbital in which the first core-electron is excited, the stronger is the red-shift of the subsequent core-excitation which can reach values up to 10 eV (see Fig. 2). This is made apparent by looking at the shifts reported in Table II for 1s(C4) → σ^*(C1–F) transitions following a C1 core-excitation to virtual orbitals of Rydberg, π^* and σ^* nature. Third, due to the localization of the virtual orbitals, subsequent core-transitions are bright only for neighboring atoms sharing a localized orbital (e.g., a σ^* orbital between two adjacent carbons). This is the reason for the absence of off-diagonal features (i.e., absence of couplings) in the 2D map of the C1–C4 couple: transitions with the largest spectral shifts [e.g., 1s(C4) → σ^*(C1–F)]involve nonoverlapping starting and arrival orbitals due to the large spatial separation of the involved carbons.

We have simulated linear and nonlinear X-ray spectroscopy of the ESCA molecule at the carbon K-edge, both at the TDDFT and at the RASSCF/RASPT2 level of theory. The two levels of theory produce rather similar results by comparing ESCA XANES spectra and 2DCXS GSB and SE contributions. Interestingly, XANES signals from different cores do overlap in the considered ∼12 eV window, producing a structured absorption spectrum in which core-specific transitions are not distinct. This is different from XPS, where the carbon core transitions are clearly distinct.

Furthermore, we showed 2DCXS maps, computed at the RASSCF/RASPT2 level, on the top of systems with reduced dimensionality (couples of carbon centers), with a restricted active space and number of states considered. Here, the effect of considering carbons at an increasing spatial distance has been explored. The physical origin behind the observed energetic shift of the ESA was related to the spatial overlap of the starting core and arrival orbitals, which is different for π^*, σ^*, and Rydberg orbitals, thereby inducing spectral shifts of varying magnitude.

Based on our findings, we envisage a number of interesting future experiments that may be performed employing different colors of pump and probe X-ray pulses.

• •

  A single-color 2DCXS may be used to detect couplings between transitions on different cores. We expect that analysis of the network of couplings and rationalization of their origin may eventually allow us to extract structural information from these kinds of maps in a way similar to 2D-NMR spectroscopy.

• •

  A two-color 2DCXS pumping at the carbon K-edge and probing at ∼50 eV higher energies may be used to detect the strongly blue-shifted ESA transitions associated with double-core excitations pertaining the same core. In this way, one can profit from a bleach-free and thus less congested signal in reconstructing the geometry.

• •

  A two-color 2DCXS pumping at the carbon K-edge and probing at ∼10 eV lower energies may be used to detect de-excitation processes from singly excited cores to valence excited states. These transitions effectively fill the core-hole with a valence electron, thus resolving the electronic structure of the manifold of occupied valence orbitals. Again, the oscillator strength is determined by the spatial overlap of the valence and core orbitals.

Our results illustrate the potential of 2DCXS to become the X-ray counterpart of 2D-NMR, with its site and structural sensitivity, and show that theoretical computations are essential for the interpretation of the congested 2D maps.

The enlargement of the active space, together with the evaluation of all the states that may contribute to the explored spectral window, and the consideration of all the possible couples of cores are necessary steps that need to be taken to move from a simplified ESCA model (perfectly suited for pedagogical purposes) to a realistic simulation of the complete four-centers 2DCXS spectra.

M.G., A.N., A.B., F.S., and S.M. acknowledge support from the Chemical Sciences, Geosciences, and Bio-sciences division, Office of Basic Energy Sciences, Office of Science, U.S., Department of Energy, through Award No. DE-SC0019484.

The SOS expressions, at t₂ = 0, for the rephasing (K_I) and nonrephasing (K_II) GSB, SE, and ESA contributions to 2DCXS maps are reported hereafter, together with their respective double-side Feynman diagrams⁹⁷ (see Fig. 6)

where ξ_ab = ω_ab − i/τ, τ being the state lifetime, assumed to be the same for all the considered states.

All the pulses were assumed to have the same polarization and to be δ-like in time. An isotropic ensemble of ESCA molecules was assumed, and an orientational averaging was carried out in the simulation of the spectra.