Explaining the $MoVO4−$ photoelectron spectrum: Rationalization of geometric and electronic structure

A multitude of studies over many decades has revealed that catalytic activity at transition metal oxide surfaces occurs at edge and point defects resulting from oxygen vacancies.^1,2 Transition metal oxides are effective redox promoters as they can readily cycle through two or more oxidation states. This characteristic is often enhanced in bimetallic systems, facilitating high degrees of reactivity and selectivity.^3–7 However, a detailed experimental study of such systems in situ is complicated by the small ratio of defect sites to stoichiometric solid resulting in low signal-to-noise ratios. Furthermore, surfaces undergo continuous structural changes such as defect formation, metal reduction, metal oxidation, reagent adsorption, reagent desorption, surface restructuring, and subsurface oxygen diffusion.^8–10 Defect sites have incomplete valencies and can exhibit strained electronic structures and atypical metal oxidation states. Thus the study of transition metal oxide clusters is largely motivated by their utility as model systems for defects on catalytic surfaces.¹¹ 

Photoelectron spectroscopy has proven to be a powerful tool for characterizing the electronic and molecular structures of metal oxide clusters. Detailed insight into ground and low-lying electronic states of the neutral species can be obtained from the resulting spectra. Observed vibrational excitations reflect the structural changes associated with changing the electron population in the orbital associated with electron detachment. However, analysis is often complicated by the number of possible spin states and structural isomers/conformers. Furthermore, as spectral congestion grows with system size, it becomes increasingly difficult to differentiate fundamental, overtone, and combination bands. Theory-based spectral simulations are a critical tool to aid assignments of species by calculation of electronic states and nuclear configurations, along with associated vibrational modes and transition intensities required for spectral simulation.^12–15 

The application of theoretical methods to systems such as metal oxide clusters is further complicated by orbital degeneracies which result in multiconfigurational electronic structure and are intrinsically connected with the observed catalytic activity. Therefore, modeling of spectra formally requires a multi-reference approach; however, application of such methods is typically impractical owing to steep computational scaling. Despite the issues in obtaining a correct description of electronic structure, a large number of studies have employed hybrid Density Functional Theory (DFT) to calculate low energy minima, Adiabatic Detachment Energies (ADEs), Vertical Detatchment Energies (VDEs), and vibrational modes in metal oxide clusters.^16–25 Development of efficient approaches for a balanced description of dynamic and static correlation remains an active research area.^26,27

An example of profound failure of DFT to appropriately describe a metal oxide cluster based on photoelectron spectral simulation is $MoVO4−$⁠. Definitive assignment of the anion and resulting neutral structures and electronic states was not possible because of an absolute lack of agreement between computation-based simulations and the well-resolved experimental spectrum, which displays a nearly vertical transition with several short but well-resolved vibrational progressions.¹⁶ Previous work has indicated that Yamaguchi’s Approximate Projection (AP) method^28–39 can correct for spin contamination errors in metal oxide clusters.²⁸ However, the relevant states of $MoVO4−$ displayed low levels of spin contamination. Correct spectral simulation is still not achieved due to the presence of close lying electronic states, low frequency vibrational modes, and geometric symmetry breaking.

In this investigation, we seek to resolve the spectral simulation issues. We therefore attempt to definitely assign the species involved in the observed photodetachment spectrum. In doing so, we reveal the behavior of transition metal oxide clusters to form neutral states accessed by two-electron transitions^40–43 and are able to rationalize the observed behavior based on a simple molecular orbital picture.

In order to examine the $MoVO4−$ potential energy surface, the same procedure was followed as our recent study on $MoNbO2−$⁠.²⁸ Following previous work,^16–20 isomers are labeled and denoted by dividing the molecule into basins, with the boundary of each basin defined by a line running through each metal center perpendicular to the metal-metal bond. In a bimetallic system, there are three basins and each basin contains oxygen atoms labeled depending on whether they are at the molybdenum terminus (X), bridging the two metals (Y), or at the vanadium terminus (Z). The nomenclature for each structure then is XYZ, where each letter gives the number of oxygens in each basin. For example, 211 indicates that there are two terminal oxygens bound to molybdenum, one bridging oxygen, and one terminal oxygen bound to vanadium. From each isomer, all possible configurations were obtained by considering the highest symmetry point group and all subgroups. Energies of each structure were then obtained for singlet, triplet, quintet, and septet spin states for the anion cluster and doublet, quartet, sextet, and octet spin states for the neutral cluster, generating a total of 432 initial geometries.

