Electronic current densities and origin-independent property densities induced by optical fields

Abstract

The interaction of a molecule with optical fields is customarily interpreted by means of induced time-dependent electric polarizabilities, magnetizabilities and mixed electric-magnetic polarizabilities. In general, these properties can be rationalized by integrals of density functions formulated in terms of induced charge and current densities. In this perspective, we focus on what has been done so far at the theoretical level, and on what can be expected to be unveiled from the topological study of suitable density functions, endowed with the fundamental requirement of origin invariance. Densities characterized by such a property can be integrated all over the configuration space to obtain electric dipole polarizability and optical rotatory power. Corresponding maps visualize domains mainly involved in the molecular response. The diagonal components of origin-independent density tensor functions that, on integration, yield corresponding electric dipole polarizability tensor of benzene, naphthalene, phenanthrene and ovalene, have been computed, confirming the ubiquitous presence of counter-polarization regions in the proximity of the atomic nuclei. They are associated with toroidal electron currents, induced by time derivative of the electric field of impinging radiation. Electron (de)localization in these systems is readily observed and estimated. The optical rotation density of the carbonyl chromophore is studied in detail. Its essential feature is the separation in quadrants of alternating sign of density about the CO bond. The presence of an extrachromophoric perturbation determines asymmetry in the extension of the quadrant distribution, thus causing optical rotation.

---

Francesco F. Summa

Francesco Ferdinando Summa graduated in chemistry at the University of Basilicata in 2018 and obtained his PhD degree from the University of Salerno in 2023, where he has now a postdoc position. His main research interest includes the topological study of static and dynamic induced current densities in molecules and the derived electric and magnetic properties.

Guglielmo Monaco

Guglielmo Monaco graduated in chemistry at the University of Naples in 1996 and obtained his PhD degree from the same University in 1999. He then moved to the University of Salerno, where he became assistant professor in 2002, and associate professor in 2020. He is interested in several aspects of computational chemistry, and particularly the study of (aromatic) molecules through the quantum mechanical current density.

Paolo Lazzeretti

Paolo Lazzeretti graduated in industrial chemistry, University of Bologna, in 1967. He became an assistant professor at the University of Modena in 1968, where he furthered his entire academic career as an associate professor of chemical spectroscopy in 1969 and professor of quantum chemistry in 1971. In 1990 he became full professor of Physical Chemistry and retired in 2013. Since 2016 he has been working as an associate professor at Istituto di Struttura della Materia, CNR Rome. He has been a visiting professor in the Universities of Lublin, College Park (Md.), USC at Los Angeles, Sheffield, Helsinki, Odense, Oulu, Paris, Buenos Aires, Bochum, Valencia. His scientific interest and activity include the formulation of theoretical methods for rationalizing and predicting response properties of molecules in the presence of electric and magnetic fields.

Riccardo Zanasi

Riccardo Zanasi graduated in chemistry from the University of Modena in 1975. After a scholarship from the CNR, in 1981 he became a researcher in Physical Chemistry at the University of Modena. In 1998 he became an associate Professor of physical chemistry at the University of Salerno, where he has been a full professor of physical chemistry since 2002. He retired at the end of 2020. His main research interests concern the theoretical definition of molecular electromagnetic properties up to the formulation of calculation methods.

---

1 Introduction

One can reasonably ask if certain regions of the electronic cloud of a molecule in the presence of electromagnetic perturbation are more involved in its reaction. In fact, the existence of moieties inducing and exalting, or diminishing, the effect of external fields and/or intramolecular electric and magnetic dipoles at the nuclei, may significantly influence the phenomenology and determine observable outcomes in response properties, e.g., electric polarizabilities, magnetizabilities and optical rotatory power.

Chemists have long studied the possibility of visualizing the influence of some substituents bestowing a particular quality, e.g., acting in a way analogous to that of a chromophore on imparting a specific color to a given compound. For instance, molecular charge density maps have customarily been employed in the physical rationalization of bonding, reactivities, and molecular properties.^1–5

Polarizability densities within simple one-electron atomic systems have been investigated by Theimer and Paul,⁶ and by Orttung and Vosooghi.⁷ An application to inert gases was considered by Oxtoby and Gelbart.⁸ Nonetheless, Sipe and Van Kranendonk analyzed some restrictions of the polarizability density idea applied to atoms and molecules in an arbitrary, nonuniform electric field.⁹

The concept of nonlocal polarizability density, first proposed by Maaskant and Oosterhoff in their theoretical study of optical rotatory power,¹⁰ has been widely adopted by Hunt et al. in several papers, dealing with van der Waals interaction,^11,12 short-range interactions on molecular dipoles, quadrupoles, and polarizabilities,¹³ infrared and Raman intensities,¹⁴ forces on nuclei in interacting molecules,¹⁵ and vibrational circular dichroism, in relation to electric-field shielding tensors.¹⁶

More particularly, the general idea of densities, functions of position r, suitable to compute molecular response properties by integration in d³r all over the molecular domain, has been analyzed in two seminal papers^17,18 by Jameson and Buckingham, who also pointed out the origin independence of some electric and magnetic property densities as a fundamental requirement for their physical usefulness.

The nuclear magnetic shielding density defined via the Biot–Savart law of classical electrodynamics¹⁹ was examined as a topical example of function,¹⁷ depending on the difference between the position of electrons and a reference nucleus, and consequently invariant of the origin of the coordinate system. Upon integration, such a function leads to the magnetic shielding of a given nucleus.

On the other hand, the density functions of electric dipole polarizability and magnetizability were shown to depend on the origin of coordinates,¹⁸ which limits their reliability and physical usefulness. Nevertheless, some authors adopted definitions relying on the concept of polarization density induced by a static electric field to calculate electric dipole polarizabilities^20,21 and hyperpolarizabilities.^22–24

More recently, it has been suggested that origin-invariant densities of electric polarizability, generated by static and dynamic electric fields, e.g., that associated with a monochromatic plane wave, can be defined via the current density vector induced by time-derivative of the electric field,^25,26 oscillating with frequency ω, carried by the wave. In fact, the polarizability density is identical to the current density tensor (CDT)²⁷ and yields a practical recipe for computations, as recently shown.²⁸

