Optimal photoelectron circular dichroism of a model chiral system Open Access

Photoelectron circular dichroism (PECD^1–3) is a fascinating chiroptical effect of a forward–backward asymmetry in the photoemission from randomly oriented chiral molecules with respect to the propagation direction of the ionizing circularly polarized light. Predicted by Ritchie⁴ about 50 years ago theoretically, it was first verified in one-photon^5,6 and then in multiphoton^7,8 ionization experiments with chiral molecules. At the present time, it is considered to be universal with respect to the photoionization regime.^9,10 Because PECD is purely an electric-dipole effect, its strength typically reaches a few percent, which is much larger than the conventional circular dichroism (CD¹¹). As a consequence, PECD can efficiently be applied to even a sub-percent enantiomeric excess determination¹² in the gas phase.^1–3 Since larger PECD strengths enable higher sensitivity in the enantiomeric excess determination, several attempts have already been made to maximize this chiral asymmetry. One of the discussed routes is to optimize the properties of the ionizing light using coherent control schemes,^13,14 which may yield a PECD well above 30%. Another very different route proposes fixing molecular orientation partly^15,16 or even fully^17,18 in space. This strategy of optimizing the properties of targets as a whole may enhance the chiral asymmetry well beyond 50%.

At present, it is understood that the magnitude of the PECD is governed by the ability of the outgoing photoelectron to probe a chiral asymmetry of the molecular potential through multiple scattering effects. Therefore, a potential route to enhance the PECD itself would be an optimization of the intrinsic properties of chiral targets to generate the strongest asymmetry of its chiral potential. For this purpose, one needs to manipulate a charge density distribution in the ion, which generates an electrostatic potential for the outgoing photoelectron. This can be achieved by changing its nuclear charges, its internuclear geometry, or the distributions of the bound electrons when, e.g., constituent atoms get exchanged or nuclear charges get merged. For instance, one may think of transforming the three terminal CH₃ groups in the small chiral molecule fenchone independently into NH₂, OH, or F, which already yields 4³ different iso-electronic chiral modifications of this molecule from three sites alone. The fundamental question immediately arises: Which of the modifications has the strongest PECD? Unfortunately, the present knowledge of the PECD phenomenon does not allow us to answer this question without enumerating all compounds computationally.

Before performing optimization of large and real chiral molecules with plenty of degrees of freedom, it is extremely important to understand the individual impacts of different inherent properties (nuclear charges, bond lengths, electron densities, etc.) on this chiroptical effect in a systematic way. In the present work, we make the first step toward this aim. In particular, we use a minimal model chiral system with a few tuneable properties and perform theoretical optimization of the chiral asymmetry of its electronic potential to achieve the largest PECD effect. Even for this relatively small theoretical problem, performing ab initio electronic structure calculations of the PECD in the whole space of free parameters is not feasible. Therefore, we employ quantum alchemy, which relies on a Taylor series expansion^19–21 and can estimate the properties of iso-electronic molecules for millions of molecules at once without enumeration. The paper is organized as follows: The present model and theoretical methods are described in Sec. II. The results of the optimization are reported and discussed in Sec. III. We conclude in Sec. IV with a brief summary and outlook.

In this section, we first discuss the presently used model chiral system (Sec. II A) and then outline theoretical methods for the electronic structure and PECD calculations (Sec. II B) and optimization (Sec. II C).

The presently used model chiral system is very similar to those introduced in our previous works.^22,23 It is built of model atoms, each represented by a positive point charge Q, which simulates a net nuclear charge being partly compensated by the inner electrons bound to this atom. In addition, we introduce an electronic cloud with a negative charge q and a spherically symmetric distribution e^−r, which is supposed to represent the outer electrons shared in the molecule. In order to be chiral, there should be at least four non-coplanar and non-equivalent atoms. Even for four atoms, there are nine independent degrees for optimization of the internuclear geometry. In the present model study, we restrict the number of geometrical optimization parameters to three. For this purpose, we place one of the atoms at the coordinate origin (Q₀, q₀) and three on the respective coordinate axes (Q_i, q_i), and assume the distances to the origin along the axis (R_i) to be free parameters. Thereby, for given distances R_i, the four atoms are always non-coplanar in a maximal way. The initial values for the geometry parameters are R_x = R_y = R_z = 3.0 a.u., which corresponds to a typical molecular bond length. Figure 1 depicts a sketch of the present chiral system and discusses quantities.

