Phase Space Electronic Structure Theory: From Diatomic Lambda-Doubling to Macroscopic Einstein–de Haas

TL;DR

This paper shows how a phase space electronic structure theory can explain both microscopic Λ-doubling and macroscopic angular momentum effects like the Einstein–de Haas phenomenon.

Contribution

A phase space method is shown to nonperturbatively capture Λ-splitting in NO without summing over states.

Findings

The phase space method recovers the Λ-doubling energy splitting of the NO molecule nearly quantitatively.

The method correctly conserves angular momentum and includes electron-rotation coupling in potential energy surfaces.

The computational cost is comparable to standard Born–Oppenheimer calculations.

Abstract

Λ-doubling of diatomic molecules is a subtle microscopic
phenomenon that has long attracted the attention of experimental groups,
insofar as rotation of molecular nuclei induces small
energetic changes in the (degenerate) electronic state.
A direct description of such a phenomenon clearly requires going beyond
the Born–Oppenheimer approximation. Here we show that a phase
space theory previously developed to capture electronic momentum and
model vibrational circular dichroismand which we have postulated
should also describe the Einstein–de Haas effect, a macroscopic
manifestation of angular momentum conservationis also able
to recover the Λ-doubling energy splitting (or Λ-splitting)
of the NO molecule nearly quantitatively and nonperturbatively (without
a sum over states). The key observation is that, by parametrizing
the electronic Hamiltonian in terms of both nuclear position (X) and nuclear momentum (P), a phase space method
yields potential energy surfaces that explicitly include the electron-rotation
coupling and correctly conserve angular momentum (which we show is
essential to capture Λ-doubling). The data presented in this
manuscript offer another small glimpse into the rich physics that
one can learn from investigating phase space potential energy surfaces E

PS
(X,P) as a function
of both nuclear position and momentum, all at a computational cost
comparable to standard Born–Oppenheimer electronic structure
calculations.

Figures

Potential energy surfaces of the ground and first excited states of NO from the phase space Hamiltonian, plotted versus the canonical nuclear angular momentum L. The blue and orange curves correspond to states whose valence Π* molecular orbitals are approximately 12[|Px⟩|↓x⟩−i|Py⟩|↑x⟩] and 12[−|Px⟩|↑x⟩+i|Py⟩|↓x⟩] , respectively. The nuclear coupling breaks the Kramers’ degeneracy of the two PESs at a given finite canonical nuclear rotation angular momentum L and form a double well with minima at L/ℏ ≈ ±0.0022.

Energy diagram of the low-energy space of NO. Energy splittings are not drawn to scale.

Λ-Splitting between the lowest e and f states and the associated α parameter for a nuclear rotation J = 1/2ℏ around one perpendicular axis. The results are calculated using the exact FCI solver (blue circle) in the STO-3G basis, the DMRG solver (purple star) which converges to the FCI accuracy in the STO-3G and 6-31G bases, and an approximate constrained Hartree–Fock solver in larger bases, from STO-3G, 6-31G, cc-pVDZ, aug-ccpVDZ, to aug-ccpVTZ (red).

Energy splitting between the e and f states as a function of total angular momentum J from the 2D model parametrized with α derived from FCI in the STO-3G basis (blue) and from cHF in the aug-ccpVTZ basis (red), both using Method III, in comparison to experiments (black). Dashed lines include only the 1-electron SOC, while solid lines include both 1-electron and 2-electron SOC (with a spin–orbit mean-field approximation , ).

a) Our original view (see Figure ) of potential energy surfaces of the ground and first excited states of NO according to a phase space electronic structure Hamiltonian as a function of the phase space canonical nuclear angular momentum L. b) The same data now labeled differently: The green and red curves can be understood as the e and f parity total electronic and nuclear states, now understood as functions of total angular momentum J (represented by the phase space L). The e/f states are the states that are observed spectroscopically.

Energy splitting between the e and f states of the 2Π3/2 manifold as a function of total angular momentum J. For this data, unlike Figure , it is well-known that the energy splitting requires a cubic (rather than linear) fit as a function of J. We compare phase space predictions, using a model parameter α derived from FCI with 1-electron SOC in the STO-3G basis with Method III (blue dots) versus experiment (black triangles). We further fit the “phase space” data to an analytic expression from standard perturbation theory, ΔE=k(J+12)(J+32)(J−12) , where k is a fitted constant (dashed lines). Note that, as in Figure above, our phase space predictions are quite accurate, which is another strong confirmation of a PS approach.

Tables

1. Coupling Parameter α Estimated in Three Ways: Method I: from ⟨L̂ a e + Ŝ a e ⟩ Evaluated on the Parity Eigenstates at L = 0; Method II: from L a at the Minimum of the PES; Method III: from Energy Splitting ΔE at L a = 1/2ℏ

Method	I	II	III
α	0.0018	0.0022	0.0022

Full text

Introduction

1

Conservation of angular momentum is a fundamental principle that governs both everyday life and microscopic dynamics. As far as everyday life is concerned, angular momentum conservation leads to the Coriolis force that drives the clockwise subtropical Pacific Ocean gyre? in the Northern Hemisphere while also allowing a rolling hoop to stay upright.? Because of angular momentum conservation, whenever one part of a system changes its state of motion, the remaining degrees of freedom must respond so that the total momentum and angular momentum are preserved?a concept we are all well acquainted with. On a microscopic level, the same logic applies to electrons and nuclei in molecules and solids: when the nuclei move, the electrons must adjust their states to the movement of the nuclei. Just like Foucault pendulums, electrons can experience Coriolis forces. ?−? ? ? ?

Unfortunately, most modern-day electronic structure theories rely on the Born–Oppenheimer approximation and ignore the direct impact of nuclear momentum on electrons.? To go beyond this limitation, we recently introduced a phase space theory that explicitly parametrizes the electronic Hamiltonian by both nuclear position and momentum and conserves the total momentum and angular momentum of electrons and nuclei. We have already demonstrated that the method can recover vibrational circular dichroism ?,? and Raman optical activity;? for model systems, the approach can also offer improved vibrational energies. ?,? Moreover, because angular momentum conservation is guaranteed, we have postulated that the approach should be useful for modeling chiral-induced spin selectivity? as well as the Einstein–de Haas (EdH) effect. ?,? The latter phenomenon remains of great interest to physicists and chemists alike; ?,? as a brief review, EdH is the phenomenon whereby changing the magnetic moment can cause a ferromagnet to rotate, and thus serves as a macroscopically large example of angular momentum conservation in solids.

