Towards an ab initio theory for metal L-edge soft X-ray spectroscopy of molecular aggregates

The Frenkel exciton model was adapted to describe X-ray absorption and resonant inelastic scattering spectra of polynuclear transition metal complexes by means of the restricted active space self-consistent field method. The proposed approach allows to substantially decrease the requirements on computational resources if compared to a full supermolecular quantum chemical treatment. This holds true, in particular, in cases where the dipole approximation to the electronic transition charge density can be applied. The computational protocol was applied to the calculation of X-ray spectra of the hemin complex, which forms dimers in aqueous solution. The aggregation effects were found to be comparable to the spectral alterations due to the replacement of the axial ligand by solvent molecules.


I. INTRODUCTION
Soft X-ray L-edge spectroscopy has become a standard technique to investigate the intricate details of the electronic structure of transition metal compounds. The most popular variants encompass the X-ray Absorption (XAS) and Resonant Inelastic X-ray Scattering (RIXS) spectroscopies, allowing to address the properties of both unoccupied and occupied valence molecular orbitals. 1,2 However, the accurate theoretical prediction of L-edge core-and valence-excited electronic states of transition metal compounds often requires taking into account multiconfiguration and spin-orbit coupling (SOC) effects. Here, the combination of the Restricted Active Space Self-Consistent Field (RASSCF) method 3 and the atomic mean-field integral approximation 4 has been proven to be a versatile computational tool. [5][6][7][8][9][10][11][12][13][14] For metal complexes, where the spectroscopically active region is rather localized, the prediction of L-edge spectra requires an active space including the 2p and all orbitals with notable metal 3d-contributions. Such a choice corresponds to account for the most important correlation terms as well as dipole allowed transitions. Further, it allows to keep the active space quite compact and the number of considered electronic states in the order of hundreds or a few thousands.
However, the treatment of systems with multiple metal centers goes beyond the current numerical capabilities because of the fast growth of the number of configurations with the size of the active space. Moreover, in this case the number of relevant electronic states scales as tens or hundreds of thousands what hinders the theoretical interpretation of the structure of multicenter systems such as molecular aggregates. To cope with such situations, in the present article, we adopted a strategy known from the theory of excitation energy transfer in molecular aggregates. 15,16 Here, the total system is decomposed into its constituent monomers such that local electronic excitations can be clearly defined, and the computational effort is substantially a) Electronic mail: sergey.bokarev@uni-rostock.de 2329-7778/2016/3(6)/062601/11 V C Author(s) 2016 3, 062601-1 reduced. Such an approach will be particularly justified in cases where the monomers forming the aggregate are held together by van der Waals forces. Considering core-hole excitations, which are rather localized at a particular metal atom, such a separation strategy might even be justified for multiple metal centers within a covalently bound complex.
In the present proof of concept study, the versatility of this exciton coupling approach is exemplified for the hemin molecule (see Fig. 1). In addition to its high biological relevance as a constituent of the hemoglobin active center, the hemin complex is interesting due to its aggregation behavior in different solvents. Staying monomeric in polar solvents like ethanol or dimethyl sulfoxide (DMSO), it forms dimers in water solution. 17 The effect of aggregation was recently addressed by means of soft X-ray Fe L-edge absorption spectroscopy in transmission (XAS) and partial fluorescence (PFY) modes as well as by off-resonant X-ray emission (XES) and RIXS on the examples of DMSO and aqueous solutions. 18,19 The general shape of the spectra for both cases was quite similar, and the pronounced difference in broadenings for RIXS as well as a 1.3 eV energy shift in off-resonant XES was taken as an indication of aggregation. These effects were, however, solely attributed to the pp-stacking, 19 although ligand p-orbitals are barely influencing local Fe 2p ! 3d transitions measured with L-edge spectroscopy. Although very recently it was shown that K-edge absorption and emission spectra show features which can be associated with p-type interactions, 20 the interaction upon electronic excitation could be, in general, more complex due to the resonant coupling of essentially degenerate electronic states of the two monomers. 15 Here, we present a first principles approach capable of quantifying the various coupling-induced effects such as energetic shifts and redistribution of oscillator strengths. In addition, the labile equilibrium between different species in a solution could be a source of spectral changes. That is why in the present paper, we have also studied the influence of the axial ligand on the X-ray spectra of the monomeric hemin. The paper is organized as follows. In Section II, we present the general theory for the exciton coupling model as adapted to the calculation of first-and second-order X-ray spectra. For the sake of simplicity, we restrict ourselves to the case of the dipole approximation. Further, the one-and two-particle approximations (TPAs) will be introduced as a means for restricting the accessible excitation space. The computational setup is detailed in Section III. In Section IV, we first discuss speciation effects for the monomeric spectra. Subsequently, dimer spectra are provided for different orientations of the two hemin monomers. Section V gives summary and conclusions.

