Electronic Currents and Anapolar Response Induced in Molecules by Monochromatic Light

: It is shown that the electric dipole- and electric quadrupole–anapole polarizabilities, denoted respectively by f (cid:48) αβ and g (cid:48) α , βγ , and the anapole magnetizability a αβ , are intrinsic properties of the electron cloud of molecules responding to optical ﬁelds. f (cid:48) αβ is a nonvanishing property of chiral and achiral compounds, whereas a αβ is suitable for enantiomer discrimination of chiral species. They can conveniently be evaluated by numerical integration, employing a formulation complementary to that provided by perturbation theory and relying on the preliminary computation of electronic current density tensors all over the molecular domain. The origin dependence of the dynamic anapolar response is rationalized via related computational techniques employing numerical integration, as well as deﬁnitions of molecular property tensors, for example, electric dipole and electric quadrupole polarizabilties and magnetizability. A preliminary application of the theory is reported for the R a enantiomer of the hydrogen peroxide molecule, evaluating tensor components of electric dipole-anapole polarizability and anapole magnetizability as functions of the dihedral angle φ ≡ ∠ H-O-O-H in the range 0 ◦ ≤ φ ≤ 180 ◦ .


Introduction
A time-independent magnetic field induces toroidal electron flow in chiral molecules [1][2][3][4][5][6][7][8], giving rise to a static toroidal moment, which can be represented by a vector odd under charge conjugation C, parity P and time reversal T. The optical fields associated with a beam of monochromatic light, oscillating with frequency ω, also induce the anapolar response, which can be rationalized in terms of dynamic anapole polarizabilities and magnetizabilities [9].
The relevance of the anapolar response cannot be overemphasized after the advent of metamaterials, containing meta-atoms endowed with strong toroidal moments, used in a wide series of promising technological applications [10].
Within the electric quadrupole approximation, two second-and one third-rank tensors, depending on ω, have been introduced [9], namely anapole-electric dipole polarizability f αβ , anapole magnetizability a αβ , and anapole-electric quadrupole polarizability g α,βγ . Corresponding quantum-mechanical definitions have been reported within the framework of time-dependent perturbation theory [11], and their magnitude, as well as fundamental properties under C, P, and T symmetry operations [12], have been discussed. In particular, the origin dependence of these tensors has been investigated, determining equations that connect values corresponding to a passive translation of the coordinate system [9].
The present paper sets out to propose a derivation of f αβ , a αβ , and g α,βγ alternative to that arrived at via perturbation theory, allowing for optional definitions of quantummechanical anapole moments in terms of electronic current densities induced by monochromatic light shining on a molecule. Instead of going along an earlier trodden path [9], our investigation aims to reconceptualize a few aspects of the anapolar response. In fact, this novel approach is expected to throw light on fundamental physical mechanisms that come into play in determining the anapolar feedback, as it preludes to future visualization of current density vector fields, toroidal flow in particular.
The paper is organized as follows. An outline of notation is presented in Section 2 and a short description of the theoretical methods employed to compute anapolar response properties is provided in Section 3, introducing current density tensors which depend on the frequency ω of a beam of monochromatic light radiating on a probe molecule. An application of theoretical methods relying on induced current densities is discussed in Section 4 and concluding remarks are made in Section 5.