In addition to geometry properties, we optimize the charge density of the model system as well. For simplicity, we fix all charges in the electronic clouds and vary only the nuclear charges. In particular, the electronic charges are fixed in a very asymmetric way to be q_x = −1, q_y = −2, q_z = −3, and q₀ = −4 (see Fig. 1). Since the model molecular ion contains ten shared electrons, ∑_iq_i = −10, the total uncompensated nuclear charge must be set to ∑_iQ_i = +11, thereby modeling a singly charged ion with an outer electron being ionized. We assume the three nuclear charges Q_x, Q_y, and Q_z to be free parameters and adjust the central charge Q₀ accordingly. The initial values of the nuclear charges are Q_x = Q_y = Q_z = +2.5, and Q₀ = +3.5. The one-particle potential, generated by this model ion for the photoelectron, is used to solve the stationary Schrödinger equation. The first five eigenstates of the Hamiltonian are assumed to be populated by ten common bound electrons, while the sixth eigenstate represents a singly occupied HOMO of this open-shell model system, which is ionized by circularly polarized light. This eigenstate is used as the initial state of the photoelectron for the time-dependent Schrödinger equation (see Sec. II B).

In the present study, we investigate the non-resonant energy-dependence of the anisotropy and dichroic parameters. We, thus, neglect autoionizing electronic resonances in the calculations. The respective transition amplitudes were computed by the time-dependent realization of the single center (TDSC) method,²² which was successfully applied to study PECD in the one- and multiphoton ionization of chiral molecules.^27–30 The working equations of the method and details on the calculations can be found in those works. In the present case, we employed single center expansion with ℓ, |m| ≤ 15, used weak exponential laser pulses $g(t)=e−(t−t0)2/τ2$ with t₀ = 3 fs and τ = 1 fs, and propagated the initial electronic wave packet for the linear and both circular polarizations during 6 fs in a radial box of 250 a.u. The pulse intensity of 10¹⁰ W/cm² was chosen to suppress the above threshold ionization by two-photon absorption. For the initial parameters of the model chiral system discussed in Sec. II A, the binding energy of its sixth electronic eigenstate is equal to 19.94 eV. The photon energies were adjusted accordingly to release photoelectrons with kinetic energies of 4, 5, 6, 7, and 8 eV. The anisotropy and dichroic parameters, computed for this initial reference set (the first line in Table I), are depicted in Fig. 2 by black stars. As one can see in the chosen energy interval, this system exhibits a sizable PECD characterized by a dichroic parameter between about +2% and +6%.

Many libraries obtain finite difference derivatives from a fixed stencil, i.e., the set of all function evaluations entering a numerical approximation. We implemented the multidimensional multivariate arbitrary order Taylor expansion from arbitrary evenly spaced finite difference stencils in the Python package nablachem.MultiTaylor³¹ using findiff.³² This library also allows us to select the mixed partial derivatives to be included selectively and performs L-BFGS-B optimization over the resulting multivariate polynomial. We always choose the maximum number of available points for stencils of lower than maximal order; e.g., first order derivatives are evaluated using all available and applicable data points rather than only two. This improves the accuracy of lower derivative orders. All monomials in the Taylor expansion for which there are not sufficiently many coefficients available are assumed to be zero. In this work, this concerns terms higher than the respective order or of mixed partial derivatives including more than two different variables. Starting from the center of the series expansion, the bounded limited memory optimizer L-BFGS-B with a budget of 1000 steps was used to find extremal points within constraints on each dimension in the six- or seven-dimensional search space of β₁(x), while β₂ was just evaluated at the resulting point. The bounds represent the unknown trust radius of the truncated Taylor expansion, since at large changes in charges Q_i or atom position R_i, instability due to the highest order polynomial term is unavoidable. While bounds for the values of β_i are known (|β₁| ≤ 1 and −1 ≤ β₂ ≤ 2), their limits for R_i, Q_i → 0 or ∞ are unknown. Because of that, Padé approximants, which often exhibit advantageous convergence properties,²⁰ could not be used.