II. THEORY
In the following, we will consider a molecular aggregate and label the monomers by M. Let us denote the electronic states of monomer M by jA M i with the ground state being jA M i ¼ jg M i. Separating the total Hamiltonian according to these states are solutions of the monomeric Schr€ odinger equations In Eq. (1), V MN is the matrix element of the Coulomb operator between the monomer states. Assuming that there is no wave function overlap between local excitations of the two monomers, the aggregate wave function can be written in a direct (Hartree) product form. This yields the matrix representation of the aggregate Hamiltonian (assuming frozen nuclei) 15 Here, the Coulomb matrix elements J MN ðA M B N ; C N D M Þ have been introduced, which can be expressed in terms of generalized monomeric total charge densities, N A M ;D M ðrÞ Indices A, B, C, and D denote different electronic states of monomers M and N. Provided that the separation between the monomers, e.g., center to center distance, jX MN j, is large as compared to the extension of the transition densities and that the coupling to charge densities can be neglected (e.g., for charge neutral systems), the transition dipole approximation can be invoked, which gives 15 where d A M D M ¼ hA M jdjD M i are the transition dipole matrix elements, with d being the electronic dipole operator. Note that as long as metal-centered transitions are considered, the transition densities are rather localized and the dipole approximation should be valid despite the close proximity of the Heme B planes. Let us specify the situation to that of a dimer (M ¼ 1, 2) as shown in Fig. 2. The different monomeric state manifolds will be denoted as ground jg M i, valence-excited, jv M i, and coreexcited, jc M i, states. Equation (3) contains couplings between all possible transitions in the dimer system. In the following, we will make use of the fact that the actual processes of interest, namely, XAS and RIXS, are of first and second order, respectively, with respect to the interaction with the external field. This suggests employing either a one-or two-particle basis. In the former, states of the type ja 1 g 2 i and jg 1 a 2 i (a ¼ v, c) are incorporated, while the latter includes, in addition, states of type ja 1 ; b 2 i (a, b ¼ v, c) and is, in principle, exact for the dimer. Note that this effectively corresponds to a CI-doubles-like treatment of the composite system with X-ray specific preselection of configurations. The respective Hamiltonian matrix is readily calculated in terms of the monomeric excitation energies and the Coulomb integrals in Eq. (3). To make the calculation computationally feasible, different approximations have been applied (cf. Fig. 2) to reduce the size of the Hamiltonian matrix as well as the number of terms in the RIXS expression (see below): (a) couplings between g $ v transitions (e.g., J 12 (v 1 g 2 ,v 2 g 1 )) and between g $ v and v $ c transitions (e.g., J 12 (v 1 c 2 ,v 2 g 1 )) were neglected, i.e., set to zero in Eq. (3), (b) no coupling between static dipoles as well as between dipoles and transition dipoles have been included, (c) a pre-selection of core-excited states to construct the basis functions has been applied according to an energy window with a width of 65r (r is the width of the Gaussian excitation pulse) around the center of the excitation pulse.
Note that in contrast to common Frenkel exciton theory not only resonant couplings have been included but we also partially account for induction and dispersion effects. 21 The use of FIG. 2. Different choices of excitonic bases and the corresponding transitions included in the dimer coupling. For the oneparticle basis (a), only transitions from or to the ground state have been included (one-particle approximation: OPA). For the two-particle basis (b), de-excitations from an arbitrary core-excited state to any other state are allowed (two-particle approximation: TPA). the one-and two-particle basis together with these conditions will be called one-(OPA) and two-particle approximation (TPA), respectively.
Upon diagonalization of the resulting Hamiltonian matrix, one obtains eigenstates (jii; jni, and jf i for initial, intermediate, and final states, respectively) and the respective transition dipole moments (d ni and d fn ) from which XAS and RIXS spectra can be calculated. Here, w(E i ) is the Boltzmann weight of the corresponding initial state jii, and e 1 (e 2 ) is the polarization of the incoming (outgoing) light. The inner sum in Eq. (7) corresponds to the matrix element of the electronic polarizability with respect to the initial and final states. For the calculation of Partial Fluorescence Yield (PFY) spectra, the 2D-RIXS spectra have been summed up over the emission energy E em , yielding an integrated amplitude over all possible outgoing photon energies in the emission energy range [E min : E max ] corresponding to a fixed excitation E exc Note that polarization effects in the XAS and RIXS spectra were taken into account as described in Ref. 22, assuming a free tumbling of the molecules in solution under the condition that the polarization is detected orthogonal to the incoming beam polarization in the laboratory frame.

