On the computation of frequency‐dependent molecular magnetizabilities via dynamical charge and current electron densities

In the computation of molecular dynamic magnetizabilities and magnetic dipole moments, three different reference points are required: (i) origin of coordinate system, (ii) origin of vector potential A , and (iii) origin of multipole expansion. This study shows that methods relying on continuous translation of origin of the current density IBrωt induced by optical magnetic fields provide an effective solution to the problem of choices (i) and (ii), in that they yield origin independent IB within the algebraic approximation, for any basis set. Frequency‐dependent magnetizabilities are also invariant with respect to (iii), as a consequence of symmetry, for a number of molecular point groups. In molecules of lower symmetries, computed magnetizabilities depend on origin of the multipole expansion. Large basis set computations carried out for water, ammonia, methane, ethane, ethylene, boranylborane, and hydroxilamine, at the DFT level, have been reported to document these statements. A comparison is made for results obtained within the conventional common origin approach for static magnetic field. Sum rules for invariance of computed properties are discussed. Representations of streamlines and stagnation graphs of dynamical current density vector field induced in the water molecule by monochromatic waves of four frequencies are displayed.


| INTRODUCTION
It is usually argued that only the first nonvanishing multipole moment is independent of the choice of origin, 1,2 that is, an n-pole is origininvariant if the n À 1-pole equals zero. For instance, the electric dipole moment of an uncharged molecule does not depend on arbitrary translations of the reference point defining probe position with respect to a coordinate system. Analogously, the electric quadrupole moment of carbon dioxide is origin independent because its electric dipole vanishes. On the contrary, the dipole moment of molecular ions depends on origin of the coordinate system, 3 nonetheless it is calculated 4 and measured, 5 assuming conventionally origin on the ion center of mass (CM).
In many other cases, measured values of physical properties for molecules responding to electromagnetic perturbations are invariant of translation. Nor do computations depend on origin of the reference frame chosen to yield corresponding theoretical estimates. At any rate, it is sometimes difficult to comprehend the meaning of the locution "origin of the coordinate system" from the context: is it origin of the multipole expansion, or that of the coordinate system of different observers placed at arbitrary reference points? By no means are they expected to coincide.
The situation is further complicated when we deal with magnetic response properties, which require specification of form and origin of the vector potential A r ð Þ. Their computational prediction is hampered by the problem of gauge dependence of results computed via approximate electronic wavefunctions. 6,7 Successful solutions have been found for static magnetizabilities 8,9 and nuclear magnetic shieldings, 10 by using basis sets of gauge including atomic orbitals (GIAO). [11][12][13] Ruud et al., 8,9 compute static magnetizabilities in the DALTON package 14 as the analytic second derivative of the energy with respect to the external magnetic field components, that is, via the magnetic Hessian evaluated at vanishing field, within coupled Hartree-Fock approximation to Rayleigh-Schrödinger perturbation theory (RSPT). 15 The advantage of the GIAO magnetic Hessian approach implemented in DALTON package 14 is that it is completely invariant of gauge origin, assumed to coincide with origin of the molecular coordinate system.
One disadvantage of this method of calculating magnetic tensors is that its implementation and application are difficult and expensive in comparison with computational schemes based on magnetic fieldinduced, origin-independent GIAO current density (GIAOCD), as shown by Keith. 16 For that reason, other authors have implemented techniques which employ electronic current density calculated via GIAO basis sets, [17][18][19] achieving substantial savings of computational effort compared with the use of magnetic Hessian within the RSPT approach, 8,9 which requires many more integral derivatives. 16 Origin-invariant current densities can also be evaluated by alternative approaches, for example, procedures allowing for a continuous set of gauge transformations (CSGT) 20 and continuous translation of origin of the current density, formally setting to zero the conventional diamagnetic term (CTOCD-DZ). 21 Nonetheless, the GIAOCD method, compared to CSGT procedures, was found to be superior for smaller basis sets. 16 At any rate, the increasing popularity of approaches employing invariant current density, either CSGT and CTOCD-DZ or GIAO, to compute magnetic-field-induced magnetic dipoles and magnetizabilities, does not seem to take into account some relevant differences with respect to the magnetic Hessian approach first pointed out by Keith 16 and examined in detail in the following, see Section 2.2, for optical perturbing fields. As a matter of fact, computational schemes using integration of the current density provide invariant static magnetizabilities only in the limit of complete basis set.
Another interesting aspect is related to the fact that paramagnetic and diamagnetic contributions to static magnetizability and nuclear magnetic shielding specified in Van Vleck 22 and Ramsey 23 theories are not uniquely defined, but exchange into one another in arbitrary gauge transformations of A, in particular in the gauge transformation equivalent to an origin shift. 24 What is remarkable, though, is that also paramagnetic parts of magnetic tensors can be separately measured: the diagonal components of magnetizability are experimentally determined by the diagonal components of the g-value tensor, within a coordinate frame with origin on the center of mass (CM), whereas paramagnetic contributions to nuclear shieldings are available from spin-rotation constants measured by microwave spectroscopy in the principal inertial axis system. 25 If values of total magnetic properties are experimentally available, then "experimental" diamagnetic contributions at origin CM are determined by difference.
At any rate, the problem of computing static magnetic properties is substantially put an end to. Notwithstanding, the theory for mixed electric dipole-magnetic dipole polarizability (MEMDP), needed for determination of optical rotation, and for frequency-dependent magnetizability (FDM), remains to some extent unresolved. 6 In fact, origin-independent estimates can be obtained for the average value of the former by adopting the dipole velocity gauge. [26][27][28] For the variational Hartree-Fock, Kohn-Sham, and MCSCF methods, translational invariance of the optical rotation, connected with the trace of MEMDP, has been achieved using London orbitals. 29,30 Origin invariant diagonal components of MEMDP can be computed via ad hoc techniques [31][32][33] and, more recently, an origin invariant rotatory power density has also been evaluated in a few chiral compounds. 34 Nonetheless, no generally accepted solution has yet been proposed to the problem of computing dynamical magnetizabilities 35,36 and related properties 37,38 via methods of time-dependent perturbation theory 39 : the extension of GIAO techniques to frequencydependent electromagnetic fields is nontrivial. No origin-independent expression of the FDM has been found in the relativistic domain. 35,36 To the best of our recollection, no visualization of dynamic current density, induced by optical magnetic fields, has so far been reported, whereas CTOCD-DZ maps started appearing. 40 The derivation of an origin-independent property which might be interpreted as "dynamic magnetizability" by Raab and de Lange 41 has been discussed, as it is not based on exact, although physically acceptable, premises. 6 In fact, it is hard to devise an experimental procedure for measuring the Raab-de Lange "ac magnetizability". 41 However, calculations have been carried out to provide theoretical estimates using approximate time-dependent density functional theory. 42 The use of time-periodic magnetic-field-dependent basis functions has also been proposed for computations of MEMDP and dynamic magnetizabilities. 43,44 The