In order to initiate the screening procedure, ab initio calculations were performed first for the following sets of parameters: Starting from the initial reference set Q_x = Q_y = Q_z = +2.5 and R_x = R_y = R_z = 3.0, each of the six parameters was varied independently in steps of ±0.05 and ±0.10, keeping the other five parameters unchanged. The respectively modified system was used to construct an ionic potential for the photoelectron, and, subsequently, anisotropy and dichroic parameters β₁ and β₂ were computed by the TDSC method for the photoelectron kinetic energies 4–8 eV. This yields 4 × 6 = 24 calculations for each of the five considered energies. In addition, in order to be able to evaluate some multivariate partial derivatives, two (out of six) parameters were simultaneously varied in steps of ±0.05 for all combinations. This yields $4×6⋅52=60$ additional calculations for each energy. The results of these calculations for 85 sets of the variational parameters (including the initial reference set) are provided in the supplementary material.

These ab initio data were used to calculate the derivatives for the Taylor series in expansion (2), which was employed to screen the search space for free parameters at which the value of dichroic parameter β₁ reaches its desired value. For reliability, the second order Taylor expansion was used for those predictions, even though terms of the third and fourth order expansions were accessible if they did not include partial derivatives in more than two or one variable, respectively. Even using the full information provided by the calculations, for certain relatively large variations of the bond lengths (R_i) and charges (Q_i), the Taylor expansions become inaccurate, and thereby the extrapolated values of the dichroic and anisotropy parameters exceed their validity intervals. Nevertheless, before diverging for larger parameter variations, the algorithm discovered (for each considered energy between 4 and 8 eV) a local maximum of the dichroic parameter β₁ that fell in the validity interval |β₁| ≤ 1, and the corresponding anisotropy parameter fulfilled the condition −1 ≤ β₂ ≤ 2. Those predicted values of β₁ and β₂ parameters were additionally checked to yield a non-negative emission probability (1) at all emission angles θ. Moreover, considering the photoelectron energy as an additional degree of freedom, we were able to extrapolate the predictions outside the chosen energy interval of 4–8 eV. The explored values of β₁ and β₂ were verified by the ab initio calculations using predicted values of the variational parameters (R_i) and (Q_i). The respective results are discussed in Sec. III A, while a detailed analysis of the present screening procedure is given in Secs. III B and III C.

We start the optimization in the chosen photoelectron energy interval of 4–8 eV and constrain our search space to maximum changes in ΔR = ΔQ = 0.1. Here, each energy is considered independently, i.e., no partial derivatives between the photoelectron kinetic energy and other parameters have been included, which allows us to evaluate the smoothness of the predictions.