III. COMPUTATIONAL DETAILS
The energies and transition dipoles for the monomeric Heme B derivatives have been calculated using geometries, which were optimized with the density functional theory (BLYP functional and the LANL2DZ basis set for iron and a 6-311 þ G(d) basis set for all other elements) using Gaussian 09. 23 Subsequently, the electronic wave function has been determined via a RASSCF calculation as implemented in the Molcas 8.0 program package 24 [25][26][27] The active space for all Heme B derivatives and for all spin configurations has been chosen to consist of three 2p orbitals in RAS1 and the five 3d orbitals in RAS2 to describe dipole allowed 2p ! 3d electronic transitions. Further, one hole was allowed in RAS1 and 5 or 6 electrons in RAS2, so that the total number of active electrons is 13 (cf. Fig. 1).
The spin-orbit coupling (SOC) was included within the LS-coupling scheme and using the atomic mean-field integral approximation as implemented in the Molcas 8.0 program package. 4,24 At this point, the calculation has been restricted to sextet (S ¼ 5/2) corresponding to the ground state spin and quartet (S ¼ 3/2) spin configurations according to the spin-orbit selection rule DS ¼ 0, 61. This resulted in 16 sextet and 174 quartet spin-free states or 792 spin-orbit coupled states, with 102 valence states covering energy range up to 7.3 eV and 690 core states with energies between 710.6 and 732.9 eV. To cover this part of the spectrum within a supermolecule-type dimer calculation, it would be required to include 61 and 34 265 sextet and quartet spin-free states, respectively, or 137 426 SOC states. This goes far beyond the present computational capabilities. Scalar relativistic effects have been taken into account using the Douglas-Kroll-Hess transformation. 28,29 Based on the local electronic wave functions, dimer states have been obtained as outlined in Sec. II. In the results presented below, we will consider the cases of one-and two-particle approximations (OPA and TPA) for the sake of comparison. The distance of 7 a.u. between monomers was chosen as representative on the basis of molecular dynamics simulations 17 of this very non-rigid system, with heme ligand planes being parallel to each other and axial ligands pointing to the outside of the dimer. To study the effect of mutual orientation, three geometrical configurations have been considered: 0 rotated (COOH groups of both monomers are on one side of the dimer), 90 rotated, and 180 rotated (COOH groups are on the opposite sides).
Spectra have been calculated according to Eqs. (6)- (8). The broadening is determined by the finite life times of the intermediate states as well as through a Gaussian lineshape of the excitation pulse. In order to conform typical experimental conditions, the latter has been chosen to have a width of 0.25 eV. Since even within the L 3 -band the lifetimes may vary, 30 the C n parameters in Eq. (7) for both monomer and dimer have been assumed to be 0.09 eV for states below 709.2 eV, 0.26 eV between 709.2 and 711.6 eV, 0.43 eV between 711.6 and 719.2 eV, and 0.61 eV above 719.2 eV as a best fit to the experimental data. 18,19 The hemin monomer has a sextet ground state, where the degeneracy of the six components is slightly lifted due to SOC and thus the states are split into three Kramers doublets. For the monomer spectra calculations in Section IV A, Eqs. (6) and (7), these six ground states are populated according to a Boltzmann distribution at a temperature of 300 K. For the case of the dimer in Section IV B, 36 combinations of the ground states are possible; here, for simplicity we considered only first four of them resulting from the lowest Kramers doublet of both monomers. This effectively corresponds to very low temperatures. For the comparison, the same Kramers pair was used for the monomer calculations as well.
For the PFY spectra, the energy interval [E min : E max ] has been chosen to range from 695 eV to 735 eV, corresponding to the 3d ! 2p radiative decay channel. An absolute energy shift of À2.8 eV has been applied to all theoretical spectra to ease the comparison with experiments.