Much in the same way, a translationally invariant density of optical rotatory power is obtained by the trace of the CDT associated with the current density vector induced by the magnetic field carried by the plane wave.^26,29 A few compounds, presenting interesting features from the chemical point of view, have been studied here.

This paper is organized as follows: the state of the art is reviewed in Section 2 by recalling the essential theoretical tools previously developed to define current density vectors and tensors needed for a computational approach to origin-invariant densities of electric dipole polarizability (EDP) and natural optical rotation, the latter related to the trace of the mixed electric dipole–magnetic dipole polarizability (MEMDP) tensor; applications are described in two separate parts of section 3, the first dedicated to EDP density, the second to MEMDP density; some final remarks concerning outlook and future work are given in section 4.

2 State of the art

Recent progress in the computation of origin invariant current densities, induced by optical fields associated with a monochromatic wave, is recalled in the following. After a brief preliminary excursus, see section 2.1, outlining the notation adopted, we focus in particular on the connection between induced electronic currents and two response properties, namely static and dynamic electric dipole polarizability, section 2.2, and optical rotatory power of chiral molecules, section 2.3.

Two density functions of static, i.e, ω = 0, and dynamic polarizability tensor α(ω) are examined in section 2.2. The former is usually taken into account^20,21,24 and depends on the origin of the reference system; accordingly it is discarded as computationally unpractical, in favor of an alternative, translationally invariant density , based on the definition of current density tensor linked up with time derivative of the optical electric field carried by monochromatic radiation.

In the following section 2.3, it is shown that the MEMDP tensor κ′(ω), introduced to interpret the optical activity of chiral molecules, can be expressed in terms of a series of density functions²⁶, which can be integrated all over the three-dimensional space to evaluate components and trace of κ′(ω). Computational approaches to the density , based on frequency-dependent electronic current densities induced by monochromatic light shone on a probe molecule, has been developed^26,29 and implemented at the time-dependent density functional theory (TD-DFT) level for the calculation of the optical rotation of chiral molecules.^30–33

The dependence of on the origin of the coordinate system has been investigated in connection with the corresponding change of . It is shown that only the trace of the density function defined via dynamic current density evaluated using continuous translation of the origin of the coordinate system, for static^34–37 and optical magnetic field,²⁵ is invariant of origin. Accordingly, this function is recommended as a tool quite useful to determine the molecular domains which determine optical activity to a major extent.

2.1 Theoretical framework

Adopting the same notation as in previous papers,^28,29,38,39 for a molecule with n electrons and N clamped nuclei, i.e., within the Born–Oppenheimer (BO) approximation,⁴⁰ charge, mass, position, canonical and angular momentum of the k-th electron are indicated, in the configuration space, by −e, m_e, r_k, p̂_k = −iħ∇_k, l̂_k = r_k × p̂_k, k = 1, 2…n. Hereafter, the SI units are used and standard tensor formalism is employed, e.g., the Einstein convention of implicit summation over two repeated Greek indices is in force. Capitals denote n-electron operators, e.g., for position, canonical and angular momentum,
then, electric and magnetic dipole operators become

[small mu, Greek, circumflex] = −eR̂, (1)

Time-dependent quantum mechanical perturbation theory⁴¹ is used to obtain expressions for the polarization charge density and current density induced in the electrons of a molecule by optical fields, assuming that the eigenvalue problem for the time-independent BO electronic Hamiltonian Ĥ⁽⁰⁾Ψ⁽⁰⁾_j = E⁽⁰⁾_jΨ⁽⁰⁾_j has been solved, determining a set of eigenfunctions Ψ⁽⁰⁾_j and corresponding energy eigenvalues E⁽⁰⁾_j. The reference (ground) state is indicated by Ψ⁽⁰⁾_a and the natural transition frequencies are ω_ja = (E⁽⁰⁾_j − E⁽⁰⁾_a)/ħ.

Within the McWeeny normalization,⁴² the n-electron density matrix

depends on electronic space–spin coordinates x_k = r_k⊗s_k, k = 1, 2,…,n, where

X₁ ≡ {x₂,…,x_n}, X = {x₁,X₁}, dX₁ ≡ {dx₂,…,dx_n}. (4)

Inserting the reference (ground) state Ψ⁽⁰⁾_a of the molecule in eqn (3) and integrating over ds₁, one gets

Thus, ρ⁽⁰⁾(r) = −eγ⁽⁰⁾(r) is the electronic charge density in the absence of perturbation.

The probability current density⁴² is defined as

in which one puts r′ = r after operating with the electronic mechanical momentum

[small pi, Greek, circumflex] = p̂ + eA. (7)

The Bloch gauge⁴³ is adopted for the vector potential A. The electron current density corresponding to (6) is obtained by multiplying by −e, i.e., J = −ej. The interaction Hamiltonian considered in the present work does not contain terms depending on electron spin, therefore the probability current density (6) includes only orbital contributions.

2.2 Electric dipole polarizability density

The time-dependent electric field E of a monochromatic plane wave with frequency ω generates a first-order polarization density²⁵ of the electron cloudwhere

To first-order in E, the electric dipole moment induced in the electron distribution is defined by the integral

whereis the electric dipole polarizability within dipole-length formalism, symmetric in the indices α–β, so that α^(R,R)_βα(ω) = α^(R,R)_αβ(ω).

The product of vector components within the integral of eqn (10), a dyadic, can be interpreted as a polarizability density tensor function

The electronic polarization density (8) is independent of origin, since 〈a|[small mu, Greek, circumflex]_β|j〉 is origin-independent, but depends on the origin of r_α. Therefore eqn (10) and (12) are unsuitable for calculations, since they do not meet the fundamental requirement pointed out by Jameson and Buckingham.¹⁸

A more practical computational recipe is arrived at by an alternative procedure.^26,29 To first-order, the time-derivative Ė induces an electronic current density