In the realm of small molecule chemistry, here we will now show that the same beyond BO physical concepts described above also naturally describe Λ-doubling. In particular, using the phase space electronic structure theory from refs ?−? ? , we will show that one naturally recovers the Λ-doubling splittings of the NO molecule quantitatively. Thus, we will argue that, whether over the smallest or largest length scales, electronic structure theory can be improved by moving to a phase space framework where the Schrodinger equation is solved in the frame of moving nuclei.

An outline of this article is as follows. First, for the sake of pedagogy, we will briefly review phase space electronic structure theory. Second, we will then discuss the physical meaning of phase space potential energy surfaces, showing that the energy surface minimum captures the strength of electron-rotation coupling for rotations around a given molecular axis. Third and finally, combined with a model of two-dimensional rotations, we will show how a phase space approach predicts the Λ-splitting of NO as observed in experiment.

As far as notation is concerned, as written below, boldface implies a 3- or 3N-dimensional vector. Operators are given hats. We will write the electronic angular momentum equivalently as L _ e _ or L ^ e ^ (depending on whether we have extra superscripts, e.g., L _ e _ ^+^, or extra subscripts, e.g., L _ x _ ^ e ^). Within phase space electronic structure theory, we label the canonical nuclear linear momentum P and the angular momentum L; however, as will be shown below, the canonical angular momentum within a semiclassical phase space theory represents the total angular momentum, so that for a system with electronic orbital and spin degrees of freedom, the L within a phase space calculation can be thought of as the total J within a fully quantum calculation. The notation should be clear from context.

Brief Review of the Phase Space Method

2

Before we address Λ-doubling for NO, let us briefly review phase space electronic structure theory. Phase space electronic structure theory provides a semiclassical means to go beyond the standard BO approximation by parametrizing the electronic Hamiltonian using both the classical nuclear position (X) and nuclear momentum (P). More than ten years ago, Shenvi wrote down the first such example,? suggesting that we diagonalize:

$$[eqn]$$

Here, X and P are classical c-numbers, Ĥ el|Φ_ I _⟩ = E _ I |Φ I ⟩ is the diagonalization of the standard electronic Hamiltonian, $[eqn]$ , where A is the atom index, is the derivative coupling vector between the many-body eigenstates |Φ I ⟩ and |Φ J _⟩, and M _ A _ is the mass of the atom A.

Unfortunately, eq is of limited value for several reasons: (i) the derivative coupling can be numerically unstable near curve crossings, (ii) the derivative coupling is ambiguously defined for degenerate BO eigenstates, and (iii) the procedure above requires two diagonalizations just for the electronic part and is computationally too demanding in practice. That being said, over the last several years, we have shown that progress can be achieved by approximating the derivative coupling vector d _ IJ _ ^ A ^ with the matrix element ⟨Φ_ I |Γ̂ A |Φ J ⟩ for a one-body operator Γ̂ A _. To ensure that our approximation is as robust as possible, we have posited that one must preserve the analogs of four important symmetry constraints that derivative coupling vectors satisfy, namely,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Equations and ? reflect translational invariance and rotational invariance of the total system, so that when nuclei are displaced translationally, the electrons are displaced with them, and when nuclei rotate, the electrons rotate with them. Equations and ? enforce that Γ̂_ A _ itself is invariant to translation and rotation.

The form of the one-body operator Γ̂_ A _ we will use that satisfies all constraints above contains three terms, Γ̂_ A _ = Γ̂_ A _ ^′^ + Γ̂_ A _ ^″^ + Γ̂_ A _ ^‴^. To write down these operators, let us define a partition of unity of space according to nuclei as

$$[eqn]$$

where σ is a locality parameter in the unit of length. The three Γ operators are then: (i) The electron translation factor Γ̂′_ A _

$$[eqn]$$

­(ii) the electron orbital rotation factor (ERF) Γ̂″_ A _

$$[eqn]$$ $$[eqn]$$

and (iii) the electron spin rotation factor Γ̂_ A _ ^‴^

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a local average position,

$$[eqn]$$

is a local version of moment of inertia partitioned to each atom, and $[eqn]$ is the 3 × 3 identity matrix. For the case of a linear molecule (as below) where K _ B _ is not invertible, we replace K _ B _ ^–1^ Ĵ _ B _ with?

$$[eqn]$$

where n 3 is the unit vector along the internuclear axis of the linear molecule. Above, ζ_ AB _ = M _ A _e^–|X _ A _ – X _ B _|^2^/8σ^2^ ^ reflects the degree of electronic/spin communication between different atoms and must be a function that decays on a length scale commensurate with the decay of the electronic density matrix.

Within our phase space electronic Hamiltonian, the Γ̂_ A _ operator above creates an approximate Hamiltonian mimicking the derivative coupling:

$$[eqn]$$

As with standard electronic structure, eq can be solved self-consistently using a standard electronic structure solver with almost the same cost as the BO Hamiltonian. But in sharp contrast to the BO Hamiltonian, we now have an explicit electron–nuclear momentum coupling term, $[eqn]$ that provides the potential to capture various beyond BO effects.

For example, when the symmetry constraints above are satisfied, the coupling term exactly captures the noninertial Coriolis effect on electrons in the rotating frame of rigid nuclei. To verify this statement, imagine that we include Γ̂_ A _ ^′^ and Γ̂_ A _ ^″^ in the Γ̂_ A _ operator (and thus satisfy a spinless version of eq; i.e., we calculate the spatial part of the electrons in the body frame while the electronic spins are in the space-fixed frame); the relevant Γ̂_ A _ operator satisfies −iℏ∑_ A _ X _ A _ × Γ̂_ A _ + L̂ ^ e ^ = 0 and conserves the total angular momentum excluding spin.? When the nuclei are rotating around the origin (we set the center of mass as the origin) with an angular velocity ω, the coupling term produced by the Γ̂ becomes

$$[eqn]$$

which is exactly the electronic Coriolis term in a rotating noninertial frame of an angular velocity ω. Previously, this result has been shown to hold generally for any type of nuclear motion, including both rotations and vibrations, in the nonlocal limit ?,? (where σ → ∞ in eq); here, we emphasize that this result holds regardless of the choice of the locality parameter in the case of pure rotations. Moreover, further including the electron spin rotation factor Γ̂_ A _ ^‴^ so as to satisfy the original eq (corresponding to a calculation of the spatial and spin components of electrons in the body frame), we find a coupling term

$$[eqn]$$

which is a spin Coriolis term.