A. Monomer spectra
Before addressing effects of the excitonic coupling in the hemin dimer, the change of the XAS, PFY, and RIXS spectra due to the different solvents need to be quantified for the case of the monomer. The following discussion in this section is based on the assumption that there are different species in solution resulting from an exchange of the chloride ligand with solvent molecules. Specifically, spectra have been calculated for various axial ligands: Cl -(original hemin), H 2 O, and DMSO, which have been used as solvents for experimental measurements, 18,19 and OHthat can be formed in an aqueous solution due to hydrolysis (Fe 3þ þ H 2 O ! FeOH 2þ þ H þ ). The XAS and PFY spectra of compounds with these ligands can be found in Fig. 3.
The XAS (Fig. 3(a)) shows a distinct sensitivity in the L 3 -edge with respect to the ligand's nature. For [Heme B-H 2 O] þ a splitting of this peak in two components at 710.0 eV and 710.9 eV can be observed and for [Heme B-DMSO] þ the lower energy feature appears as a shoulder at 710.1 eV. In contrast, for [Heme B-Cl] 0 and [Heme B-OH] 0 , there is a single main L 3 -peak only. In comparison to XAS, the PFY spectra ( Fig. 3(b)) are less sensitive to the nature of the ligand but show an intensity enhancement of all features for excitation energies above 716 eV. The differences between XAS and PFY are caused by inelastic features that become more pronounced for higher excitation energies. An explanation for this fact will be given after the interpretation of 1D-RIXS spectra below (see Fig. 4).
Due to the multi-configurational nature of the core-excited states, the excited 2p electron is mostly evenly distributed over the 3d orbitals. However, the most prominent transitions in XAS are due to 2p ! d x 2 Ày 2 and d z 2 excitations. Interestingly, the splitting in XAS of [Heme B-H 2 O] þ is due to the energetic lowering of the 2p ! d z 2 transitions with respect to the 2p ! d x 2 Ày 2 ones upon change of axial ligand from Clto H 2 O.
2D-RIXS spectra have been obtained for all Heme B derivatives. Exemplarily, Fig. 4(a) shows the 2D-RIXS spectrum for the [Heme B-Cl] 0 case. However, since the analysis of 2D spectra is rather difficult and the differences between the species are not very pronounced in the 2D presentations, one-dimensional cuts of the RIXS spectra will be analyzed for [Heme B-Cl] 0 and [Heme B-H 2 O] þ . Five excitation energies, belonging to distinct spectroscopic features, were selected as shown in the left panel of Fig. 4(b).
The 1D-RIXS spectra show a prominent elastic peak for lower excitation energies, whose intensity decreases upon increasing the excitation energy. This behavior can be rationalized as follows. Below 712 eV core-excited sextet states are dominating the spectrum. Due to the spin selection rules, the preferred emission is to the sextet ground state, thus yielding an intense elastic peak. For excitation energies above 712 eV, the core-excited states are mostly of quartet type. Here, the most intense emission is the relaxation from a core-excited state to a valence-excited state with a high quartet contribution. The elastic peak corresponds to a relaxation to the (sextet) ground state, which is spin forbidden and therefore less intense than the inelastic features.
The 1D-RIXS spectra for [Heme B-Cl] 0 and [Heme B-H 2 O] þ differ mostly in the inelastic peaks, whereas the elastic peak has a comparable intensity in both spectra. Among the inelastic features, the peaks with loss energies E em À E exc of À2.5 eV, À3.6 eV, À4.7 eV, and À7.0 eV are most prominent and labeled by 1-4 in Fig. 4. All inelastic features in the RIXS spectra are