| QUANTUM MECHANICAL APPROACH TO INDUCED ELECTRIC AND MAGNETIC MOMENTS
Using SI units and tensor notation of the previous papers, 54,55 operators for electric and magnetic, dipole and quadrupole moments of n electrons in a molecule with N nuclei, within the multipolar Bloch gauge 57 for vector and scalar potentials, are defined as follows: In these relationships, position, canonical momentum and angular momentum of the kth electron are indicated by r k , b p k ¼ Àiℏ= k and b l k ¼ r k Â b p k , respectively.
In the presence of a magnetic field B, the operator for the per- Let us introduce the general definition of n-electron densitymatrix within the McWeeny normalization 58 denoting by Ψ X ð Þ a wave function, which depends on space-spin elec- Thus, integrating over ds 1 , one gets from Equation (2.7)   59 the product of electric quadrupole operator with the electric field gradient associated with a monochromatic wave is also considered. 27,59 Accordingly, the first-order Hamiltonian can be cast in Within EQA and to first order in the electric field, electric field gradient, and magnetic field, six charge densities are considered to account for perturbation effects 54 which are conveniently expressed via vectors and second-rank tensors obtained by differentiation. For the first line of Equation (2.11), Terms appearing on the second line of (2.11) are defined in the supplementary information.
The polarization charge density (2.11) is invariant of origin in passive and active translations 54 if the hypervirial momentum is satisfied. Within the algebraic approximation, 61 this condition is met only for virtually complete basis sets.
To first-order, the induced electronic current density, 27 contains a series of terms, defined by contributions corresponding to the perturbing fields, is calculated relying on the Van Vleck common origin (CO) approach. 22 The total electronic current density (2.17) is invariant to passive and active translations 54 in the case of exact and optimal variational wavefunctions. 60 Within the algebraic approximation, 61 is invariant under a translation of the coordinate system for any basis set. 40,54 The term defines the difference between DZ and CO current densities, 54 in the ideal case of electronic wavefunctions fulfilling the off-diagonal hypervirial relationship (2.15). 60 Within the algebraic approximation, 61 this condition is exactly satisfied only for complete basis sets.
It is expedient to rewrite Equations (2.17)- (2.19) in terms of second-and third-rank current density tensors (CDT), 21,62 thus, for the first line of (2.17),

| Induced electric and magnetic dipole moments
Two basic relations can be used to compute the induced electric dipole moment,      Tables 1-3, the conversion factor to SI units is 1 a.u. of magnetizability e 2 a 2 0 =m e = 7:8910366008 Â 10 À29 JT À2 , from CODATA compilation. 66 approximately satisfied in any finite basis set calculation (gaugeless CTOCD-DZ or GIAO), which means that static magnetizability components calculated via numerical integration of (2.35) are not origin independent, unless other conditions are satisfied, for example, for symmetry reasons, vide infra.
As regards dynamic magnetic dipole moment (2.30) and related magnetizability, it is preliminarly observed that, in general, the last term of (2.35) does not vanish. However, relying on this definition to compute frequency-dependent magnetizabilities, the integrals which contribute to the last term are, for the CO approach, 54 and, for the CTOCD-DZ approach, 54 Eventually, within CO approach, the definitions arrived at via time-dependent perturbation theory are for the paramagnetic part and for the diamagnetic part of the magnetizability. The DZ contribution replacing the conventional CO diamagnetic term is given by the relation 52,69 It is worth recalling that static magnetizabilities (2.35) are computed by differentiating (with respect to the components B α and B β of the magnetic field) the energy ð2:51Þ