Throughout this article, we will include all Γ̂′, Γ̂″and Γ̂‴in Γ̂ unless otherwise noted. As we will show below, the features above enable us to describe the Λ-doubling effect – an effect traditionally considered inaccessible directly within a BO approximation.

Phase Space in the Case of Degeneracy

3

One of the most interesting facets of phase space electronic structure theory is the behavior of potential energy surfaces in the presence of degeneracy. In this context, a phase space approach to electronic structure can break symmetries that are protected within BO theory. In particular, with the BO framework, every radical molecule (with one unpaired electron) must be part of a degenerate electronic subspace (usually a doublet) according to Kramers’ theorem.? That being said, however, according to a phase space description, the two electronic states may interact differently with nuclear motion (that breaks time reversal symmetry) and lose their energetic degeneracy. As a result, even for the smallest molecule, one finds multiple minima in P-space; see ref ?.

Multiple Minima on the Potential Energy Surfaces

3.1

In Figure, we plot the two lowest phase space surfaces from solving eq for the molecule NO in a STO-3G basis. On the x-axis, our order parameter is the canonical nuclear rotation angular momentum L = ∑_ A _ X _ A _ × P _ A _. Note that, whereas we find a doubly degenerate ground state PES under the BO framework, PS theory now splits into two shifted, approximately parabolic phase space PESs. As will be discussed below in section, the new broken-symmetry eigenstates have different reflection symmetries in the case of diatomics. Note, however, that the two PESs are still degenerate at zero canonical nuclear momentum and have a mirror symmetry with respect to L = 0, all due to time reversal symmetry in the combined electronic and nuclear space.

Potential energy surfaces of the ground and first excited states of NO from the phase space Hamiltonian, plotted versus the canonical nuclear angular momentum L. The blue and orange curves correspond to states whose valence Π molecular orbitals are approximately 12[|Px⟩|↓x⟩−i|Py⟩|↑x⟩] and 12[−|Px⟩|↑x⟩+i|Py⟩|↓x⟩] , respectively. The nuclear coupling breaks the Kramers’ degeneracy of the two PESs at a given finite canonical nuclear rotation angular momentum L and form a double well with minima at L/ℏ ≈ ±0.0022.*

In order to understand the physical meaning of the two split PESs and of their unique minima along the canonical nuclear rotational angular momentum L, one must recognize that the variable P, which would appear to be the nuclear linear momentum within the Born–Oppenheimer framework, represents instead the combined nuclear and electronic momentum.? This representation can be shown by applying the total momentum operator P̂ + p̂ = ∑_ A _ P̂ _ A _ + ∑_ i _ p̂ _ i _ to a coupled electronic and nuclear wavefunction Ψ­(X, r) = ∑_ n ψ n (X) ϕ n _(X, r) in the BO formalism

$$[eqn]$$

Here, we have used the fact that the electronic wavefunction ϕ_ n _(X, r) depends only on the relative electron–nucleus positions,

$$[eqn]$$

Thus, the nuclear momentum evaluated on the nuclear wavefunction ∑_ n _ ⟨ψ_ n |P̂|ψ n _⟩ in fact represents the total momentum of the electrons and nuclei ⟨Ψ|(P̂ + p̂)|Ψ⟩. In fact, when one rewrites the Schrodinger equation in the adiabatic basis and replaces

$$[eqn]$$

where d̂ _ A _ is the derivative coupling operator, one can intuitively think of d̂ _ A _ as incorporating the electronic momentum. ?,? Thus, while P̂ _ A _ represents the total electronic + nuclear momentum, P̂ _ A _ – iℏd̂ _ A _ represents the nuclear momentum. Analogous statements apply to angular momentum as well as linear momentum: the canonical nuclear rotational angular momentum L in fact represents the total angular momentum J = R + L ^ e ^ + S ^ e ^ where R is the nuclear rotational angular momentum, L ^ e ^ and S ^ e ^ are the electronic orbital and spin angular momentum.?

The facts above form the basis for phase space electronic structure theory. As described above, within PS theory, we regard P̂ _ A _ as a classical object P _ A , but we replace d̂ _ A _ with Γ̂ A _; thus, at the end of the day, the semiclassical quantity

$$[eqn]$$

represents the nuclear momentum and is known as the kinetic linear momentum. At this point, the physical meaning of the multiple minima of E PS becomes clear. Namely, energy minima can arise only when the kinetic momentum of each nucleus is zero:

$$[eqn]$$

Thus, at the minimum, the value of P _ A _ ^min^ is dictated by the different (nonzero) electronic momenta that can be captured by the different electronic states in a degenerate subspace, P _ A _ ^min^ = iℏ⟨Γ̂_ A ⟩. Again, analogous statements hold for angular momentum, L ^min^ _ A _ = iℏX _ A _ × ⟨Γ̂ A _⟩, which is more rich because these angular momenta couple with each other and with the environment as magnetic dipoles.

Extracting Rotational Energies

3.2

Within BO theory, the simplest means of extracting vibrational and rotational energies is very well explored. ?−? ? Namely, one finds the local minimum geometry X 0 of E _ BO _(X 0) in nuclear configuration space, one separates rotations from translations by looking for zeros of the Hessian, one makes a rigid rotor approximation for the rotations, and one makes a harmonic approximation for the vibrations.

How should one accomplish this task within a phase space perspective, given the existence of multiple minima in E _ PS _(X, P)? With reference to eq above, the simplest means is to diagonalize the electronic Hamiltonian, and then both make a harmonic approximation for Ĥ el around X 0 and fix the value of Γ̂ at the Born–Oppenheimer equilibrium geometry X 0 (with P = 0):

$$[eqn]$$

Equation is obviously quite reminiscent of more advanced theories of molecular dynamics that separate rotations and vibrations; for example, see ref ?. For this article (with a focus on the spectrum of NO), we will focus exclusively on rotations; vibrations and vibrational–rotational coupling will be treated within a forthcoming article. With the interpretation above, the rotational kinetic energy of a molecular system is written as

$$[eqn]$$

so that one might conjecture a natural nuclear-electronic coupling of the form

$$[eqn]$$

The operator evaluated in eq is often termed the sum of the L-uncoupling and S-uncoupling operators.? In the context of ground state diatomics, as we will show below (section), such an electron–nuclear coupling exactly gives rise to the Λ-doubling effect.

Λ-Doubling and NO

4