All calculations were carried out using a local development version of Gaussian.⁴⁴ The employed model chemistry was B3PW91 with the SDDPlusTZ basis set, which is comprised of the SDD relativistic pseudopotentials on the molybdenum and vanadium centers to describe core electrons, an augmented Stuttgart/Dresden (SDD) basis set for valence metal orbitals, and aug-cc-pVTZ to describe orbitals on the oxygen atoms. Analysis of the initial and final Kohn-Sham determinants was carried out to ensure that they corresponded to a stable solution.⁴⁵ Harmonic vibrational modes were computed at the same level of theory used for geometry optimization, and these modes were used in the simulation of Franck-Condon (FC) profiles.^46,47 The spectra were simulated at 100 K using the FC approximation with Dushinsky rotations.^13–15 Gaussian widths were added to the stick spectra with widths of 0.0062 eV (50 cm⁻¹) and 0.0124 eV (100 cm⁻¹).

In order to reproduce the observed photoelectron spectrum, the models used in the simulations must produce accurate relative spin-state energies, geometries, and force constants. Thus the discussion below is broken into three main sections. We first discuss the identification of the global ground state, which allows us to identify available spin-channels in the photodetachment process. From this, we are able to compute the VDE and ADE resulting from the associated transitions. Second, we discuss the issue of very low energy barriers between different conformational isomers in both $MoVO4$ anion and neutral states and how the lack of a well defined minimum energy structure complicates the calculation of FC factors. Simulation of vibrational progressions accounting for the presence of multiple low energy barriers demonstrates assignment validity. Finally, we discuss the presence of multiple broken-symmetry electronic structure solutions that suggest the presence of close lying states and use a Molecular Orbital (MO) picture along with Natural Ionization Orbital (NIO) analysis⁴⁸ to determine the relevant states for comparison with experiment.

In previous work by Mann et al., it was found that two anion structures have nearly identical energies and that either could be the global anion minimum.¹⁶ These two structures were the 220 ⁵B₁ state and the 121 ³A′ state which they calculated to be only 0.06 eV apart. Our extensive search (using a slightly different model chemistry than employed in the original study) did not yield any additional candidate structures. Mulliken spin density analysis and Natural Bond Orbital (NBO) analysis indicated that both unpaired electrons in the ³A′ structure were localized to the vanadium center occupying 3d-like orbitals. An additional electronic state was identified for the 121 structure, ³A″, in which both unpaired electrons occupy orbitals localized on molybdenum. However, this heavily spin contaminated state can be dismissed as irrelevant to this study due to its significantly high energy (+0.558 eV relative to 220 ⁵B₁).

In contrast to previous results using the B3LYP functional, our calculations using B3PW91 predicted the 121 ³A′ state to be 0.271 eV higher in energy than the 220 ⁵B₁ state (Table I). Previous work on $FeO2$ has indicated the average deviation of B3PW91 energies from multi-reference configuration interaction (MRCI) calculations to be on the order of 0.2 eV.⁴⁹ Assuming transferability of this result to $MoVO4−$⁠, the B3PW91 results indicate that 220 ⁵B₁ is the global minimum. While the 220 ⁵B₁ structure is nearly spin pure, the 121 ³A′ state exhibits mild spin contamination as determined by the value of $⟨S^2⟩$⁠. These low levels of spin contamination indicate that the broken symmetry density functional scheme is suitable, in contrast to other recent studies involving other transition binary metal oxide clusters.²⁸ As both B3PW91 and B3LYP studies found no other structures within 0.39 eV of the global minimum, we focus on 220 ⁵B₁ and 121 ³A′ states for the remainder of this discussion.