conveniently expressed by the CDTwhich provides a preferable definition of origin-independent spatial density of dynamic electric polarizability,²⁶sincewhere α^(R,P)_βα(ω) is the electric dipole polarizability in mixed dipole length–dipole velocity formalism. It is symmetric in the indices α–β, and identical to that in the (R,R) formalism (11) only if the off-diagonal hypervirial relationshipis fulfilled. In finite (small) basis set calculations α^(R,P)_βα ≠ α^(R,P)_αβ.

The previous considerations apply in the case of dynamic polarizability; in the static case one simply puts ω = 0 in eqn (14).

Multiplying the origin-independent electric dipole polarizability density tensor by E(ω,t), we obtain

which defines the same induced dipole density vector as (10) if (17) is satisfied. This can also be interpreted, by virtue of eqn (13), as a polarization current density generated by the time-derivative of the oscillating electric field of the radiation.

2.3 Electric dipole–magnetic dipole polarizability density

The first-order perturbed polarization density induced by time-derivative of the dynamic magnetic field is written as^25,26,29where

The first-order electric dipole moment, induced in the electron distribution by time derivative of the time-dependent magnetic field Ḃ of a monochromatic plane wave with frequency ω, is given by^25,26,29

whereis the MEMDP in dipole length-angular momentum (R,L) formalism. This is identical to that in the dipole velocity–angular momentum (P,L) formalismif the off-diagonal hypervirial theorem (17) is satisfied.

The trace of the MEMDP tensor, either (22) or (23), is related to the angle of natural optical rotation through the Rosenfeld equation.^44,45 In SI units,

where μ₀ is the permeability of free space, L is the path length, and is the number density of optically active molecules.⁴⁶ Most commonly used is the specific rotation power^47,48where the molecular parameter . Eqn (25) provides the specific rotation in the usual deg [dm g cm⁻³]⁻¹ units, when the radiation wavelength λ is in cm, the molecular mass M is in g mol⁻¹ and β is in cm⁴. Optical rotations calculated in the dipole-length gauge, i.e., using the trace of MEMDP (22), are origin dependent for approximate calculation methods, whereas in the velocity gauge (23) they are not, even using truncated basis sets. It is usually assumed that the closer optical rotation estimates in length and velocity gauges, the higher the quality of the computation.

A MEMDP density can be defined for the (R,L) formalism via the integrand function in (21)

along with a second definition for the (P,L) formalism, which can be obtained from the electronic current density induced by the oscillating magnetic field²⁷where the CDT, for a fixed origin of the magnetic dipole operator (2), is the sum of two contributions,a paramagnetic term given byand a diamagnetic term determined by the electronic charge density in the absence of perturbation

Indeed, allowing equality

within the paramagnetic term (29), after some algebraic eqn (28) is cast in the formwhereis the MEMDP density in the (P,L) formalism, since spatial integration of (33) leads to eqn (23), i.e.,

MEMDP density (26) depends on the origin of r_α and both densities (26) and (33) depend on the origin of L̂_β. Therefore, their use is impractical for computational purposes. In particular, in view of an optical rotation density, the trace of both densities is also origin dependent and of no practical interest, see ref. 29 for further details.

An origin-invariant definition is arrived at by considering the CDT associated with the frequency-dependent and origin-invariant current density obtained via the CTOCD-DZ computational scheme for static^34–37 and optical magnetic field²⁵

J^B_DZ(r,ω,t) ≡ I^B(r,ω,t) = J^B_p(r,ω,t) + J^B_Δ(r,ω,t). (35)

The corresponding CDT is defined as

In eqn (35) a new term J^B_Δ(r,ω,t) replaces the conventional diamagnetic contribution. The corresponding CDT is given by

Substituting (31) in (37), after some manipulation we obtain the trace

where the MEMDP density for the (L,P) formalism is given by²⁹

It is worth noting that spatial integration of eqn (39) leads to the MEMDP tensor κ′^(L,P) that is simply transposed with respect to κ′^(P,L) given in eqn (23).

Now, eqn (38) suggests that a suitable definition of origin-independent MEMDP density is

whose trace can also be obtained as

In the static limit, i.e., ω → 0, eqn (40) is indeterminate. However, in this limit, the trace (41) vanishes, as can be inferred viaeqn (33) and (39). At any rate, eqn (40) defines an origin-independent MEMDP tensor density for any r all over the molecular domain, regardless of the approximation used for the calculation, which is an essential requisite postulated by Jameson and Buckingham^17,18 for computational purposes. We remark that definition (40) is slightly different from that introduced in ref. 29, since for ω = 0 has been removed. This improved definition finds his analogue for the integrated property in the modified velocity gauge (MVG) formulation of Pedersen et al.⁴⁹ over the origin-independent velocity gauge of Grimme et al.⁵⁰ At any rate, the integral in d³r of the trace of the second (static) term on the r.h.s. of eqn (40) vanishes for a complete basis set. Thus, the trace (41) can be used to obtain an origin-independent specific rotation power density, by means of an equation similar to (25), i.e.,

where β(r) is now a molecular density function equal to . Spatial integration of the specific rotation power density (42) provides exactly the same optical rotation obtained adopting the velocity gauge.

3 Applications

Densities of the electric dipole polarizability (EDPD) and specific rotation power (SRPD) were implemented at the time-dependent Hartree–Fock (TD-HF) and density functional theory (TD-DFT) level within the SYSMOIC program package.^51,52 Full implementation details can be found in ref. 28, 29 and 38, along with a number of preliminary examples, which were presented for a few molecules for both EDPD and SRPD. One of the main features of the implementation consists of transforming transition amplitudes from the time-dependent calculation into perturbed molecular orbitals, which permits orbital decomposition of the results, which are origin-independent for both densities. We believe that this last feature is particularly interesting because it represents a practical analysis tool to deepen the study of the structure–property relationship, especially in the case of the optical rotation, a topic recently addressed for the length gauge.⁵³

In the following, we review the main results obtained so far and present new cases of study, which illustrate the usefulness of EDPD and SRPD in sections 3.1 and 3.2 respectively.

3.1 Electric dipole polarizability density (EDPD)

As discussed in Section 2.2, only the electric dipole polarizability density (15) is quite practical for computational purposes.