Outline of Notation and Theoretical Methods
Within the Born-Oppenheimer (BO) approximation [13], for a molecule with n electrons and N clamped nuclei, charge, mass, position and canonical momentum of the k-th electron are indicated, in the configuration space, by −e, m e , r k ,p k , k = 1, 2, . . . , n, using boldface letters for electronic vector operators. Analogous quantities for nucleus I are Z I e, M I , R I , and so forth, for I = 1, 2, . . . , N.
The imaginary unit is represented by a Roman i. Throughout this paper, SI units are used and standard tensor formalism is employed; for example, the Einstein convention of implicit summation over two repeated Greek indices is applied. The third-rank Levi-Civita pseudotensor is indicated by αβγ . Capitals denote n-electron vector operators; for example, for position, canonical and angular momentum, except for the electric and magnetic dipole, Electronic tensor operators are expressed, specifying the components; for example, for the second-rank electric quadrupole operator for the diamagnetic contribution to magnetizability, for the magnetic quadrupole operator in Hermitian form and for the paramagnetic contribution to the anapole operator A "perturbed" operator, linearly dependent on an external magnetic field B, is defined via the relation,â where the diamagnetic contribution is given bŷ Expressions for the polarization charge density and current density induced in the electrons of a molecule by a monochromatic plane wave are obtained by time-dependent quantum mechanical perturbation theory [11], assuming that the eigenvalue problem for the time-independent BO electronic HamiltonianĤ (0) Ψ j has been solved to determine a set of eigenfunctions Ψ A general definition of n-electron density matrices is used throughout this work allowing for the McWeeny normalization [14], They are functions of electronic space-spin coordinates x k = r k ⊗ s k , k = 1, 2, . . . , n, where thus, integrating over ds 1 , one gets from Equation (9), for the reference (ground) state Ψ a of the molecule. The probability current density [14] is obtained from Equations (9)-(11) for the density matrix, in which one puts r = r after operating with the electronic mechanical momentum adopting the Bloch gauge [15] for the vector potential A. The corresponding charge current density is obtained multiplying by −e, that is, J = −ej. The interaction Hamiltonian considered in the present work does not contain terms depending on electron spin, therefore the probability current density Equation (12) includes only orbital contributions.