Further evidence that 220 ⁵B₁ is an unsuitable candidate for the initial state of the observed experimental spectrum can be obtained by comparing experimental and theoretical spectroscopic parameters, which indicate that transitions from the 220 ⁵B₁ state are not consistent with the observed spectrum. Table II shows the computed VDE and ADE for all spin channels accessible from the lowest energy ground states. The 220 neutral ground state is the ⁴B₁ state, and the ADE for the 220 ⁵B₁ → 220 ⁴B₁ + e⁻ transition (2.189 eV) is smaller than the experimentally determined value (2.66 eV). The strong peak intensity at the progression onset of the X and A bands in the experimental spectrum suggests that there is a small structural difference between anion and neutral geometries. The calculated relaxation energy on the neutral energy surface is an order of magnitude larger than the experiment (0.364 eV vs 0.06 eV), resulting in a much broader simulated FC profile than observed experimentally. The transition energy of the alternative spin channel, involving the 220 ⁵B₁ → 220 ⁶B₁ + e⁻ transition, is significantly higher.

In the case of the transitions initiated from 121 ³A′, two spin channels are energetically accessible: 121 ²A″ and 121 ⁴A″. Three different 121 ²A′ minima were found, one of which is not accessible from the anion through a one-electron transition, while two 121 ⁴A′ minima were located, one of which is not one-electron accessible. The lowest energy states of both the doublet and quartet are not accessible through one-electron transitions. Analysis of the electronic structure of the species involved in these transitions indicated that the relevant transitions are the 121 ³A′ → 121 ²A″ (1et) + e⁻ and 121 ³A′ → 121 ⁴A″ (2et) + e⁻ transitions (vide infra). In this section, we discuss the structures and properties of the clusters that must be accounted for in order to simulate the photodetachment spectrum before returning to the nature of the electronic states of these transitions.

Fig. 1 shows the structures involved in the lowest energy transitions. The 220 ⁵B₁ structure has C_2ν symmetry and is well defined assuming a single potential energy well around the minimum. The 121 ³A′ structure has C_s symmetry and, in contrast to the 220 ⁵B₁ structure, has flexible terminal oxygen atoms with low energy barriers to cis-trans isomerization across the plane formed by the two metal centers and two bridging oxygens. Mo–V–O are collinear and the motion associated with this angle bending is 104 cm⁻¹, while the V–Mo–O bend is out of the plane with frequency 50 cm⁻¹. The barrier for isomerization of the V–Mo–O bend is very small at just 0.42 kcal mol⁻¹ (146 cm⁻¹), well below the zero point energy.

The bending of the terminal oxygens to break C_2ν symmetry can be understood by considering the ligand environment of each metal. The oxygen atoms around Mo and V form an effective $MX3$ environment in a trigonal planar orientation. The ligand field splitting of the metal d-shell in such an environment gives rise to a degenerate d_xy and $dx2−y2$ pair being higher in energy than the degenerate d_xz and d_yz set and the $dz2$ orbital.^50,51 Depending on the nature of the ligand, the order of the d_xz and d_yz pair and the $dz2$ orbital will change. σ donating ligands favor the d_xz and d_yz pair while π donating ligands will give the $dz2$ orbital as the most stable in the d manifold. The presence of Jahn-Teller distortion, achieved by bending of the terminal oxygen out of the plane, will depend on the order of the d_xz/d_yz pair and the $dz2$ orbital. If the metal center is d¹, the symmetry distortion will occur when the d_xz and d_yz pair is lowest in energy; if the $dz2$ orbital is lowest in energy, only d² metals will undergo symmetry distortion. As the oxygen ligands are both σ (p_z) and π (p_x, p_y) donating, the exact order of these lowest energy d orbitals depends upon the particular geometry and the ratio of ligand electron donation types to the metal center.