So far, we have studied a series of atoms and simple molecules in terms of static and dynamic electric dipole polarizability density²⁸ and in terms of electric dipole polarization current density.³⁸ What emerged can be summarized as follows: (i) despite an always positive polarizability in the transparent region far from absorption, the dipole moment density (18) is found counter-oriented with respect to the induced electric field in regions close to the nuclei of a molecule, with the exception of He and H atoms, as for example, in H₂, LiH, and BeH₂; (ii) in terms of current induced by time derivative of the electric field, these regions of counter-polarization correspond to the inner portion of toroidal currents, which are found ubiquitous in the proximity of all the nuclei in a molecule, except hydrogen, for any orientation of the inducing field; and (iii) quite obviously, the higher the electron density, the higher is the polarizability density.

As a further application, we considered the possibility of describing electron (de)localization in terms of EDPD, focusing in particular on aromatic molecules. As shown previously,⁵⁴ab initio calculations of π current density maps and dipole polarizability contributions confirm the expected correlation for monocycles, but show the limitations of polarizability as a measure of aromaticity – indeed, this is not our aim. On the other hand, it is reasonable to think about a connection between electron density and EDPD and this is what we preferred to look at carefully.

Four molecules have been taken into account: benzene, naphthalene, phenanthrene, and ovalene. Their geometries were optimized at the HF/6-31G(d) level, assuming standard Cartesian axis orientations. Origin independent EDPDs (15) were calculated using the TD-HF method⁵⁵ and the 6-31G(d) basis set.^56,57 For all molecules a radiation wavelength of 633 nm was used, which is sufficiently far away from any absorption. Superimposed streamlines on EDPD plots represent trajectories of the current density induced by Ė, see eqn (13), at t = 0. Terms in-plane and out-of-plane will be used to indicate the diagonal tensor components whose indices correspond to directions parallel and perpendicular to the molecular plane, respectively.

3.1.1 Benzene. As a progenitor of aromatic compounds, benzene is expected to display an EDPD that exemplifies at best the common perception of electronic delocalization people have. For planar aromatics, this is best observed over (or below) the molecular plane, that is, in regions where the π-electron density is higher. In Fig. 1 the average value of EDPD, together with the in-plane components, calculated at points covering a plane parallel to the molecular plane 1 a.u. from it, are displayed as diverging color maps. The out-plane component is negligible, in comparison with the in-plane components, and is not shown. As can be observed, the average value gives a meaningful representation of the electron delocalization, while the in-plane components document strong polarization parallel to the molecular plane. It is interesting to see how the two in-plane components are enhanced in different regions, containing C–C bonds on opposite sides of the molecule. As shown later on, this motif is a common feature of the in-plane components and we will refer to it just as an opposing arrangement.

Fig. 1 Origin-independent electric dipole polarizability density of benzene, calculated at the TD-HF/6-31g(d) level of theory, for a radiation wavelength of 633 nm, on a plane located at 1 a.u. below the molecular plane. Left ; center and right in-plane components. Here and in the following, the EDPD is represented using diverging color maps. Values in a.u. are: 0.6 dark red; 0 grey; −0.6 dark blue.

No regions of negative polarization density can be seen in Fig. 1. On moving to the molecular plane, see Fig. 2, a region of counter-polarization appears around each carbon atom, irrespective of the direction of the inducing electric field. For the two in-plane components, the extension of the negative regions follows the opposing arrangement. Of course, on the molecular plane the π-electron density vanishes, therefore the electric polarization of the σ-electrons can be observed. It is mainly localized around C–H bonds, four of them in one direction and the other two in the other direction, in accordance with the opposing arrangement. In terms of the current density tensor, multiplying by the time derivative of the electric field, see eqn (13), the current density vector induced by an oriented radiation can be obtained. In the benzene molecule, the corresponding vector fields induced by a unitary Ė are superimposed upon the in-plane components in the central and right plots of Fig. 2.

Fig. 2 EDPD of benzene over the molecular plane. Left ; center and right in-plane components. Trajectories (in magenta) of the current density (13) induced by the time derivative of the electric field are superimposed to the in-plane components. See caption to Fig. 1 for other details.

It is worth noting that the section of toroidal flow in each counter-polarization region. In many cases, for the benzene and the other molecules, they have the shape of a topological circumference formed by a pair of heteroclinic streamlines that connect two conjugated saddle-nodes.³⁸ From the saddle nodes, other trajectories emerge, delimiting the vector field. In other cases, the topological circumferences open up and the toroidal flow runs over more than one nucleus. Counter-polarization regions and toroidal currents are strictly related.³⁸

3.1.2 Naphthalene. Discussion of naphthalene EDPD goes much in the same way as before. The EDPD over a plane located 1 a.u. below the molecular plane is shown in Fig. 3. The electron delocalization can be readily observed as a diffuse polarization cloud extending all over the carbon network. Also in this case, the in-plane components follow the opposing arrangement, located on different sets of opposite C–C bonds. Looking carefully at , one observes that the density is not uniformly distributed, in particular it is higher about C_α–C_β bonds. On the molecular plane, see Fig. 4, C–H bonds show the larger polarization in the opposing arrangement fashion. The same occurs in counter-polarization regions close to carbon nuclei. As can be clearly noticed, sections of toroidal flow, represented by topological circumferences of streamlines connecting conjugated saddle-nodes, incorporate perfectly the regions of negative polarization.

Fig. 3 EDPD of naphthalene over a plane located at 1 a.u. below the molecular plane. Left ; center and right in-plane components. See caption to Fig. 1 for other details.

Fig. 4 EDPD in-plane components of naphthalene over the molecular plane. See captions to Fig. 1 and 2 for other details.

3.1.3 Phenanthrene. A commonly accepted description of the phenanthrene molecule consists of two Clar's sextets located on the external rings and a double bond connecting carbons 9 and 10. Such interpretation is in accord with the polarizability density shown in the left plot of Fig. 5. The double bond is clearly visible as a site of major polarization density. However, the density distribution over the two external rings is not perfectly uniform as in benzene. In particular, the density in correspondence of C₂–C₃ (C₆–C₇) is a bit lower with respect to the other bonds. Also in this case, the in-plane components show an opposite arrangement, with the one showing the alternating bonds more intense than the other. Moving on the molecular plane, see. Fig. 6, the opposite arrangement of polarizability density located about C–H bonds is observed. Also the bridge C–C bond linking the two external rings is characterized by strong density. The counter-polarization regions go in parallel, together with the sections of ubiquitous toroidal patterns.³⁸