The values of the dichroic and anisotropy parameters predicted for the optimization result set 1 from Table I are depicted in Fig. 2 by red dashed curves. The top panel of this figure predicts a maximal value of $β1e(4eV)=24.1%$⁠. However, the ab initio calculations yield a somewhat smaller value of the dichroic parameter for that configuration, and this energy $β1c(4eV)=16.1%$ (red triangles). As one can see, the predicted values of β₁ are overestimated and those for β₂ are underestimated (cf. the red dashed curve and red triangles in Fig. 2), while the overall trend is well-reproduced, as exemplified by the correct identification of a turning point in β₂ at about 5 eV. Following the computed trend toward smaller kinetic energies, we find a maximal value of about $β1c(1.94eV)=25.3%$ at the kinetic energy of 1.94 eV. Here and below, we list spline-interpolated extreme values of the computed dichroic parameter at non-integer energy values. In comparison, the largest value for β₁ for all ab initio data points in the finite difference stencil is 9%. Since there is no formal convergence guarantee for the expansion of many properties, including β_i, this numerical evidence is important as it indicates that both the dichroic and anisotropy parameters can be optimized via perturbative treatment and are sufficiently smooth and differentiable to allow for this gradient-based optimization. Moreover, we find that the predictions in both parameters β_i yield clear trends and smooth curves along the photoelectron energy, even though each energy value of 4–8 eV has been predicted independently. This makes numerical effects an unlikely cause for the overestimation of the magnitude of β_i but rather points toward relevant higher order terms.

In the second step, we now consider photoelectron energy as a free parameter and search for optimal β₁ values within the extended interval of 2–10 eV. To this end, we have to include mixed partial derivatives between the photoelectron energy and all other parameters, i.e., we only build one model for the whole energetic, spatial, and electronic domain covered by the stencil.

From the results of set 1, one might expect the largest prediction value for energies below 4 eV. In addition, the optimization on the Taylor expansion yields a new point, set 2 (see Table I), which is predicted to have the largest dichroic parameter value at the boundary of the search interval, at 2 eV. For this system, the perturbation expects an almost linear dependency of the dichroic parameter on the photoelectron energy; see the blue dotted curves in Fig. 2(a). For the domain of 4–8 eV, we can compare this to the individual results, where five independent models predict the same configuration in set 2 without the inclusion of any mixed partial derivatives with the energy. Doing so would yield the same linear trend, although with a steeper slope for β₁ and even the wrong sign of the curvature for β₂. Since either method overestimates the parameter magnitude, this points toward those mixed partial derivatives containing valuable information about the interaction between the energy parameter and the spatial and electronic degrees of freedom. This will be important to consider when optimizing β_i for real molecules.

Extending the configuration of set 2 (see Table I) toward even lower photoelectron energies does not reproduce the local maximum in β₁ from the ab initio calculations around 2 eV, so second order mixed derivatives in the photoelectron energies are not sufficient. Note that while fourth-order derivatives in energy alone are available (since the data set contains the reference system for five different energy values), this would only apply to a frozen spatial and electronic configuration. Once those dimensions are allowed to change, their interplay can only be captured up to second derivatives in most cases, given our particular stencil. The optimal predicted value of the dichroic parameter is $β1e(2eV)=32.7%$ [see Fig. 2(a) and Table I], while the ab initio calculations yield a somewhat smaller value of $β1c(2eV)=25.8%$⁠, and the computed dichroic parameter reaches its maximal value of $β1c(1.56eV)=26.3%$ at the lower energy of 1.56 eV. As shown in Fig. 2(b), predictions of β₂ are of comparable accuracy for both set 1 and set 2. Given the perturbative nature of the model, some error from higher order terms is expected, but trends are well-preserved in the optimization procedure: just as predicted, the value of the dichroic parameter computed for the parameter set 1 at 4 eV is indeed slightly larger than that computed for set 2, and, simultaneously, its value computed for the parameter set 2 at 2 eV is slightly larger than that computed for set 1 [see Fig. 2(a)].

Finally, we search for configurations where the forward–backward asymmetry in the photoemission changes its sign, i.e., where β₁ is negative. The optimizer suggests set 3 (see Table I) as the point of minimal β₁, where the predicted dichroic parameter reaches $β1e(2eV)=−5.1%$⁠. Ab initio calculations confirm this negative value and yield even stronger negative asymmetry of $β1c(1.7eV)=−14.7%$⁠. While this is in line with the optimization goal, the prediction accuracy for both parameters in this task is worse than in the previous tasks. Potentially, a smaller displacement in the finite difference stencil may improve the derivative calculation, as the noise of the ab initio calculations upon small perturbations during the validation of the prediction was found to be negligible.