B. Dimer spectra
Below we will focus mainly on the L 3 -edge since it is more structured and less subject to the lifetime broadening than the L 2 -edge. The OPA does not show any notable difference between monomer and dimer spectra both for XAS and RIXS for all three orientations of the molecules. The reason for the minute differences is that the coupling of transition densities is rather weak for core excitations as compared to valence ones. This can be attributed to the intensities of the metal 2p ! 3d and 3d ! 3d transitions relevant for L-edge X-ray spectra that are lower (due to smaller radial overlap and dipole selection rules) than those of the p ! p * and n ! p * transitions usually discussed in the case of organic dyes. Moreover, the OPA exciton basis by construction should be appropriate only for the first-order XAS spectra (one-photon transitions). For the second-order RIXS, it is natural to take the ja 1 c 2 i (a 1 ¼ g 1 , v 1 ) type of basis functions into account, since one needs to describe the two-photon g ! c ! a transitions which are additionally interfering with each other (see Eq. (7)). Indeed for RIXS spectra, in contrast to the OPA exciton basis, the TPA basis predicts aggregation effects that are up to one order of magnitude larger (not shown). It should be noted, however, that even for TPA the form of XAS essentially does not change upon dimerization. Distinct fingerprints of dimerization can be seen in RIXS, evidencing that absorption spectroscopy should be less sensitive to aggregation than RIXS. Therefore, in the following, we discuss only TPA results for RIXS spectra.
As mentioned in Section III, to reduce the size of the Hamiltonian matrix (Eq. (3)) and the number of terms in the innermost sum in Eq. (7), we have limited the basis to those states jci which are within a 61.25 eV energy window around the prominent absorption features. This allowed to reduce the rank of Hamiltonian block to be diagonalized from 140 964 to about 22 000 (depending on the excitation energy). Such an approximation is justified by the finite width of the excitation pulse as well as the fast Lorentzian decay of interference terms in Eq. (7) with the energy separation between radiative channels. Further, recall that only four degenerate initial states were taken into account for the dimer case. For the purposes of comparison, the corresponding monomer spectra include only two degenerate initial states. Note that as far as the monomer is concerned, the differences between two and six initial states are small as compared to the dimerization effects.
The 1D-RIXS spectra in Fig. 5 show notable differences between the monomer and dimer cases: a prominent increase of intensities of elastic bands (especially for the 180 orientation) is observed, whereas the energy shifts are not larger than 0.1 eV. The spectra demonstrate a distinct orientation dependence, with the largest differences from the monomer case being observed for 0 and the smallest for 180 . For 0 , the overall intensity mostly increases, for 90 slightly decreases, and for 180 stays almost intact. However, the high density of states hinders a detailed analysis in terms of, e.g., orientations of transition dipole moments. Qualitatively, one can say, that apart from the elastic features, the largest deviations are observed for d-orbitals having an out-of-plane zcomponent (d xz , d yz , and d z 2 ). The "in-plane" transitions are affected only for the 0 orientation.
To summarize, although the L-edge transition moments for the monomers are quite small if compared to valence excitations of organic dyes, there is an effect of dimerization that can be seen in the RIXS spectra. Moreover, there appears to be a pronounced dependence on conformation. Given the fact that in a solution the dimer system is floppy and interconversion between conformations is rather likely, 17 a direct comparison with experiment would require, e.g., a computationally demanding molecular dynamics based sampling of spectra.