Fig. 5 EDPD of phenanthrene over a plane located at 1 a.u. below the molecular plane. Left ; center and right in-plane components. See caption to Fig. 1 for other details.

Fig. 6 EDPD in-plane components of phenanthrene over the molecular plane. See captions to Fig. 1 and 2 for further details.

3.1.4 Ovalene. Prompted by the results obtained for phenanthrene, we have considered one more system characterized by having more Clar's sextets and C=C double bonds, i.e., ovalene, in which a maximum of four Clar's sextets and as many C=C double bonds can be allocated. Computed EDPD of ovalene is shown in Fig. 7. Within the plot of , on the left of the figure, one can clearly discern the four C=C double bonds located at the top-right, top-left, bottom-right, and bottom-left corners. Clar's sextets are hidden a little by the non uniform density distribution, see rings at the top, bottom, center-right, and center-left. They are visible with greater clarity in the right plot of Fig. 7 relative to one of the two in-plane components, which are characterized by a slightly more complex, though not less intriguing, opposing arrangement.

Fig. 7 EDPD of ovalene over a plane located at 1 a.u. below the molecular plane. Left ; center and right in-plane components. See caption to Fig. 1 for further details.

Over the molecular plane, ovalene displays the EDPD in-plane components shown in Fig. 8. The opposing arrangement of the polarizability density mainly located around C–H bonds is clearly observable. Some C–C bonds also contribute. The omnipresence of counter-polarization regions, embedded within domains of toroidal currents, is further verified also in ovalene, particularly in the vicinity of the innermost carbon atoms.

Fig. 8 EDPD in-plane components of ovalene over the molecular plane. See captions to Fig. 1 and 2 for further details.

3.2 Specific rotation power density (SRPD)

As discussed above in Section 2.3, only the specific rotation power density (42), connected with the trace (41) of the MEMDP density tensor, and defined by the dynamic origin-invariant current density (35), is suitable for computational purposes.

In order to emphasize the main features of SRPD, we have recently presented some results for the hydrogen peroxide molecule,²⁹ reporting SRPD maps which visualize , at λ = 355 nm for a number of chiral conformations, varying the HO–OH torsional angle θ between 0° and 180°. We observed that SRPD is mainly located in the vicinity of oxygen atoms, with a conspicuous alternation of sign. For the two achiral conformations at θ = 0° and θ = 180°, the symmetry of density is consistent with the obviously vanishing optical rotation. In particular, SRPD changes sign by reflection through a symmetry plane, over which it vanishes at all points. For the intermediate chiral conformations, the absence of symmetry planes gives rise to positive and negative regions which do not cancel one another.

In summary, for this simple example, the SRPD provided by the CDT evaluated via the CTOCD-DZ approach is a function characterized by high intensity peaks of opposite sign in the proximity of atoms, whose symmetry is clearly connected with the integrated property. It provides more precise understanding of how the absence of symmetry planes gives rise to optical rotation.

As a next step, we first focus on the optical activity generated by the presence of a light absorbing unit, such as an inherently achiral carbonyl chromophore, which cannot exhibit any optical rotatory power by itself without the intervention of an extrachromophoric perturbation. The analysis that we can carry out, considering the SRPD, is quite interesting. Passing from the highly symmetrical system constituted by the carbon monoxide molecule to increasingly less symmetrical systems like formaldehyde and acetaldehyde, the latter in two conformations, one achiral and the other chiral, we can study in detail how SRPD changes.

In carbon monoxide, each spatial point of the molecular space belongs to a σ_v symmetry plane of the C_∞v symmetry point group, therefore [Φ]_λ(r) is zero at every r.

In formaldehyde, two σ_v planes of symmetry of the C_2v symmetry point group, divide the SRPD scalar field in four regions of alternating sign. This can be observed on the top panel of Fig. 9 where SRPD changes sign on passing from one side to another of each plane, for a pair of iso-surface values (±100 deg[dm g cm⁻³]⁻¹a₀⁻³), computed by the TD-DFT method.^58,59 The calculation was carried out adopting the B3LYP functional^60,61 and the aug-cc-pVTZ,^62,63 for the radiation wavelength of 589.3 nm (sodium D-line). Obviously, the optical activity vanishes also in this case, i.e., , the four regions equivalent and opposed in sign by symmetry. However, it is interesting to observe that each of the four regions integrates, at this level of theory, to 6895 deg[dm g cm⁻³]⁻¹ in absolute value, that is, a very large estimate, compared with usual specific optical rotations.

Fig. 9 Specific rotation power density, [Φ]_D(r): top formaldehyde; middle acetaldehyde with achiral conformation; bottom acetaldehyde with chiral conformation. Calculated at B3LYP/aug-cc-pVTZ level of theory. Iso-surfaces are ±100 deg [dm g cm⁻³]⁻¹a₀⁻³, red/blue positive/negative.

The achiral conformation of acetaldehyde is characterized by one plane of symmetry. Despite that, as in the case of formaldehyde, the SRPD field around the carbonyl group is partitioned into four regions of alternating sign, which are in absolute value equivalent by symmetry only by reflection through the unique symmetry plane, see plot in the middle of Fig. 9. Consequently, there is no optical activity, since the positive contribution to the specific optical rotation (+8516 deg [dm g cm⁻³]⁻¹) cancels out exactly the negative one (−8516 deg [dm g cm⁻³]⁻¹).

A chiral conformation of acetaldehyde was set by rotating the methyl group until one of the protons is at right angles to the carbonyl plane. The plot of calculated SRPD is displayed in the bottom panel of Fig. 9. Remarkably, the loss of symmetry does not change what appears to be an essential feature of the carbonyl group represented by the separation of the SRPD in quadrants, which are no longer equivalent in absolute value by symmetry.