Strictly speaking (i.e., without the Born–Oppenheimer approximation), for a diatomic molecule, as for any molecule, the only good quantum numbers are J, the total angular momentum quantum number, and its lab-frame z-axis projection M; each energy eigenstate can be simultaneously labeled by one J and one M. The Λ-doubling effect occurs when there are two degenerate electronic configurations with positive and negative angular momentum projections Λ along the internuclear axis according to the BO electronic Hamiltonian. As will be seen below, this degeneracy is lost due to electron-rotation coupling, leading to a splitting in the rotational spectrum.

Standard View of Λ-Doubling

4.1

When the nuclear rotational coupling and spin–orbit coupling (SOC) are small compared to electrostatic interactions, these objects are usually treated as perturbations on top of the nonrelativistic BO Hamiltonian H 0, i.e.

$$[eqn]$$

In the case of a relatively strong SOC and weak rotational coupling (Hund’s case (a)), one expands the wavefunction in terms of a specific set of the eigenstates of the nonrelativistic BO Hamiltonian, where each function is labeled by the additional approximately good quantum numbers: the electronic orbital angular momentum projection along the internuclear axis Λ, the electronic spin angular momentum S, the latter’s projection along the internuclear axis Σ, and Ω = Λ

• Σ which is also the total J projection along the internuclear axis (because any nuclear rotation angular momentum R has no projection along the internuclear axis). Within the BO approximation, each basis function is then a product of the electronic wavefunction, solved in the molecular frame defined by the nuclei, and a rotational wavefunction that describes the orientations of the molecule, including both nuclei and electrons, within the lab frame,

$$[eqn]$$

where J is the total electronic plus nuclear angular momentum and M its projection along the lab-frame z-axis, and we define

$$[eqn]$$

where D is the Wigner D-matrix.? Note that J and M are the only good quantum numbers for an isolated system that has rotational invariance. Additionally included in eq is an energy label n (n = X, A, B, ...) and a separable vibrational wavefunction |ν⟩, but these parts will be omitted in the following to simplify the discussion.

In the case of the NO molecule, the lowest energy states are primarily made up of the ^2^Π_1/2_ (Λ = ±1, S = 1/2, Σ = ∓1/2, Ω = ±1/2) manifold; these states are fairly well separated from the higher energy ^2^Π_3/2_ (Λ = ±1, S = 1/2, Σ = ±1/2, Ω = ±3/2) states due to the strong SOC. As a heteronuclear diatomic molecule with C ∞v _ point group symmetry, NO has a “parity symmetry,” or more precisely, a reflection symmetry with respect to any plane that contains its internuclear axis. Any eigenstate should preserve such symmetry (or in the case of degeneracy, one can choose eigenstates that preserve such symmetry), so within the ^2^Π_1/2 space, the basis functions come in pairs to form parity eigenstates

$$[eqn]$$

Here, p ^±^ refers to the parity eigenstate with total parity ±(−1)^ J−S ^ for electronic Π terms. For the half-integer J values and $[eqn]$ considered here, the |p ^+^⟩ state with parity $[eqn]$ is further labeled e, while the |p ^–^⟩ state with parity $[eqn]$ is labeled f. A′ and A″ have also been used for labeling based on the reflection symmetry of the electronic spatial coordinates.?

Now, if we consider the effect of nuclear rotations as in eq, the ^2^Π manifold interacts with higher-lying ^2^Σ states. (Importantly, the Σ label here signifies Λ = 0; in an unfortunate aspect of small molecule notation, Σ has two different meanings and the Σ here is different from the Σ in eq which refers to the internuclear axis projection of the spin angular momentum.) For NO, the interaction with high-lying ^2^Σ states involves the off-diagonal L-uncoupling operator,?

$$[eqn]$$

where B = 1/2I is the rotational constant in the equilibrium geometry (divided by the factor of ℏ^2^), and the off-diagonal SOC operator, A/2­(L̂ _ e _ ^+^ Ŝ _ e _ ^–^ + L̂ _ e _ ^–^ Ŝ _ e _ ^+^) where A is the spin–orbit constant (also divided by the factor of ℏ^2^). Altogether, this approach generates a sum-over-states first-order corrected ground-state wavefunction

$$[eqn]$$

where |0⟩ is a BO ground state |^2^Π_1/2_, J, M, p ^±^⟩ and each |i⟩ is an excited state |^2^Σ_1/2_, J, M, p ^±^⟩ with the same parity as |0⟩. The final sum-over-states second-order energy correction is

$$[eqn]$$

Evaluating eq is difficult as one needs to include enough states for convergence (which is many). Finally, note that the S-uncoupling operator in eq also couples together the ^2^Π_1/2_ and ^2^Π_3/2_ manifolds for J ≥ 3/2. A careful diagonalization within the combined space of ^2^Π_1/2_ and ^2^Π_3/2_ manifolds ?,? yields the following result for the Λ-splitting of the ^2^Π_1/2_ manifold

$$[eqn]$$

where Y = A/B, $[eqn]$ , and

$$[eqn]$$ $$[eqn]$$

which are calculated by explicitly summing over all important ^2^Σ excited states.

This concludes our brief review of the standard, well-established formal theory of Λ-doubling. ?,?,? Note that below we will compare phase space theory directly to experimental data (rather than to eqs and ?), but the analysis above will be helpful for reasons of interpretation. Note also that the derivation of eq ignores the parity dependency of the energy denominator, which would require a much more involved analysis to include.?

A 1D Model Hamiltonian Compatible with Phase

Space Electronic Structure Theory

4.2

Equations–? above represent the standard, perturbative approach to Λ-doubling that one finds in the literature. We would now like to show that this result can also be recapitulated quite naturally (and without a sum over states) using PS theory; the PS approach is not based on perturbation theory.

In line with section above, the simplest model Hamiltonian that can capture coupled nuclear rotation-electronic effects is a standard rigid rotor that can be rotated with angle θ around one axis (hereby the name “1D” model). We let θ be the rotational angle of interest and we imagine that there are two electronic states of interest (as in the case of NO described above in Figure) with ±α total electronic angular momentum. In the spirit of eq, the kinetic momentum is $[eqn]$ , which leads to the following (trivially separable) Schrodinger equation:

$$[eqn]$$

Given the need for periodicity, the solutions to this Schrodinger equation are ψ_ n _(θ) = exp­(inθ) with energies

$$[eqn]$$

where n takes on half-integer numbers for systems with odd numbers of electrons and integer numbers for systems with even numbers of electrons. The difference between these energetic solutions is

$$[eqn]$$