The neutral structures associated with the transitions reflect the complexity observed in the anions from which they originate. The lowest energy vibrational modes in the 220 ⁴B₁ structure are 108 cm⁻¹ and 229 cm⁻¹ associated with wagging and bending of the two terminal oxygens. As with the 220 ⁵B₁ anion ground state, the 220 ⁴B₁ structure is described well by a single harmonic well. The three 121 doublet minima all show out-of-plane bending of both terminal oxygen atoms, including the Mo–V–O angle that is in the plane of the metals and bridging oxygens in the 121 ³A′ anion structure. In the 121 ²A″ (1et) structure, which was identified as the relevant state for comparison with experiment, the V–Mo–O bend has a barrier at 180° of 5.81 kcal mol⁻¹ (harmonic frequency of the associated vibrational mode is 39 cm⁻¹). While this barrier is relatively high, comparison with the 121 ³A′ potential energy surface (PES) shows that it lies along the valley connecting the two minima in the 121 ³A′ structure and so is likely to be obtained as the initial structure on the neutral state after photodetachment—the FC structure (Fig. 2, center). Thus 121 ²A″ (1et) is a double well system with two possible equivalent paths (Fig. 2, left). Similar double well systems have been identified and studied in tungsten oxide clusters where the presence of the double-well in the excited anion state renders FC based simulations inadequate and some form of conformational averaging must be used instead.¹⁹ A similar pattern is observed on the 121 ⁴A″ (2et) PES, although the two equivalent minima lie along the Mo–V–O bending coordinate orthogonal to the double-well coordinate in the 121 ²A″ (1et) structure (Fig. 2, right). The 121 ⁴A″ Mo–V–O coordinate has a barrier of 2.65 kcal mol⁻¹ with a vibrational frequency of 116 cm⁻¹. There is a very small barrier to inversion of the V–Mo–O angle of just 0.01(47) kcal mol⁻¹ to a trans minimum which lies 0.01(13) kcal mol⁻¹ above the cis minimum. Given that the energy difference is so small, inversion of this coordinate can be considered barrierless.

Given the low barriers for inversion, we considered the possibility that differences between the experimental and theoretical ADE and VDE values may stem from the fact that the minimum energy structures do not give a correct representation of the thermally accessible ensemble of structures. While the 220 ⁵B₁ → 220 ⁴B₁ + e⁻ simulation is relatively straightforward, as previously discussed, the presence of multiple/flat minima in the 121 ³A′ → 121 ⁴A″ (2et) + e⁻ and 121 ³A′ → 121 ²A″ (1et) + e⁻ transitions complicates the simulation of FC progressions as they assume the structures are described by a single harmonic well. In order to simulate these species, we follow the protocol used by Waller et al.¹⁹ described in Sec. II whereby we optimize to the critical point closest to the Boltzmann weighted average of the 121 ³A′ structures and ignore imaginary and very low frequency modes.

Figure 3 shows a comparison of the calculated FC spectrum with experiment for the three likely transitions identified above. The simulated spectra are shown with added Gaussian widths of 50 cm⁻¹ (blue) and 100 cm⁻¹ (red). The simulations were carried out at a temperature of 100 K. Clearly the 121 ³A′ → 121 ²A″ (1et) + e⁻ transition gives good agreement with the experimental A band. In particular, the profile shows a very intense peak at the onset consistent with experiment, the progressions have similar peak spacings, and the higher energy bands are broader due to an increase in the density of transitions. The 220 ⁵B₁ → 220 ⁴B₁ + e⁻ and 121 ³A′ → 121 ⁴A″ (2et) + e⁻ transitions are both broad with similar peak spacing in vibrational progressions and hence either can be considered to be in reasonable agreement with the X band. The fact that the 121 ³A′ → 121 ⁴A″ (2et) + e⁻ transition is a two-electron transition provides an excellent explanation of the high intensity difference between the A and X bands. Therefore, as the one-electron 220 ⁵B₁ → 220 ⁴B₁ + e⁻ transition should be of comparable intensity to the 121 ³A′ → 121 ²A″ + e⁻ transition, we use the relative intensities of the A and X bands to tentatively exclude the 220 ⁵B₁ → 220 ⁴B₁ + e⁻ transition, along with experimental evidence that both transitions arise from a common anion.¹⁶ Furthermore, the 220 ⁵B₁ → 220 ⁴B₁ + e⁻ transition is beyond the typical 0.3 eV error range for ADE and VDE energies.^18,52 Thus we assign the X band to the 121 ³A′ → 121 ⁴A″ (2et) + e⁻ transition.

To obtain a direct comparison with experimental results, the simulated 121 ³A′ → 121 ⁴A″ (2et) + e⁻ and 121 ³A′ → 121 ²A″ (1et) + e⁻ transitions were fitted to experiment minimizing the mean squared error of a combined spectrum with a spectral Gaussian linewidth of 50 cm⁻¹ (Fig. 4). The resulting spectrum was obtained with a mean squared error of 0.005. The 121 ³A′ → 121 ²A″ (1et) + e⁻ and 121 ³A′ → 121 ⁴A″ (2et) + e⁻ transitions were shifted by 0.316 and 0.294 eV and scaled by 0.812 and 0.525, respectively. The composite profile gives good agreement between the experimental and simulated progressions and lends weight to the transition assignments made above.