The decomposition into origin-independent orbital contributions provides an explanation of this separation, which recalls some aspects of the famous octant rule.⁶⁴ Indeed, according to the results of the TD-DFT calculation, the first excited state is a virtual transition from the HOMO to the LUMO (n → π*). Moreover, the orbital decomposition analysis shows that the HOMO provides the dominant contribution to SRPD, since plots calculated using only the HOMO contribution are indistinguishable from that calculated ones using the full set of MOs. Now, the HOMO shape is shown in Fig. 10 and, as can be observed, it presents a nodal surface which bisects the H₃ C–C–H angle and is perpendicular to the plane of the carbonyl group. Such a nodal surface provides, together with the nodal surface of the π* LUMO, the observed separation in quadrants of SRPD. A similar discussion based on transition densities has also been reported in ref. 32.

Fig. 10 The HOMO of acetaldehyde, B3LYP/aug-cc-pVTZ.

Looking at the two signed surfaces in the bottom panel of Fig. 9, it is not easy to notice a predominance of one over the other. However, spatial integration of the SRPD leads to a non vanishing negative specific rotatory power, [ϕ]_589.3 = −69.2 deg [dm g cm⁻³]⁻¹, resulting from a negative contribution which is only a little bigger than the positive contribution, i.e., −8314.7 and +8245.5 deg [dm g cm⁻³]⁻¹ respectively.

The remarkable ratio of approximately two orders of magnitude between individual contributions and total specific rotatory power is commonly found also in other systems. Therefore, it may well account for the known difficulties to obtain accurate theoretical optical rotations,^30,65 as even small effects, for example, interaction with the solvent^66,67 or molecular vibration,^68–70 can produce considerable variations of the final result, which also may easily change sign.

However, to check the quality of the results exposed above, we report in Table 1 the specific rotatory power of the chiral conformation of acetaldehyde, calculated for three commonly used wavelengths, using basis sets of increasing quality at the B3LYP level of theory; GIAO predictions, along with optical rotations determined within modified-velocity and length gauges, are listed. Considering the GIAO results obtained adopting the aug-cc-pV5Z basis set as the best of the series, it is interesting to observe the good performance of the modified-velocity and length gauges already using the TZ basis set.

Table 1 Specific rotatory power (deg [dm g cm⁻³]⁻¹) of a chiral conformation of acetaldehyde calculated at the B3LYP level for three commonly used radiation frequencies adopting aug-cc-pVXZ basis sets of increasing accuracy

λ (nm)	Basis set	GIAO	Modified-velocity	Length
633	DZ	−52.7	−49.7	−54.8
	TZ	−53.6	−57.4	−56.3
	QZ	−55.2	−56.4	−56.3
	5Z	−56.1	—	—
589.3	DZ	−63.7	−59.7	−66.1
	TZ	−64.8	−69.2	−68.0
	QZ	−66.6	−68.1	−68.0
	5Z	−67.7	—	—
355	DZ	−414.7	−349.1	−417.5
	TZ	−412.0	−439.1	−433.8
	QZ	−428.4	−439.1	−439.8
	5Z	−436.8	—	—

---

3.2.1 (1S,4S)-Norbornenone. Of course, due to rapid rotation of the methyl group, none of the chiral conformations of acetaldehyde is accessible. Therefore, we have extended our study to a real chiral system, in which the achiral C=O chromophore determines its optical rotation. Among many possible choices, we have considered (1S,4S)-norbornenone as a suitable candidate for its large optical rotation and for a number of other interesting facts briefly summarized as follows. Owing to substantially similar solution phase optical rotations measured at the sodium D-line for different solvents,^71–73 it was assumed for a long time that the optical rotation of this molecule is very little affected by solvent effects. The presumed lack of solvent effect was somehow hidden also by the apparently good agreement found with the B3LYP/aug-cc-pVDZ estimate of the OR obtained for the free molecule.³⁰ However, as reported successively,⁷⁴ this was found to be not true. Notably, the experimental rotatory power exhibited by isolated (vapor-phase) molecules is much lower in magnitude and in good accord with coupled cluster predictions computed using both the origin-independent dipole-velocity and origin-dependent dipole length representations.^74,75 Last but not least, the constructive interaction of the two chromophores (C=O and C=C) has been invoked to explain the otherwise unusually large specific rotation of (1S,4S)-norbornenone.⁷⁶

In the light of what has been discussed above for the simple acetaldehyde case, we have now the possibility to analyze where the large specific rotation of (1S,4S)-norbornenone originates in terms of SRPD and, in particular, where major contributions are located and how many they are.

For the calculation of the MEMDP we have used the B3LYP/aug-cc-pVDZ combination of functional and basis set as adopted by Stephens et al. in ref. 30. The result is shown in Fig. 11. As can be observed, there is an amazing partition of the SRPD in quadrants, quite similar to that studied before for much simpler systems. This makes us think that we are really observing the essential features of MEMDP of the carbonyl chromophore. Moreover, as can be readily noted in Fig. 11, the predominance of blue sectors over the red ones well agrees with the negative value of [ϕ]_D, estimated via spatial integration as large as −1195 deg [dm g cm⁻³]⁻¹, in accord with previous DFT calculations.^30,74,75 Once again, we observe very large positive and negative contributions, which are calculated to be 4704 and −5899 deg [dm g cm⁻³]⁻¹. This explains why even small variations of big contributions lead to large optical rotations, which turn out to be very sensitive to the environment, i.e. free molecule vs. solution phase, and the calculation method. Norbornenone contains a second chromophore, i.e., a C=C double bond, which, seemingly, does not play a substantial role in determining the total optical rotation density at the sodium D-line, apart from the possibility of providing an influence as an extrachromophoric perturbing agent.⁷⁶

Fig. 11 Specific rotation power density of (1S,4S)-norbornenone calculated at the B3LYP/aug-cc-pVDZ level of theory for λ = 589.3 nm. Iso-surfaces are ±100 deg [dm g cm⁻³]⁻¹a₀⁻³, red/blue positive/negative.

4 Conclusions