This expression for 1D Λ-splitting can be further substantiated by comparing with the perturbative (non-phase space) result for a rotation J in 1D. Without loss of generality, we look at a finite J along the a-axis. We label the eigenstates of the standard electronic Hamiltonian (without SOC) as |0⟩, |1⟩, |2⟩, ...; in cases of degeneracy between any two eigenstates, these states are rotated to form eigenstates with reflection symmetry. The Λ-splitting due to the second-order coupling between the SOC and L-uncoupling, which dominates in the ^2^Π_1/2_ manifold in molecules with strong SOC like NO, according to eq, is

$$[eqn]$$

where |Ψ_SOC_ ^(1)^⟩ is the first-order corrected wavefunction due to SOC when no nuclear rotation is considered

$$[eqn]$$

Within a perturbative treatment, |Ψ_SOC_ ^(1)^⟩ should be symmetric under reflection. To a good approximation, we can also write eq above as

$$[eqn]$$

where ⟨L̂ _ a _ ^ e ^⟩_ J _ a =0 ≈ ⟨Ψ_SOC_ ^(1)^|L̂ _ a _ ^ e ^|Ψ_SOC_ ^(1)^⟩ is the expectation value of L̂ _ a _ ^ e ^ with respect to a relativistic (SOC-including) BO eigenstate (again without any nuclear rotation). Above we have used the fact that ⟨0|L̂ _ a _ ^ e ^|0⟩ = 0.

Equation matches eq above if we set α = ⟨L̂ _ a _ ^ e ^⟩_ J _ a =0.

A 2D Model Hamiltonian Compatible with Phase

Space Electronic Structure Theory

4.3

Next, let us imagine using phase space electronic structure theory to describe Λ-doubling in a full three-dimensional context, beyond the 1D model above (which is trivially solvable). As we will see, a multidimensional model cannot be solved analytically and requires either a large matrix diagonalization or approximations.

To begin our multidimensional analysis, let c denote the molecular axis of the diatomic molecule (from N to O). Due to the C _∞v _ symmetry of the NO molecule around the internuclear c-axis, any choice of an axis perpendicular to c is equivalent. We will pick the a- and b-axes associated with χ = 0, ?,? where χ is the last of the Euler angles (ϕ, θ, χ) that describes the molecular orientation in the z–y–z convention. As can be found in many textbooks, ?,? the quantum mechanical Hamiltonian for the nuclear kinetic energy associated with rotations of a linear molecule

$$[eqn]$$

where I = I _ a _ = I _ b _, is isomorphic to the Hamiltonian

$$[eqn]$$

The Hamiltonian Ĥ rot ^iso^ has one more degree of freedom than the true Hamiltonian in eq; namely, χ remains as an independent variable in the definition of Ĵ. For an isomophormism to hold, we must diagonalize eq within a basis that imposes the condition that J _ c _ = J _ c _ ^ e ^, so as to remove all extraneous solutions; for more details, see refs ? and ?.

In the context of phase space electronic structure theory, we will use a variant of the above isomorphic Hamiltonian where the nuclei are classical. It is then clear that the ground state phase space potential energy surfaces for a diatomic can be constructed by minimizing the following Hamiltonian:

$$[eqn]$$

or in short hand (remembering that I _ c _ = 0):

$$[eqn]$$

Now, the first and third terms are merely constant energy shifts for the lowest two states, and thus we focus on only the coupling term. Given that we are working with a system with two electronic states, we will use the form of the Pauli matrices (σ_ a , σ b , σ c _) for a spin 1/2 system and write

$$[eqn]$$

where

$$[eqn]$$

and ∥Ĵ∥ denotes the 2-norm (square root of the sum of squares) of operator Ĵ. At this point, we can rewrite the coupling term as a matrix in the bases of |Λ = 1, S = 1/2, Σ = −1/2⟩ and |Λ = −1, S = 1/2, Σ = 1/2⟩ electronic states. Because we find the matrices of Ĵ _ a _ ^ e ^ and of Ĵ _ b _ ^ e ^ are real off-diagonal and complex off-diagonal in such bases with a certain phase choice, the coupling term becomes

$$[eqn]$$

and evaluate the energy splitting by diagonalizing the Hamiltonian

$$[eqn]$$

In the last step, when evaluating J ^2^, we have requantized J so as to be the total angular momentum operator Ĵ. The energy splitting becomes

$$[eqn]$$

Equation is our most reliable result (fulfilling the expectation of eq) and can be checked against experiment; in order to model Λ-doubling, in principle all we require are the matrices for Ĵ _ a _ ^ e ^ and Ĵ _ b _ ^ e ^, which can be calculated in several ways through an ab initio phase space framework.

Ab Initio Results

5

Let us now address the necessary calculations. Below, we fixed the bond length of NO to be the experimental value? 1.154 Å and ignored any length dependence on vibrational or centrifugal distortion. Except where otherwise noted (see below), we treat SOC with the one-electron Breit–Pauli Hamiltonian? (not including the two-electron spin-same-orbit and spin-other-orbit coupling). In our study of NO, as a sanity check, the first thing we calculated was the energy splitting between the ^2^Π_1/2_ and ^2^Π_3/2_ manifold at L = 0 due to SOC. For these calculations, we have run full configuration interaction (FCI) calculations for absolute certainty. The low-energy eigenstates of NO, illustrated in Figure, are split by strong SOC into the ground state ^2^Π_1/2_ manifold and a ^2^Π_3/2_ manifold that lies 143 cm^–1^ higher, in reasonable agreement with the experimental separation? of 123 cm^–1^ with about 16% overestimation. Note that this splitting is far bigger than the Λ-splitting; see Figure.

Energy diagram of the low-energy space of NO. Energy splittings are not drawn to scale.

Next, let us address the Λ-splitting. In line with the discussion above, we will focus on the splitting of the ground doublet in the ^2^Π_1/2_ manifold. As mentioned above, if we wish to use the approach in section, there are two ways to extract the key objects, L ^ e ^ and S ^ e ^ (or really just L ^ e ^ in practice), from a phase space calculation (and see Table below.).

• Method I: On the one hand, one can evaluate the two complete tensors in the low-energy manifold at L = 0; i.e., one computes (L̂ ^ e ^ + Ŝ ^ e ^)|_ L=0_, as suggested in eq above and discussed also in eq. Thereafter, one can set $[eqn]$ . Alternatively, for the NO problem, if one rotates to a basis with reflectional symmetry in the b–c plane, L̂ _ a _ ^ e ^+ Ŝ _ a _ ^ e ^ will be diagonal and so we can set α = |⟨L̂ _ a _ ^ e ^ + Ŝ _ a _ ^ e ^⟩|_ L=0_, where the bracket ⟨⟩ indicates evaluation on the ground state.
• Method II: On the other hand, if one seeks to intuitively extract rotational energies exclusively from the ground state phase space PES (without any further electronic calculations), one can alternatively find the PES minimum in phase space (L γ ^min^) as a function of rotation in each γ-direction. Because of eqs and ? above, one is guaranteed that