Again, three 121 ²A″ minima were located: (i) [B] 121 ²A″ (1et) with two unpaired electrons on vanadium and a single anti-parallel spin electron on molybdenum (Fig. 5); (ii) [C] 121 ²A″ (2et) with two parallel electrons on molybdenum and a single anti-parallel electron on vanadium (Fig. 5(c)); and (iii) [D] 121 ²A″ (1et2) with two spin paired electrons on molybdenum and a single electron on vanadium (Fig. 5(d)). [C] 121 ²A″ (2et) was found to be lowest in energy by 0.088 eV compared to [B] 121 ²A″ (1et), while [D] 121 ²A″ (1et2) was higher in energy by 0.135 eV. Attempts to find a further doublet state with α and β V electrons and an α Mo electron collapsed to the [D] 121 ²A″ (1et2) state and was thus discarded as being too high in energy. NIO analysis⁴⁸ indicated that transition from [A] 121 ³A′ to [B] 121²A″ (1et) and [D] 121 ²A″ (1et2) states results from one-electron transitions from molybdenum and vanadium d-orbitals, respectively. On the other hand, transition to the [C] 121 ²A″ (2et) state is accessed via a two-electron transition (Fig. 6).

There are two 121 ⁴A″ minima: (1) [E] 121 ⁴A″ (1et) which retains two unpaired electrons on vanadium and has a single parallel spin electron on molybdenum (Fig. 5(e)); and (2) [F] 121 ⁴A″ (2et) which has the opposite electronic configuration with two unpaired electrons on molybdenum and an unpaired parallel spin electron on vanadium (Fig. 5(f)). As with the doublet, NIO analysis confirmed that the transition to the [F] 121 ⁴A″ (2et) structure is a two-electron transition while the [E] 121 ⁴A″ (1et) structure results from a one-electron transition (Fig. 7).⁵³ As shown in Table II, [F] 121 ⁴A″ (2et) is lower in energy than [E] 121 ⁴A″ (1et) by 0.099 eV. We note that attempts to compute VDEs to the ⁴A″ (2et) state collapsed to the ⁴A″ (1et) electronic structure. Exploration of the PES located a transition state between the [E] 121 ⁴A″ (1et) and [F] 121 ⁴A″ (2et) structures corresponding to lateral motion of the bridging oxygens between the two metal centers only 0.025 eV above the [E] 121 ⁴A″ (1et) structure and 0.124 eV above the [F] 121 ⁴A″ (2et) structure. The transition structure that interconverts the two electronic states implies the existence of a crossing of the two states of interest. Furthermore, our calculations suggest that the intersystem crossing point is quite low in energy and we expect interconversion between these two states to occur relatively rapidly on the experimental time scale.

Fig. 5 illustrates how the available detachment processes can be interpreted in terms of different valence bond structures and how the ability to transition between electronic states adiabatically can be rationalized using a simple orbital picture. Analysis of [C] 121 ²A″ (2et), [D] 121 ²A″ (1et2), and [F] 121 ⁴A″ (2et) structures involves Mo(IV) and V(IV) oxidation states, while [B] ²A″ (1et) and [E] ⁴A″ (1et) involve Mo(V) and V(III) oxidation states. Thus the metals are either d¹ or d² and are able to couple magnetically through bridging oxygen p-orbitals. In the case of the 121 ⁴A″ states, the electronic structure of the [F] ⁴A″ (2et) structure (Fig. 5(f)) is obtained from the [E] ⁴A″ (1et) structure (Fig. 5(e)) by transfer of an electron from vanadium to molybdenum. As all electrons in the metal d-orbitals have the same spin, spin-conserving electron transfer likely occurs without significant energy penalty, coupled through nuclear vibrations of the bridging oxygen centers. The lateral position of the bridging oxygens indicates the position of the transferred electron. The low energy barrier and the lower energy of the minimum likely result in the [F] ⁴A″ (2et) structure (Fig. 5(f)) to be favored in this case. After photodetachment, it is expected, therefore, that the [F] ⁴A″ (2et) structure corresponding to the two-electron transition would be present, though, spectroscopic measurements would be significantly less intense than for one-electron transitions.