Two origin-independent polarizability density tensor functions have been defined and evaluated for a number of molecules. They are as follows: (i) the time-dependent electric dipole polarizability density (EDPD) and (ii) the average time-dependent mixed electric dipole–magnetic dipole polarizability density (MEMDPD). They have been employed to rationalize origin-independent properties induced by optical fields, thus paving the way for topological studies aimed at discovering domains of molecular space and functional groups which influence electromagnetic response to a higher extent.

The origin-independence of these densities constitutes an essential requisite of their physical meaning, quite recommendable for computational purposes.

EDPD coincides with a current density tensor fully characterizing the electron flow induced by time derivative of the oscillating electric field of the radiation. Its integral equals the dynamic electric dipole polarizability, therefore EDPD is effectively an electric dipole polarizability density. In the limit of zero frequency, it defines the static polarizability density.

EDPD is shown here to be useful for detecting regions of electron distribution in which polarization density is larger and where counter-polarization densities appear, in connection with a peculiar shape of induced currents, referred to as toroidal flow. Despite the always positive value of total static electric dipole polarizability far from absorption, counter-polarization in regions close to the atomic nuclei of a molecule, except protons, is ubiquitous and always deeply connected with toroidal electron currents.

Application to aromatic hydrocarbons confirms the ability of EDPD to indicate where electron delocalization occurs and to distinguish aromatic rings and localized double bonds.

Much remains to be done. In perspective, we would like to put in a to-do list: (i) the need for more clearly differentiating various kinds of aromatic rings, which was partially achieved in the molecules examined, but still deserves to be investigated at the quantitative level and (ii) extending the topological study of polarizability density to different kinds of functional groups and verifying if their features are transferable from one molecule to another.

Also MEMDPD is itself a current density tensor connected with the electronic currents induced by the oscillating magnetic field of the radiation. It vanishes in the limit of zero frequency. Its trace is related to the optical rotation density, or, more conveniently, to the specific rotation power density (SRPD).