$$[eqn]$$

Thus, given that L̂ _ a _ ^ e ^+ Ŝ _ a _ ^ e ^ will be diagonal in the ground state with the proper symmetry (as noted for Method I above), we can set α = L γ ^min^. Note that, for a nonlinear molecule with multiple different moments of inertia, this method would require several different optimizations to recover all three rotational directions. Note also that, in practice, the phase space minimum L γ ^min^ is very small and close to zero for NO, so that ⟨L̂ γ ^ e ^ + Ŝ γ ^ e ^⟩_ L γ ^min^ _ ≈ ⟨L̂ γ ^ e ^ + Ŝ γ ^ e ^⟩_ L=0_.

**1: Coupling Parameter α Estimated in Three Ways: Method I: from ⟨L̂

a

e

• Ŝ

a

e ⟩ Evaluated on the Parity Eigenstates at L = 0; Method II: from L

a at the Minimum of the PES; Method III: from Energy Splitting ΔE at L

a = 1/2ℏ**

Finally, note that there is also a third and direct semiclassical approach.

• Method III: We can simply run a PS calculation with L = 1/2ℏ and obtain the Λ-splitting associated with a rotation around one axis directly from the energy gap between two PESs (Figure). α can then be fit from the energy gap using eq.

In practice, we find that all three of these approaches yield very similar results, though our results for Method II are perhaps the least accurate because of the need for a high resolution, accurate potential energy minimum. Note also that, for the NO molecule, the matrix element α is exclusively from the L̂ ^ e ^ operator; the spin component Ŝ ^ e ^ is negligible.

Estimate of Coupling From PS calculations

5.1

During the course of our investigation, our initial approach was to calculate the Λ-splitting directly by running a phase space electronic structure problem for a value of L = ℏ/2 (unless otherwise noted, and then directly measuring the gap). We performed several such Method III calculations:

• (i)We ran an FCI calculation in the minimal basis STO-3G.
• (ii)We solved the same Hamiltonian using the density matrix renormalization group (DMRG) solver? with enough bond dimension to converge to the FCI accuracy in a slightly larger basis 6-31G, after benchmarking against FCI in the STO-3G basis. For these calculations, the ground and first three excited states were solved using a state-averaged DMRG approach.? Bond dimensions of 500 and 3000 were used for the STO-3G and 6-31G basis sets, respectively, to converge the first excitation energy to about 10^–12^ and 10^–9^ Ha precision. Unlike the calculations for (i) and (iii), for better numerical accuracy when evaluating a small energy gap, we invoked a linear approximation for the phase space energies and ran calculations with L = 3ℏ instead of ℏ/2 and rescaled the final energy splitting by 1/6; the assumption that the Λ-splitting scales approximately linearly with J for small J follows from the parabolic nature of Figure.
• (iii)In order to reach a full basis limit, we developed and implemented a constrained Hartree–Fock (cHF) ?−? ? ? approach which avoids the incorrect broken symmetry solutions with S _ c _ ^ e ^ = ±1/2 and L _ c _ ^ e ^ = ∓1 as found by standard HF theory. Here, our constraint penalizes finite S ^ e ^ vectors by adding an energy penalty term ΔE = −λ|S ^ e ^|^2^ and therefore yields the correct parity eigenstates with S _ c _ ^ e ^, L _ c _ ^ e ^ ≈ 0. We applied this cHF method first with λ = 0.1, and then relaxed the converged solution with λ = 0 until convergence again. To obtain the first excited state with the opposite parity, we used the time reversal of the ground state electronic wavefunction as the initial guess, converged the cHF calculation with λ = 0.1, and then relaxed the solution with λ = 0.

All of our Method III results are shown in Figure. In the minimal basis, the first excitation energy from DMRG agrees with the FCI result within 0.5% of error, confirming that the linear approximation of Λ-splitting holds quite accurately and that the DMRG result in the 6-31G basis provides an accurate estimate of the FCI result. Relative to the FCI and DMRG benchmarks, in both STO-3G and 6-31G basis sets, cHF consistently overestimates the energy splitting by about 20%. Nevertheless, cHF provides an efficient means to approach the complete basis set limit. The cHF results show that the basis-set incompleteness error in the minimal basis is small. Therefore, it is not unreasonable to expect that the minimal basis FCI result provides a good estimate of the FCI energy in the complete-basis limit. Thus, below we will use the FCI energy splitting in the STO-3G basis to estimate the α parameter in a 1D version of the model Hamiltonian eq (“Method III” in Table).

Λ-Splitting between the lowest e and f states and the associated α parameter for a nuclear rotation J = 1/2ℏ around one perpendicular axis. The results are calculated using the exact FCI solver (blue circle) in the STO-3G basis, the DMRG solver (purple star) which converges to the FCI accuracy in the STO-3G and 6-31G bases, and an approximate constrained Hartree–Fock solver in larger bases, from STO-3G, 6-31G, cc-pVDZ, aug-ccpVDZ, to aug-ccpVTZ (red).

As a complement to the Method III results, we also followed Method II above and scanned the FCI phase space PES looking for a minimum so as to determine the α parameter based on eq. We find that the PES minima are at L _ a _ = ±0.0022, and therefore, α = ∥Ĵ _ a _ ^ e ^∥ = 0.0022 (“Method II” in Table). Lastly, we also followed Method I above and computed the matrices of L̂ ^ e ^ and Ŝ ^ e ^ for the two lowest energy eigenstates at L = 0. This yields α = ∥Ĵ _ a _ ^ e ^∥ = 0.0018 (“Method I” in Table). Thus, in total, we can estimate the splitting in three ways and they all yield consistent results.

Before concluding this section, we note that using Γ̂ = Γ̂′ + Γ̂″is enough; the absolute coupling is largely unchanged by adding Γ̂‴. Both Γ̂′and Γ̂″ contribute to the energy splitting.

The Splitting of the 2Π1/2 Manifold as a Function of J

5.2

Having shown that we can evaluate the Λ-splitting between the e/f states for J = 1/2 for a “clamped” rotation around one axis in several different consistently ways, we will now use the value of α found above and plug into eq and make predictions for rotations in the full two dimensions and for larger J values. In principle, one could use the 1D model in eq, but we will use the 2D model in eq (as one can always argue about the validity of a 1D rotational model). That being said, in practice, the only difference is a factor of J versus $[eqn]$ .