Current Density Approach to Anapolar Response Properties and Their Origin Dependence
The time derivative of the electric field associated with a beam of monochromatic light impinging on a molecule induces an electronic current density JĖ(r) [16]. A corresponding current density tensor (CDT), to first order in the electric field, is obtained by differentiation, Allowing for the long-wavelength approximation, only electric dipole terms are retained within the Goeppert-Mayer interaction Hamiltonian [17,18]. The electric field of the monochromatic wave is assumed spatially uniform, and Equation (14) defines a quantity independent of the origin of the coordinate system.
The mixed anapole-electric dipole polarizability induced by the time derivative of the electric field of monochromatic radiation can be cast in the form [9] The interesting feature of Equation (16), and of analogous relationships (34) and (41), is that they provide computational recipes alternative to the definitions determined by Rayleigh-Schrödinger perturbation theory [9] by expressing them as integrals of current density tensors. Thus anapole polarizabilities and magnetizability can be calculated by numerical integration using standard methods [19,20]. In this work, spatial integration of the density functions has been performed using the Becke scheme, Ref. [21] adopting 131 angular points for the Lebedev's quadrature of 59th order of accuracy [22] and 131 radial points for the Gauss-Chebyshev radial quadrature of the second-kind [23].
According to Equation (16), the behavior of the f αβ polarizability under the fundamental symmetries is different from that of the anapole magnetizability a αβ [5]. Whereas the former is even under C, P and T, the latter is Cand T-even and P-odd, therefore only a αβ can be used for chiral discrimination [9].
Since the CDT in Equation (14) is invariant under a translation of coordinates, the origin dependence of f αβ is only due to the prefactor within brackets in Equation (16), which will change in a passive transformation of reference system, corresponding to an arbitrary shift d of the origin. In fact, a simple calculation allowing for the second identity of (16) shows that two different observers, at the origin of refer-ence systems K and K , would calculate, or measure, different anapole-electric dipole polarizabilities related by The fourth term on the r.h.s. of Equation (18) is soon evaluated via Equation (90) of Ref. [24], as the argument (14) of the integral defines an origin-independent density of dynamic electric dipole polarizability in the (R, P) gauge, identical to that in the (R, R) and (P, P) gauges if the hypervirial theorem [25][26][27] a is satisfied [28]. Calculation of the other integrals appearing in Equation (18) is greatly facilitated by employing a general and quite effective relationship for the moments of total electronic current density [24]. According to Equation (82) of that paper, to first-order in the perturbing fields of the monochromatic wave, Within the algebraic approximation [29], right-and left-hand sides of Equation (23) are expected to yield different numerical results, unless the basis set is close to completeness, or flexible enough to provide accurate representation ofR,P andL at the same time, and to satisfy hypervirial theorems [25][26][27] to a good extent. In fact, the current density tensors in Equations (14), (26) and (27) are expressed in terms ofp, whereas ∆ μ αβ is usually evaluated within the dipole length gauge in Equation (28). Thus the degree to which the equality (23) is fulfilled provides a measure of basis set quality.
This relationship is applied in Equation (18) by taking into account terms that depend on the same perturbing fields at both sides. Therefore, within the electric quadrupole approximation [16], it is properly used relaxing the long wave-length assumption, that is, supposing that JĖ β (r, ω) is not origin independent, since it varies according to the change of electric field at different points within the molecular region, At any rate, in Equations (23) and (24), both the electric field gradient ∇ γ E β ≡ E γβ and the magnetic field of the impinging radiation are assumed to be spatially uniform all over the molecular domain. They are coupled with the electric quadrupole and magnetic dipole operators, respectively, within the interaction Hamiltonian. This amounts to considering also the current density vector fields induced by the gradient ∇Ė and the magnetic field B, that is, J in agreement with the electric quadrupole approximation. The current density terms appearing in Equation (25) are taken into account introducing the corresponding CDTs [16], and On the r.h.s. of Equation (23), the electric quadrupole induced in the electronic cloud is expressed, within the electric quadrupole approximation, in terms of intrinsic tensor properties; see Appendix A, An analogous expression is used for the induced magnetic dipole: From Equations (18) and (23), one finds hence the relationship connecting the anapole polarizability f αβ for the different origins becomes the same as that arrived at via a more complicated procedure [9], allowing for the translation of operators (6) and (7). Thus the origin dependence of f αβ is fully determined by the prefactor in round brackets in Equation (16).
The mixed anapole-electric quadrupole polarizability induced by the time derivative of the gradient E γβ associated with the monochromatic beam [9] can be expressed in the form: even under C and T and odd under P, hence suitable for chiral recognition. A procedure allowing for the second identity in Equation (34), analogous to that leading to (18), gives The integral in the fourth term on the r.h.s. of Equation (35) corresponds to the electric dipole-electric quadrupole polarizability in the dipole velocity gauge [24]: identical to that in the dipole length gauge, reported in the Appendix A, Equation (A1), if the hypervirial theorem (22) is fulfilled. From Equations (23), (28) and (29) one finds (r − r )JĖ βγ α d 3 r = α βγ, α (r ) − δα D δ,βγ (r ), (37) (r δ − r δ )JĖ βγ δ d 3 r = α βγ,δδ (r ), (r α − r α )JĖ βγ δ d 3 r = α βγ,αδ (r ) − α δ D δ,βγ (r ); then the relationship connecting the anapole-electric quadrupole polarizability g α,βγ for the different origins becomes: Thus the prefactor within round brackets in Equation (34) yields only a few terms linear and quadratic in the d shift, as shown by comparison with Equation (96) of Ref. [9]. Missing contributions would be recovered by coupling the prefactor to terms arising from a change of origin of the JĖ βγ α CDT. In fact, only the total first order current density (25) is origin independent, whereas the three contributions on the r.h.s. interchange in a passive transformation of the coordinate system [16].
The anapole magnetizability is evaluated via a magnetic-field induced current density tensor (27), and a straightforward calculation using the first identity in Equation (41) for different origins yields The fifth term on the r.h.s. is expressed via the magnetic CDT (27), which is related to a density of electric dipole-magnetic dipole polarizability [24], and the other integrals are calculated allowing for Equations (23), (28) and (29). Thus one finds then the relationship connecting the anapole magnetizability a γδ for the different origins becomes Therefore, for the anapole magnetizability (41), only a few terms of those appearing in Equation (93) of Ref. [9] are recovered, for the reasons explained in the discussion following Equation (35).
For static anapole magnetizabilties, Equation (44) becomes [5] a γδ r = a γδ r + 1 2 αβγ ξ δα d β , thus the trace a γγ is invariant of the origin for ω = 0. For a freely tumbling chiral molecule in disordered media, the isotropic anapole moment induced in the direction of the static magnetic field is ∆ â α = (1/3)a γγ B α , whereas the isotropic magnetic dipole moment induced by the curl of a static magnetic field B, that is,