As discussed in Sec. III A, the models built from 85 single point calculations can be used to optimize the PECD parameters β_i. This is computationally feasible now because the models are much faster to evaluate than one ab initio calculation. We can now use these models to draw statistical conclusions about which terms typically contribute strongly to the values of β_i: All terms that typically do not contribute significantly are good candidates to neglect to reduce the number of high-quality ab initio calculations, as this allows us to reduce the stencil if some derivatives are not needed. Going to larger molecules, where the number of derivative terms increases due to the total number of sites, this is crucial to maintain computational feasibility.

Figure 3 shows the impact of the individual derivatives grouped by their order and combination of electronic, spatial, and charge dimensions as obtained from one million model estimations. This allows us to compare which interaction terms are relevant for the PECD parameters β_i. We first notice that higher orders typically contribute substantially less: First order terms typically contribute 10%–35% to the prediction value, second order terms 5%–10%, and higher orders typically contribute less than 3%. This is reassuring, as this motivates the truncation of the series expansion after a few orders. For design applications, minor deviations from the true value are acceptable since the optimization can be performed using the model and then validated using a few accurate ab initio calculations.

For both parameters, the electronic derivatives, i.e., the dependency on the photoelectron energy to first order, are by far the largest contribution. This means that future applications to molecules will surely require data from multiple photoelectron energy values. It is interesting to note that this term loses importance with higher orders quite quickly. Since there is only one input dimension related to the photoelectron energy, all even order contributions can either be positive or negative, as is evident from Fig. 3. Even though the distributions have a long tail, this tail is particularly thin. This likely means that the inclusion of third order energy derivatives is sufficient.

For the charge parameters, two orders cover almost all effects. It is interesting to see that for β₂, all second order derivatives of the charges are positive. This means that the model geometry is a local minimum in β₂. For β₁, a highly asymmetric distribution in the second order contributions indicates the same effect for part of the charge dimensions. Higher order terms have negligible contributions. Therefore, second order derivatives are likely sufficient for molecular design.

For the spatial contributions, only first order terms are significant, which is particularly helpful for molecular design since this means that the response of β_i due to the changed geometry of the design target is comparably easy to capture with a linear response. This is fortunate because there are three times as many spatial dimensions as there are charge dimensions to consider. Since obtaining the approximate minimum energy geometry for molecules is much cheaper than obtaining new PECD data via ab initio calculations, this offers a clear path toward a quick assessment of the molecular PECD parameters.

Since the mixed terms that couple different kinds of degrees of freedom may exhibit correlations, we analyzed their impact from ablation, i.e., dropping all groups of derivative contributions from the full set until the mean absolute error for a quantity reached 20% of the average quantity value. We expect this to be sufficiently accurate for iterative optimization and in line with the accuracy estimation from the first part of the results. We find that the second order spatial and charge degrees of freedom mix with the photoelectron energy but not among each other for both parameters. In the third order, only linear response, spatial, and charge parameters are coupled with the photoelectron energy. For β₂, the linear response of the nuclear charges couples strongly with the second order photoelectron energy response.