In Figure, we compare the Λ-splitting hereby estimated with the experimental measurement for small half-integer J. Here, we plot FCI results using α estimated from the energy splitting in the minimal basis (Method III). When we only include the 1-electron SOC, these results are expected to moderately overestimate the splitting due to the overestimated SOC discussed at the beginning of section. We also plot cHF results in the complete basis set limit, which is also expected to overestimate the experimental answer, due to both the overestimation of SOC and the mean-field approximation. According to Figure, we find that if we include only the 1-electron SOC, the FCI results (“FCI (1e SOC)” in Figure) overestimates the splitting by 7–25% and the cHF results (“cHF (1e SOC)”) overestimate the splitting by 36–59%. If we further include the 2-electron SOC under the spin–orbit mean-field approximation, the FCI results (“FCI (1e+2e SOC)”) underestimates the splitting by 34–23% and the cHF results (“cHF (1e+2e SOC)”) deviate by only 4–9%. It would appear that, if we include the full SOC tensor and work in a big enough basis, our semiclassical approach can be nearly quantitative–though further testing (and comparisons with fully quantum-mechanical calculations) will be necessary to say much more.

Energy splitting between the e and f states as a function of total angular momentum J from the 2D model parametrized with α derived from FCI in the STO-3G basis (blue) and from cHF in the aug-ccpVTZ basis (red), both using Method III, in comparison to experiments (black). Dashed lines include only the 1-electron SOC, while solid lines include both 1-electron and 2-electron SOC (with a spin–orbit mean-field approximation , ).

Discussion

6

Kramers’ Theorem

6.1

At this point, we have shown that a phase space interpretation of electronic structure theory automatically yields a meaningful Λ-splitting. Before we interpret this result, a few words are appropriate with regard to how our finding above relates to Kramers’ theorem. Kramers’ theorem states that for any time-reversal-invariant system with an odd number of electrons (such as NO), if we ignore nuclear spin, then every energy eigenstate is at least doubly degenerate, as each eigenstate |ψ⟩ can be paired with its orthogonal time-reversal pair T|ψ⟩. This statement holds when considering either (i) the electronic Hamiltonian within BO theory or (ii) the total electronic + nuclear Hamiltonian. For our purposes [focusing on (ii)], it is crucial to emphasize that the parity states (eq) that are associated with Λ-splitting contain both electronic and rotational components. In what follows below, we will follow the phase convention of the time-reversal operator T used by Condon and Shortley? and by Hougen?

$$[eqn]$$

where the second quantum number denotes either lab-z or body-c axis projection of J. The same convention is used for |S, Σ⟩ and |L, Λ⟩. Of course, in tandem with eq, one must have

$$[eqn]$$

Let us now consider the implications of Kramers’ theorem for the wavefunctions of the low-energy manifold of NO. If we apply the time reversal operator to the electronic component, focusing on the dominant ^2^Π_1/2_ contribution to the ground state, one finds

$$[eqn]$$

Next, if one applies time reversal to the rotational component, one finds

$$[eqn]$$

This result follows from eqs and ? and the well-known identities (for a half integer M):

$$[eqn]$$ $$[eqn]$$

As a result, if one applies the time reversal operator to the total (nuclear + electronic wavefunction in eq), one obtains

$$[eqn]$$ $$[eqn]$$

The same conclusion holds for the ^2^Σ components of the low-energy state wavefunctions: the time reversal yields a state with the same parity but opposite M, apart from a phase factor. In the end, the time-reversed state is orthogonal to and degenerate with the original state; after all, the energy must be invariant to the sign of M (which is really a statement about time-reversal symmetry). Of course, without an external field, there can be further symmetry. For example, without an external electric field, rotational invariance would imply not just 2-fold degeneracy, but actually 2J + 1-fold degeneracy. That being said, in the presence of an external electric field, 2-fold degeneracy still remains (between +M and −M) for an odd spin system. In order to violate Kramers’ theorem, one must break time-reversal symmetry by, e.g., applying a magnetic field.

A Phase Space Visualization of the e/f States

6.2

Finally, how should we interpret the result above, where we have found a reasonably accurate estimate of Λ-doubling with a quite simple calculation? The phase space framework above paints a very intuitive picture for visualizing the e and f parity (total electronic and nuclear) wavefunctions in eq. To establish such a picture, let us explore the phase space PES for a rotation around the a axis, e.g., L _ a _ > 0. We rewrite the dominant components of the low-energy eigenstate wavefunctions eq as

$$[eqn]$$

where

$$[eqn]$$

are the two nuclear rotational states whose Ĵ _ a _ expectation values are

$$[eqn]$$

and

$$[eqn]$$

are the two electronic eigenstates that are symmetric with respect to the reflection through the plane bc, up to a phase factor. According to eq, when the molecule has a positive angular momentum around the a-axis, the electronic ground state is dominantly |ζ^+^⟩; when the molecule has a negative angular momentum around the a-axis, the electronic ground state is dominantly |ζ^–^⟩. This fact implies that, whereas the crossing phase space PESs in Figure (copied in Figurea) are defined (and labeled) by their constant electronic character, the observed e/f states arise and are defined (and labeled) by their relative energies as in Figureb. In the latter case, the e/f states are not so much discrete states but rather arise as part of a continuous potential energy surface that must be requantized (as in section above).

a) Our original view (see Figure ) of potential energy surfaces of the ground and first excited states of NO according to a phase space electronic structure Hamiltonian as a function of the phase space canonical nuclear angular momentum L. b) The same data now labeled differently: The green and red curves can be understood as the e and f parity total electronic and nuclear states, now understood as functions of total angular momentum J (represented by the phase space L). The e/f states are the states that are observed spectroscopically.

From a quantum chemistry perspective, if one considers the two BO electronic states of NO in one orientation as the degenerate “diabatic states” of interest, then the nuclear motion breaks this degeneracy, leading to the ”adiabatic states” of interest in Figure. However, the e/f states have identical energies at L = 0 because of the form of the splitting in eq, so that there are two different labeling schemes – Figurea versus Figureb. In a sense, one can argue that NO has a “conical intersection” ?−? ?