Results and Discussion
A preliminary test of the theoretical scheme outlined in Sections 2 and 3 has been conducted on the hydrogen peroxide molecule, by computing the anapole-electric dipole polarizability f αβ , Equation (16), and anapole magnetizability a αβ , Equation (41). The origin dependence of the former has been investigated via Equations (18) and (33), which account for all the terms reported in a previous study [9].
Calculations have been carried out for the R a enantiomer of H 2 O 2 shown in Figure 1, evaluating both f αβ and a αβ tensor components as functions of the dihedral angle φ = ∠ H-O-O-H, at the time dependent Hartree-Fock (TD-HF) level of theory, equivalent to random-phase approximation (RPA), adopting a fairly large basis set, that is, the uncontracted d-aug-cc-pVQZ [30][31][32] on hydrogen atoms and d-aug-cc-pVTZ [30][31][32] on oxygen atoms. Both basis sets have been downloaded from BSE [33,34]. In our implementation, outlined in detail in a previous reference [35], transition amplitudes S j and T j and corresponding transition energies have been obtained by means of TD-HF calculations. The complete procedure for computing frequency dependent anapole polarizabilities f αβ and anapole magnetizabilities a αβ has been coded within the freely available SYSMOIC program package [36]. The results obtained via numerical integration [21] of Equations (16)       Quite remarkably, the curves representing isotropic anapole magnetizability in Figures 4 and 5, for the frequencies λ = 355 nm and λ = 589.3 nm, vanish in the proximity of a dihedral angle φ ≈ 100 • . The presence of an accidental null value of a αα for a geometrically chiral conformation of hydrogen peroxide, also observed in previous studies of its static anapole magnetizability [37], is due to a different sign of diagonal components: a xx is positive and quite large for φ ≈ 100 • , offsetting the algebraic sum of a yy and a zz .
A test of accuracy of computations is provided by results obtained via numerical integration, reported in Table 1. The quality of the basis sets adopted, and of its ability to give quite reliable representations of operatorsR,P andL, is evident by comparison of electric dipole polarizabilities α αβ in different gauges (R, R), (P, R) and (P, P), as well as electric dipole-magnetic dipole polarizabilities κ αβ in (R, L) and (P, L) formalisms: computed tensor components are nearly coincident to three or more significant figures. A similar conclusion is arrived at by inspection of Table 2, reporting electric dipoleelectric quadrupole polarizability within dipole length (R) and dipole velocity (P) formalisms, Equations (36) and (A1). A further proof of the near Hartree-Fock accuracy of the present calculations is obtained by comparing results from either side of Equation (23), fulfilling the off-diagonal hypervirial relations (22) A broad discussion about the off-diagonal hypervirial relationship, concerning their fulfillment within incomplete basis sets, can be found in Ref. [28]. Quite small differences among (R, R), (R, P) and (P, P) gauges should only be imputed to basis set quality and not to the numerical integration method.

Concluding Remarks and Outlook
The numerical results, arrived at in the present study within the theoretical scheme outlined in Section 3 and the discussion in Section 4, show that frequency dependent anapolar properties of chiral and achiral molecules responding to monochromatic radiation can be easily computed by simple relationships expressed in terms of induced electronic current densities.
Analogous procedures were found to be applicable to study the origin dependence of anapole polarizabilities and magnetizability.
The method proposed and implemented in the present work is an alternative to the conventional formulation based on perturbation theory. It appears to be practical to apply and is computationally more advantageous, since it avails itself of numerical integration, and is quite easy to perform whenever three-dimensional grids of values are available for the current density vector fields induced by the magnetic field, the time derivative of the electric field and of the electric field gradient associated with the impinging radiation. Visualizations of these current densities, which will be available in future investigations, are expected to yield very useful tools for the rationalization of mechanisms coming into play and giving rise to the anapolar response of electrons in a molecule. Effects due to electron correlation will also be investigated in future work.
Funding: Financial support from MIUR is gratefully acknowledged.

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