As far as SRPD is concerned, we remark that it can be dissected in origin-independent orbital contributions and that it integrates to obtain MVG optical rotations. It can be used in conjunction with a variety of quantum mechanical methods, coupled cluster in particular,⁴⁹ to analyze the optical rotation of chiral molecules without recourse to London orbitals. The perspective is to extend the study to other achiral chromophores, e.g., double and triple C–C bonds, cyano groups, aromatic rings, etc., to cite quite a few. Moreover, a very interesting case worth studying would be that of chromophores that are themselves chiral, as, for instance, helicenes and compounds containing non-planar diene moieties. Last but not least, a very important goal to achieve would be that of predicting in which way a chiral perturbation acts on the SRPD topology.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Financial support from MUR (Ministero dell'Università e della Ricerca), FARB 2020 and FARB 2021, is gratefully acknowledged.

References

1. R. F. W. Bader, in Int. Rev. Sci. – Theoretical chemistry, ed. A. D. Buckingham and C. A. Coulson, Butterworths, London, Boston, 1975, p. 43.
2. R. F. W. Bader, Atoms in molecules: a quantum theory, Clarendon Press, Oxford University Press, Oxford [England], New York, 1994.
3. P. Popelier, P. Smith and A. Hinchliffe, Quantum topological Atoms, Royal Society of Chemistry, United Kingdom, 2002, vol. 8, pp. 391–448.
4. A. Morgenstern, T. R. Wilson and M. E. Eberhart, J. Phys. Chem. A, 2017, 121, 4341–4351.
5. A. D. Buckingham, Advances in Chemical Physics, John Wiley & Sons, Inc., 2007, pp. 107–142.
6. O. Theimer and R. Paul, J. Chem. Phys., 1965, 42, 2508–2517.
7. W. H. Orttung and D. Vosooghi, J. Phys. Chem., 1983, 87, 1432–1437.
8. D. W. Oxtoby and W. M. Gelbart, Mol. Phys., 1975, 30, 535–547.
9. J. Sipe and J. Van Kranendonk, Mol. Phys., 1978, 35, 1579–1584.
10. W. Maaskant and L. Oosterhoff, Mol. Phys., 1964, 8, 319–344.
11. K. L. C. Hunt, J. Chem. Phys., 1983, 78, 6149–6155.
12. K. L. C. Hunt and J. E. Bohr, J. Chem. Phys., 1986, 84, 6141–6150.
13. K. L. C. Hunt, J. Chem. Phys., 1984, 80, 393–407.
14. K. L. C. Hunt, J. Chem. Phys., 1989, 90, 4909–4915.
15. K. L. C. Hunt and Y. Q. Liang, J. Chem. Phys., 1991, 95, 2549–2559.
16. K. L. C. Hunt and R. A. Harris, J. Chem. Phys., 1991, 94, 6995–7002.
17. C. J. Jameson and A. D. Buckingham, J. Phys. Chem., 1979, 83, 3366–3371.
18. C. J. Jameson and A. D. Buckingham, J. Chem. Phys., 1980, 73, 5684.
19. J. D. Jackson, Classical electrodynamics, Wiley, New York, 3rd edn, 1999.
20. A. Alparone, Sci. World J., 2013, 2013, 1–9.
21. N. Otero, C. Van Alsenoy, C. Pouchan and P. Karamanis, J. Comput. Chem., 2015, 36, 1831–1843.
22. M. Nakano, I. Shigemoto, S. Yamada and K. Yamaguchi, J. Chem. Phys., 1995, 103, 4175–4191.
23. S. Yamada, M. Nakano, I. Shigemoto, S. Kiribayashi and K. Yamaguchi, Chem. Phys. Lett., 1997, 267, 445–451.
24. M. Nakano, H. Fukui, T. Minami, K. Yoneda, Y. Shigeta, R. Kishi, B. Champagne, E. Botek, T. Kubo, K. Ohta and K. Kamada, Theor. Chem. Acc., 2011, 130, 711–724.
25. P. Lazzeretti, J. Chem. Phys., 2018, 149, 154106.
26. P. Lazzeretti, J. Chem. Phys., 2019, 150, 184117.
27. P. Lazzeretti, J. Chem. Phys., 2018, 148, 134109.
28. F. F. Summa, G. Monaco, P. Lazzeretti and R. Zanasi, J. Phys. Chem. Lett., 2021, 12, 8855–8864.
29. F. F. Summa, G. Monaco, R. Zanasi, S. Pelloni and P. Lazzeretti, Molecules, 2021, 26, 4195.
30. P. J. Stephens, F. J. Devlin, J. R. Cheeseman and M. J. Frisch, J. Phys. Chem. A, 2001, 105, 5356–5371.
31. J. Autschbach, S. Patchkovskii, T. Ziegler, S. J. A. van Gisbergen and E. J. Baerends, J. Chem. Phys., 2002, 117, 581–592.
32. J. Autschbach, Chirality, 2009, 21, E116–E152.
33. J. Autschbach, L. Nitsch-Velasquez and M. Rudolph, in Time-Dependent Density Functional Response Theory for Electronic Chiroptical Properties of Chiral Molecules, ed. R. Naaman, D. N. Beratan and D. Waldeck, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011, pp. 1–98.
34. P. Lazzeretti, M. Malagoli and R. Zanasi, Chem. Phys. Lett., 1994, 220, 299–304.
35. P. Lazzeretti, Theor. Chem. Acc., 2012, 131, 1222.
36. R. Zanasi, P. Lazzeretti, M. Malagoli and F. Piccinini, J. Chem. Phys., 1995, 102, 7150.
37. R. Zanasi, J. Chem. Phys., 1996, 105, 1460.
38. F. F. Summa, G. Monaco, R. Zanasi and P. Lazzeretti, J. Chem. Phys., 2022, 156, 054106.
39. F. F. Summa, G. Monaco, R. Zanasi and P. Lazzeretti, J. Chem. Phys., 2022, 156, 154105.
40. M. Born and R. Oppenheimer, Ann. Phys., 1927, 389, 457–484.
41. P. W. Langhoff, S. T. Epstein and M. Karplus, Rev. Mod. Phys., 1972, 44, 602–644.
42. R. McWeeny, Methods of molecular quantum mechanics, Academic Press, London, 2nd edn, 1992.
43. F. Bloch, in W. Heisenberg und die Physik Unserer Zeit, ed. F. Bopp, Friedr. Wieveg & Son, Braunschweig, 1961, pp. 93–102.
44. L. Rosenfeld, Z. Phys., 1929, 52, 161–174.
45. L. Rosenfeld, Theory of Electrons, North-Holland Publishing Company, 1951.
46. L. D. Barron and A. D. Buckingham, Annu. Rev. Phys. Chem., 1975, 26, 381–396.
47. E. U. Condon, Rev. Mod. Phys., 1937, 9, 432–457.
48. P. Stephens, F. Devlin, J. Cheeseman and M. Frisch, Chirality, 2002, 14, 288–296.
49. T. B. Pedersen, H. Koch, L. Boman and A. M. S. de Merás, Chem. Phys. Lett., 2004, 393, 319–326.
50. S. Grimme, F. Furche and R. Ahlrichs, Chem. Phys. Lett., 2002, 361, 321–328.
51. G. Monaco, F. F. Summa and R. Zanasi, J. Chem. Inf. Model., 2020, 61, 270–283.
52. G. Monaco, F. F. Summa and R. Zanasi, SYSMOIC Version 3.3, 2022, https://SYSMOIC.chem.unisa.it/MANUAL/.
53. T. Balduf and M. Caricato, Theor. Chem. Acc., 2023, 142, 11.
54. P. Fowler and A. Soncini, Chem. Phys. Lett., 2004, 383, 507–511.
55. P. Jørgensen and J. Linderberg, Int. J. Quantum Chem., 1970, 4, 587–602.
56. W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 56, 2257–2261.
57. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
58. M. E. Casida, Recent Advances in Computational Chemistry, WORLD SCIENTIFIC, 1995, vol. 1, pp. 155–192.
59. M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem., 2012, 63, 287–323.
60. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
61. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
62. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007.
63. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
64. W. Moffitt, R. B. Woodward, A. Moscowitz, W. Klyne and C. Djerassi, J. Am. Chem. Soc., 1961, 83, 4013–4018.
65. P. Stephens, D. McCann, J. Cheeseman and M. Frisch, Chirality, 2005, 17, S52–S64.
66. Y. Kumata, J. Furukawa and T. Fueno, Bull. Chem. Soc. Jpn., 1970, 43, 3920–3921.
67. B. Mennucci, J. Tomasi, R. Cammi, J. R. Cheeseman, M. J. Frisch, F. J. Devlin, S. Gabriel and P. J. Stephens, J. Phys. Chem. A, 2002, 106, 6102–6113.
68. J. Kongsted, T. B. Pedersen, M. Strange, A. Osted, A. E. Hansen, K. V. Mikkelsen, F. Pawlowski, P. Jørgensen and C. Hättig, Chem. Phys. Lett., 2005, 401, 385–392.
69. K. Ruud and R. Zanasi, Angew. Chem., Int. Ed., 2005, 44, 3594–3596.
70. J. Kongsted, T. B. Pedersen, L. Jensen, A. E. Hansen and K. V. Mikkelsen, J. Am. Chem. Soc., 2005, 128, 976–982.
71. K. Mislow and J. G. Berger, J. Am. Chem. Soc., 1962, 84, 1956–1961.
72. J. M. Janusz, L. J. Gardiner and J. A. Berson, J. Am. Chem. Soc., 1977, 99, 8509–8510.
73. D. A. Lightner, J. K. Gawronski and T. D. Bouman, J. Am. Chem. Soc., 1980, 102, 5749–5754.
74. P. Lahiri, K. B. Wiberg, P. H. Vaccaro, M. Caricato and T. D. Crawford, Angew. Chem., Int. Ed., 2013, 53, 1386–1389.
75. T. D. Crawford and P. J. Stephens, J. Phys. Chem. A, 2008, 112, 1339–1345.
76. M. Caricato, P. H. Vaccaro, T. D. Crawford, K. B. Wiberg and P. Lahiri, J. Phys. Chem. A, 2014, 118, 4863–4871.

---

Footnote

† All authors contributed equally to this work.

---