It is possible that another triplet ground state exists with triplet-coupling of Mo electrons forming an (overall) triplet with the two V electrons being singlet-coupled. This would allow direct access to the [F] 121 ⁴A″ (2et) state upon photodetachment. However, we were unable to locate such an electronic structure, and although it is not possible to conclusively state that this alternative triplet state is not present, within the constraints of the tractable and available theory at our disposal this appears to be a reasonable conclusion. This conclusion is further supported by the experimental photoelectron spectrum which would not display the weaker X band intensity if the transition was one-electron accessible.

As shown in Fig. 5, the antiferromagnetically coupled metal centers facilitate a superexchange mechanism in the ²A″ structures. Similar to the situation observed in the quartet, [B] 121 ²A″ (1et) (Fig. 5(b)) and [D] 121 ²A″ (1et2) (Fig. 5(d)) structures lie on the same PES connected by a spin-conserving electron transfer. [C] 121 ²A″ (2et) (Fig. 5(c)) is the lowest energy structure; however, there is no direct, one-electron pathway to concurrently detach an electron and rearrange the electronic structure. Instead, such a change requires the exchange of electron density between metal centers from the lowest VDE state [B] ²A″ (1et). Owing to the extensive electron rearrangement occurring in response to electron detachment, we expect the [C] ²A″ (2et) state is not readily accessed. The [B] 121 ²A″ (1et) structure (Fig. 5(b)) corresponding to a one-electron transition is therefore the relevant minimum even though the [C] 121 ²A″ (2et) (Fig. 5(c)) structure is lower in energy.

To summarize, analysis of the different accessible spin channels suggests that transition from the anion ³A′ to the neutral ²A″ state is a one-electron transition, while transition to the neutral ⁴A″ state arises from a two-electron transition. The different selection rules for one-electron and two-electron transitions provide an explanation of the difference in intensity between the experimental X and A bands.

The transitions responsible for the observed FC profile of $MoVO4−$ were assigned as the 121 ³A′ → 121 ⁴A″ + e⁻ transition for the X band and the 121 ³A′ → 121 ²A″ + e⁻ transition for the A band. A potential third transition, resulting from 220 ⁵B₁ → 220 ⁴B₁ + e⁻, was excluded due to experimental evidence suggesting that the X and A bands arise from a common isomer and discrepancies between computed VDE and ADE and observed experimental values. Although the 220 ⁵B₁ structure was computed to be lower in energy that the 121 ³A′ structure, further justification for the exclusion of the 220 ⁵B₁ → 220 ⁴B₁ + e⁻ transition is based on the tendency of hybrid functionals to overstabilize higher spin states.

Discrepancies between simulated FC profiles in previous reports and the experimental photoelectron spectrum have been shown to result from significant geometric distortion along either very low frequency vibrational modes or over low energy barriers resulting from Jahn-Teller distortions. The presence of flat or multiple minima on the PES causes FC simulations, which assume a single well-defined minimum in both the ground and excited states, to be invalid and yield poor agreement with experimental and simulated FC profiles. In order to correctly simulate FC profiles, Boltzmann weighted structural averaging was used as well as exclusion of imaginary and low frequency vibrational modes from the FC simulation.

A further consideration is the presence of a number of broken symmetry solutions in the neutral states that complicate the assignment of the observed transitions. Broken symmetry density functional theory calculations were able to provide both the energy and electron distribution of these states. Within a single-reference approach such as DFT, transition structures may qualitatively indicate regions where different states approach each other and adiabatic transitions between states can occur. In $MoVO4−$⁠, the pertinent transition vector involves the lateral shifting of the bridging oxygens and is associated with electron transfer between the metal centers. Rationalizating electron transfer mechanisms using an orbital analysis, it was possible to identify the 121 ³A′ → 121 ⁴A″ + e⁻ transition as a two-electron process while the 121 ³A′ → 121 ²A″ + e⁻ transition was determined to be a one-electron event. Accounting for differences in transition types gave agreement with the observed experimental intensity difference between X and A bands.

See supplementary material for geometries, energies, and spin densities of computed anion and neutral $MoVO4$ structures as well as for the geometry used to compute natural ionization orbitals on the 121 ³A′ → 121⁴A″ (2et) transition.

We gratefully acknowledge financial support from the University of California, Merced, and computing time on the Multi-Environment Computer for Exploration and Discovery (MERCED) cluster which is supported by National Science Foundation Grant No. ACI-1429783.