V. CONCLUSIONS
The Frenkel exciton model is usually applied to describe aggregation effects on spectra in the visible domain as well as excitation energy transfer in molecular aggregates. 16 In the present contribution, its basic idea has been adopted to the computation of the core-level spectra of multi-center transition metal compounds. Thereby, the RASSCF-based protocol for treating multiconfiguration and spin-orbit coupling effects has been extended to multi-center systems, which are not accessible by the standard protocol due to computational limitations. While in molecular aggregates only valence excitations are of relevance and often the treatment can be reduced to monomeric two level systems, the description of core-level spectra of transition metals requires taking into account a large number of possible transitions. This renders the interpretation, e.g., in terms of a few coupled transition dipoles to become essentially impossible. Further, in the standard exciton theory, one usually classifies the collective excitations as zero-, one-, two-exciton states, etc. This is particularly useful in the context of (non)linear spectroscopy. 31 In the present case, however, such a classification is not very useful due to the multitude of possible excitations. Still, different approximations derive from the used aggregate basis. Here, we discussed the oneand two-particle basis, the latter being exact for the dimer, but an approximation for larger aggregates. It turned out that, similar to standard exciton theory, the one-particle basis is suitable for describing one-photon processes like XAS only. The two-photon RIXS requires taking into account two-particle excitations.
In general, due to the quite small transition dipole moments for the 2p ! 3d excitations, if compared to the valence p ! p * and n ! p * transitions of organic dyes, the effect of aggregation on XAS spectra will be rather small. As far as the RIXS spectrum is concerned, Eq. (7) points to a dependence on the electronic polarizability, making it more sensitive to the fine peculiarities of electronic structure. Here, even weak intermonomer contributions will determine the aggregation effect on the spectrum.
The developed protocol was applied to the hemin system forming dimers in water solution, while staying monomeric in other polar solvents. Remarkably, ligand coordination in various solvents has been shown to have a pronounced influence for XAS and, in particular, RIXS spectra. This is an important result, since solvent effects have not been considered in the previous experimental studies of hemin X-ray spectra. [18][19][20] In the present work, it was found that coordination and aggregation effects on RIXS spectra are of similar magnitude. This could make the unequivocal assignment of aggregation induced features difficult. However, in this proof-of-principle study, a direct comparison with experiment was not attempted. First, due to construction of the active space the pp stacking effect, discussed, e.g., in Ref. 19, was not taken into account. Second, it was found that the aggregation induced changes in the RIXS spectrum are depending on the mutual orientation of the monomers in the dimer. Due to the flexibility of the dimer structure in solution, a more accurate description would require a rather time consuming sampling of different conformations, i.e., by combining molecular dynamics with the present RIXS calculation.
The present Frenkel exciton-like approach to the X-ray spectroscopy of multi-center systems should be particularly suitable to describe coupled highly local core excitations of weakly bound van der Waals complexes. To include situations with more extended electron densities, the dipole approximation has to be replaced by a more accurate calculation of transition densities. This could be achieved using standard tools for integration of Gaussian or Slater-type orbitals. For situations where covalent bond formation is of importance or where the monomers are ferromagnetically and antiferromagnetically coupled, taking into account the exchange contribution would be mandatory. Finally, the present approach has been developed for a molecular dimer only. However, the extension to larger aggregates is straightforward although eventually bound to the applicability of a certain n-particle basis, with n being small enough to accommodate current computational resources.
The Frenkel exciton approach has been extensively used in the context of nonlinear spectroscopy and dynamics of dye aggregates. The present adaption to the X-ray regime, in principle, facilitates similar investigation for core-level excitations. Thus, upcoming ultrafast spectroscopic techniques in the X-ray regime 32,33 could be a target for future advancement of the present approach.