| Symmetry considerations
The magnetic CDT, as a function of position J Bβ α r, ω ð Þ, and the κ 0 αβ tensor are parity-odd and time-even, the Levi-Civita ε αβγ transforms as a This relation is consistent with (2.27) and (2.30), which give Equation (2.54) confirms that for D 2 symmetry the induced magnetic dipole is invariant of origin, since it is parallel to the induced electric dipole: they transform respectively as a rotation ℛ α and a translation t α belonging to the same irreducible representation.

| CHANGE OF REFERENCE POINT
The theory of polarization charge densities and induced current densities is valid under an arbitrary translation of axes, referred to as a passive transformation, 54,56 in which the molecule is maintained in the same fixed position in space. In fact, no restriction on the choice of the coordinate system for a description of these physical quantities is acceptable, since a preferential reference frame does not exist.
In other words, the representations of ρ 1 ð Þ r, ω, t ð Þ and J 1 ð Þ r, ω,t ð Þ given by two different (virtual) observers, each using his own coordinate system, K 0 and K 00 , with origins r 0 and r 00 at distance d, see Equation (2.32) and Figure 1, are both admissible. This means that the results of measurements and calculations of induced charge and current density, operated within physically equivalent K 0 and K 00 , are identical.
Instead, an active translation 54,56 keeps the axes of a fixed K standing and displaces position vectors r, scalar polarization density ρ 1 ð Þ and current density vector J 1 ð Þ , shifting them by a vector d. It is assumed that such a displacement is carried out by translation operators b T d ð Þ and b t d ð Þ, acting on different domains (in practice, the same operator), in such a way that for any function f r ð Þ or vector V r ð Þ, Þr  F I G U R E 1 For any point at distance r from a physically relevant reference point, for example, the origin of the multipole expansion, in the coordinate system K, the descriptions of induced charge and current densities of two observers in systems K 0 and K 00 , with origins r 0 and r 00 respectively, separated by a distance d, are equivalent.
Values of the induced charge density (2.11) at a point on a contour level are the same for both observers. The same conclusion is arrived at for the current density (2.17). The gauge transformation associated to the passive translation (2.32) is induced by the gauge function for a time independent magnetic field B, that is,  Dynamic magnetizability components computed for two origins, separated by 5, 3, and 1 bohr in the x, y and z directions respectively, are slightly different for compounds such as water, ammonia and hydroxilamine. In addition, difference increases on increasing ω below the first resonance frequency, as clearly documented in Figure 2.
Accordingly, we confirm that, for molecules with symmetry other than (2.52) and (2.53), the last term of (2.35), that is, (2.51), does not vanish.

| STREAMLINES AND STAGNATION GRAPHS OF CTOCD-DZ CURRENT DENSITY
Plots of asymptotic orbits of origin-independent CTOCD-DZ current density I B r, ω, t ð Þ, Equation (2.19), on the H O H σ v 0 yz ð Þ plane of the water molecule, are displayed in Figure 3. The magnetic symmetry group characterizing the electron motion, specified within Schönflies and International notation,  Two real eigenvalues of the local Jacobian matrix, with the same magnitude and opposite sign, characterize the saddle regime.
All the panels A À D in Figure 3    Static and dynamic magnetic dipole moments and molecular magnetizabilities are efficiently computed, employing high quality basis sets, by CTOCD-DZ approaches as integrals of origin-independent current density distributions induced in the electrons by timeindependent and optical magnetic fields. Provided that a grid of accurate current density values is available all over the molecular domain, numerical integration is comparatively easier than alternative methods using GIAO basis sets within the RSPT formulation, which require cumbersome calculation of matrix elements for various operators.
Computational schemes based on integration of dynamical current density vector fields, obtained by extended basis sets, are also expected to yield reliable estimates of several response tensors, for example, electric dipole and mixed electric dipole-magnetic dipole polarizabilities-virtually the same as those from conventional procedures of perturbation theory. At any level of approximation used to compute electronic wave functions, large deviations from the limit value of a given response property would only depend on insufficient basis set quality and size.
The appealing aspect of the CTOCD-DZ method is that representation of quantum mechanical electronic current density fields, via streamline and modulus maps, singular points and stagnation graphs, is easy to implement and in general very useful to rationalize molecular response to the optical fields of monochromatic light using the approaches of classical electromagnetism.

ACKNOWLEDGMENTS
Financial support from FARB 2019 and FARB 2020 is gratefully acknowledged. Open Access Funding provided by Universita degli Studi di Salerno within the CRUI-CARE Agreement.

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