Because each derivative coefficient requires one more ab initio calculation, it is desirable to define the minimal set of derivatives that is needed to obtain a quantum alchemy prediction with predictive power for β₁. This is amplified in the case of treating real molecules with N sites to change, since there are 3N spatial and N charge degrees of freedom to consider. Each of the derivative terms has a different scaling with N: There is always a constant number of derivatives of order k with respect to the photoelectron energy, but, e.g., (3N)!/(3N − k)! many with respect to spatial degrees of freedom. Assuming a molecule to have ten atoms that could be altered, Fig. 4 shows the ablation results if individual terms are neglected in relation to the cost it would take to obtain these terms. There is a striking general trend that more expensive terms also contribute less—a strong argument for neglecting them to improve computational efficiency. Leaving out all terms in order of their computational cost until a mean error in the prediction of 20% is reached leaves 11 (mixed) derivatives, which for the system in this work require 29 ab initio calculations rather than the 85 done initially. For larger non-linear molecules with N atoms to change, there are 3N coefficients for all R1 terms: N for Q1, 3N for R1E1, N² for Q2, N for Q1E1, 3N for R1E2, N for Q1E2, and 1 each for E0 to E4. Subtracting six coefficients to account for rotation and translation yields N² + 12N − 1 calculations to determine all coefficients. Note that the results for the system in this work deviate from this formula since each of its atoms has only one spatial degree of freedom.

This formula allows us to estimate the molecular size for which this perturbative approach is more efficient compared to brute force enumeration of the target search space. For energies, density matrix elements, orbital eigenvalues, and one-electron properties in real molecules, quantum alchemy typically²⁰ converges for at least ΔZ = ±1. Here, this would allow three elements per atom, which can be varied during the search. For N changing sites, this allows for 3^N targets to be estimated from the same reference perturbations. This search space clearly grows faster than the quadratic number of derivatives that need to be determined. With four atoms, the quantum alchemy approach becomes faster, and with nine atoms, it is already 100 times faster. It is worth noting that for larger systems, the overall electronic structure’s response to individual local changes becomes less pronounced. This would likely increase the convergence of the perturbation, which is in favor of the method.

In the present work, bond lengths and electron density distribution parameters of the model chiral system were optimized to achieve the desired values of the PECD. The screening of the variational space was performed by the Taylor series expansion algorithm using a relatively small set of ab initio computed data around the reference configuration. We were able to increase the dichroic parameter β₁ by 3 to 4 times as compared to the reference geometry and, for some photoelectron energies, to achieve its value of around 26%. Keeping in mind that the dichroic parameter emerges on a typical scale of a few percent,^1–3 the obtained huge value is rather unusual for randomly oriented targets, where sizable contributions from different molecular orientations (on the order of 20%–50%^15–18) compensate each other almost completely, resulting in a much smaller effect. In addition, we were able to predict the parameter set at which the PECD changes its sign.

The results allowed us to dissect the contributions of different physical interactions of spatial and electronic degrees of freedom with the charge distributions, and that analysis offers clear information about which derivatives are likely to be of importance for molecular design to optimize the PECD parameters. We found that all expensive high-order contributions have only little impact on the overall prediction, and, therefore, the perturbative approach offers a computationally advantageous setting for systems with many atoms in the design space. This is particularly interesting as—unlike other electronic properties of molecules studied thus far—all degrees of freedom couple strongly for PECD parameters. This means that molecular design likely requires an iterative procedure where changes to the composition and geometry are suggested repeatedly in order to obtain large changes. The current work, however, showed that the one-step optimization already achieved a tremendous increase in the dichroic parameter.

In addition to the investigation of molecules, future work might explore whether theoretical limit expressions can be derived for β_i, which could then be used in Padé approximants rather than Taylor series, or whether variable transformations allow for the disentanglement of the interplay of all degrees of freedom. Both routes would improve predictive accuracy, allowing further reduction of the computational cost. Automatic differentiation as applied to other ab initio calculations^33–36 may be a promising alternative to finite difference evaluation, further reducing computational costs for design applications.

The supplementary material for this article includes the full set of the ab initio β1 and β2 parameters, computed at different initial sets of the variational parameters and used for the present optimization procedure.

A.N.A. and P.V.D. acknowledge the support from the Deutsche Forschungsgemeinschaft (DFG)—Project No. 328961117—SFB 1319 ELCH (Extreme light for sensing and driving molecular chirality, subproject C1).

The authors have no conflicts to disclose.

Guido F. von Rudorff: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Anton N. Artemyev: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Boris M. Lagutin: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Philipp V. Demekhin: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

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