in momentum space, which provides a new interpretation of Λ-doubling (and a different view of the magnetic monopole derivative couplings discussed in refs ? and ?). Note that for the case of a diatomic molecule (as in eq), the conical intersection in phase space has codimension two (because only L _ a _ and L _ b _ enter the rotational Hamiltonian); more generally, in line with ref ?, we would expect the conical intersection around L = 0 to have codimension three. Note also that our analysis above of eq is effectively quantitative insofar as a purely quantum calculation using $[eqn]$ will necessarily lead to an energy splitting proportional to a factor $[eqn]$ , which is in agreement with eq above.

The 2Π3/2 Manifold

6.3

During the review process for this manuscript, a reviewer asked a very interesting question related to the scaling above. Namely, the reviewer asked us if a phase space method would recover the correct form for the Λ-splitting within the ^2^Π_3/2_ manifold of NO. The question above is interesting for two reasons. First, the origin of the Λ-splitting for the ^2^Π_1/2_ manifold is different than the origin of Λ-splitting for the ^2^Π_3/2_ manifold. For the former, the dominant contribution is the double perturbation of H SOC ^(1)^ and the Coriolis potential H ^(1)^ Cor. By contrast, for the latter, the lambda splitting has a different origin; here, the effect arises exclusively from the Coriolis potential; i.e., a second-order perturbative treatment in H Cor ^(1)^ is required. As a result of this difference, the lambda splitting scales as J for the former (see Figure) and cubically for the latter. Obviously, one must wonder if a phase space approach can recover this difference? After all, if one were to follow Methods (i) and (ii) in section, i.e., evaluating α at P = 0, and then diagonalize a matrix as in section, one would presumably find a splitting that is always linear in J.

Second, a keen observer will note that, in order to recover the Λ-doubling results in Figure above, all that is required is that (i) one include SOC within a BO calculation and then (ii) one couple together the nuclear and electronic angular momentum through the L̂ ^ e ^ operator, as in eq above. In other words, one can in fact recover the Λ-doubling splitting of the ^2^Π_1/2_ manifold without really ever invoking the phase space electronic structure Hamiltonian in section; if one wishes to avoid exploring phase space potential energy surfaces, one can model Λ-doubling directly with the relativistic BO states by including L̂ ^ e ^ alone (and without building the Γ̂ operator). That being said, just considering the scaling argument above, such an approach is clearly not appropriate for the ^2^Π_3/2_ manifold. For the latter case, the only possibility is to invoke Method (iii), fully explore the phase space potential energy surfaces, and then investigate the energy gap at L = nℏ/2. As a result, such an approach will inevitably require that we build Γ̂ for P ≠ 0 (i.e., not just L̂ ^ e ^) and whether or not PS can recover the correct Λ-splitting for the ^2^Π_3/2_ manifold is in fact a very strenuous test of the method.

As shown in Figure, phase space theory does yield an accurate Λ-splitting as compared with experiment. Note that, in order to resolve the tiny absolute energy scale (on the order of 10^–10^ Hartree), we required a tight convergence threshold for the diagonalization. In Figure, we also fit our data to an analytical energy expression from standard perturbation theory (valid under slow rotation and strong SOC),? namely, $[eqn]$ cm^–1^ (see Figure, blue dashed curve). Our results show a good fit, with a R ^2^ value of 0.9995. Altogether, the data in Figure again demonstrate the strong accuracy that can be found by solving electronic structure problems within a phase space framework.

Energy splitting between the e and f states of the 2Π3/2 manifold as a function of total angular momentum J. For this data, unlike Figure , it is well-known that the energy splitting requires a cubic (rather than linear) fit as a function of J. We compare phase space predictions, using a model parameter α derived from FCI with 1-electron SOC in the STO-3G basis with Method III (blue dots) versus experiment (black triangles). We further fit the “phase space” data to an analytic expression from standard perturbation theory, ΔE=k(J+12)(J+32)(J−12) , where k is a fitted constant (dashed lines). Note that, as in Figure above, our phase space predictions are quite accurate, which is another strong confirmation of a PS approach.

Conclusions and Outlook

7

We have demonstrated that a simple phase space perspective on electronic structure theory provides a natural framework to capture the Λ-doubling in diatomic molecules. Thus, in line with previous discussions, ?,? we assert that the same beyond-BO physics that underlies one of the smallest microscopic diatomic phenomena is also responsible for the macroscopic Einstein–de Haas effect ?,?,? and, perhaps, chiral-induced spin selectivity. ?−? ? Interestingly, the major result here has been the nonadiabatic splitting ΔE between different states of different angular momentum, where dynamically ℏ/ΔE corresponds to the tunneling time between degenerate states. In the context of the Einstein–de Haas effect, one is less interested in parity eigenstates (e.g., e/f), but rather in a basis with maximally aligned and opposite spin directions along the magnetic field direction. Nevertheless, the approach above should be able to isolate the tunneling time for a spin (or a collection of spins) to change orientation. In other words, whereas most treatments of Einstein–de Haas phenomena use angular momentum conservation to predict the steady state nuclear rotational angular momentum that arises upon a total change in electronic spin angular momentum, the present approach based on phase space electronic structure theory will be able to go further and calculate transient effects: How long does it take for the spins to align? For the metal to reach a steady state of rotation? And of course, in the long run, how much energy is lost if we allow for frictional losses from nuclear vibrations?

With regard to this last point, although in this work we have restricted our attention to electron coupling to purely rotational nuclear motion, the most interesting physics undoubtedly lies ahead as we consider larger molecules and materials, and will need to model both vibrations and rotationsand their coupled effect on electrons ?−? ? all on an equal footing. That being said, in the future, if we seek to model larger molecules or materials, it will be essential to efficiently build the Γ̂ operators, which will require some numerical tricks on our end and a more efficient implementation. Our results here were calculated in PySCF, ?,? and our group also has a preliminary code in the Q-Chem package.? Our strong feeling is that, given the success of the current manuscript in tandem with previous results in VCD? and ROA,? we will soon be equipped to explore a very broad set of questions: what are the roles of chiral phonons in determining spin-dependent electron transfer rates?? How should we best quantitatively describe spin relaxation?? Can we definitively prove the nature of the fundamental mechanisms behind different magnetic field effects?? Lastly, this manuscript has suggested that whenever the ground state has a degeneracy built on different electronic angular momenta, nuclear motion will lead to a small splitting. As such, might not the equations above (especially eq, which describes a gap opening due to off-diagonal electron–phonon coupling and is analogous to a Su–Schrieffer–Heeger (SSH)/Peierls Hamiltonian in the case of electron hopping ?,? ) also apply to a metal and be a starting point for superconductivity calculations? When we move beyond a simple diatomic molecule, we can expect the number and the nature of interesting problems involving angular-momentum transfer to grow very rapidly.