Evaluation of Molecular Polarizability and of Intensity Carrying Modes Contributions in Circular Dichroism Spectroscopies

We re-examine the theory of electronic and vibrational circular dichroism spectroscopy in terms of the formalism of frequency-dependent molecular polarizabilities. We show the link between Fermi’s gold rule in circular dichroism and the trace of the complex electric dipole–magnetic dipole polarizability. We introduce the C++ code polar to compute the molecular polarizability complex tensors from quantum chemistry outputs, thus simulating straightforwardly UV-visible absorption (UV-Vis)/electronic circular dichroism (ECD) spectra, and infrared (IR)/vibrational circular dichroism (VCD) spectra. We validate the theory and the code by referring to literature data of a large group of chiral molecules, showing the remarkable accuracy of density functional theory (DFT) methods. We anticipate the application of this methodology to the interpretation of vibrational spectra in various measurement conditions, even in presence of metal surfaces with plasmonic properties. Our theoretical developments aim, in the long run, at embedding the quantum-mechanical details of the chiroptical spectroscopic response of a molecule into the simulation of the electromagnetic field distribution at the surface of plasmonic devices. Such simulations are also instrumental to the interpretation of the experimental spectra measured from devices designed to enhance chiroptical interactions by the surface plasmon resonance of metal nanostructures.


Introduction
The steady and exciting development of quantum chemistry and, notably, density functional theory (DFT) methods has greatly enhanced our capabilities in understanding molecular spectra [1]. Current quantum chemistry computer codes can be used very reliably to assign the features observed in molecular spectra to electronic or vibrational transitions, with increasingly accurate levels of theory [2][3][4][5][6][7][8][9]. Notable advances in tools available for computational spectroscopy have appeared in the recent literature [10][11][12][13][14][15][16][17][18][19]. These are commonly based on existing quantum chemistry engines, and their development often enhances the accuracy of the existing approaches. For instance, the harmonic approximation can be improved to deal with anharmonic effects [20][21][22][23]. Within the rich and evolving scenario of computational spectroscopy, it is interesting to discuss approaches complementary to those most often employed to simulate spectra, which have been in use for a long time, and basically consist in weighted sums of selected lineshape functions [24] (ultimately, the weights of such sums are related to transition moments computed by quantum chemistry-see, for instance, [25]). In fact, one of the aims of the developments of our work is the embedding of the quantum-mechanical details of the chiroptical spectroscopic response of a molecule (which requires its complex electric dipole-magnetic dipole polarizability) within the simulation of the electromagnetic field distribution at the surface of plasmonic devices. Such simulations are instrumental to the interpretation of the experimental spectra measured from devices designed to enhance chiroptical interactions by the surface plasmon resonance of metal nanostructures [26]. The role of molecular polarizability is well recognized in the field of linear and non-linear spectroscopy, and led to the development of insightful computational tools for non-linear optics (e.g., [27][28][29]). Interestingly, polarizability was also introduced in the context of electronic CD aiming to calculate CD spectra of large molecules from the polarizability characterization of fragments [30].
We have investigated polarizability from the principles of quantum mechanics and light-matter interaction [31] adopting the notation used in [24]. The paper is organized as follows. We first illustrate the connection of CD with the polarizability tensor to obtain the relevant equations which relate the quantum mechanical transition rate to the trace of molecular polarizability (i.e., Equation (9)). Thereafter we show the use of the polar code, which has been written to implement such equations, for both electronic and vibrational transitions, based on the transition matrix elements obtained from quantum chemistry codes such as Gaussian [5]. For the benchmark of the electronic polarizability we consider a set of rigid molecules, thus avoiding conformational issues (which are irrelevant to the present study). Hence, π-conjugated molecules offer a rich playground, with reliable reference theoretical and experimental data available for comparison. For vibrational polarizability, we have considered another set of small and rigid chiral molecules, which have been very well characterized in the VCD literature, both experimentally and theoretically. In the final part of this work, we show how the Intensity Carrying Modes, which were pioneered by Torii et al. in the context of vibrational polarizability and IR spectroscopy [32,33], can be smoothly extended to VCD spectroscopy as well.
Before describing our results, let us introduce the needed complex polarizabilities α ee , α em , and α mm . This notation resembles the one adopted in theoretical investigations of antennas (e.g., [34]). It is also a notation that fits computer coding, and has fostered the development of the polar program.
In the context of chiroptical spectroscopy, the α em tensor is often denotedG [35]. Similarly, the α mm tensor represents the magnetic susceptibilityχ [35,36]; the electric polarizability α ee is typically written asα [35]. We remark that the electric dipole-magnetic dipole polarizability α em was repeatedly written as sum over states in the literature [35,[37][38][39][40], even though with different purposes and different notation. In the present work, we neglect the dipole-quadrupole contribution (denotedÃ [35]), which is relevant in presence of anisotropic systems and/or gradients of the electric field. Such contribution, which can be relevant, e.g., in plasmonics, can be evaluated in perspective by a natural extension of the proposed approach. The proposed method to simulate molecular spectra is essentially a post-processing of outputs from standard quantum chemistry packages (e.g., Gaussian in the present case, but interfaces to other packages can be introduced straightforwardly). This approach is general, and it is meant to support the simulation of spectra in selected experimental conditions without the need of re-computing the molecular response from quantum chemistry. This also alleviates modifications at level of the source code of the adopted quantum chemistry engine. Moreover, the calculation of the polarizability elements is a crucial step in the implementation of effective chiral media within nano-optical electromagnetic simulations.

Results
The theoretical description of CD phenomena entails both electric and magnetic effects. In a detailed review [24], Schellman presented the theoretical description of circular dichroism and optical rotation for molecules in solution state. The central role is played by the perturbation operator V ± associated to left-handed (−) and right-handed (+) circularly polarized light, which is defined through the electric (µ) and the magnetic (m) dipole operators. By making use of the results illustrated in Ref. [24], for a given transition between states 0 and k, the expectation value of such operators is the following: where E 0 and B 0 are the electric and magnetic field amplitudes of the impinging light wave, µ 0k = 0|µ|k is the electric transition dipole moment, and m 0k = 0|m|k is the magnetic transition dipole moment. I represents the light intensity, c is the speed of light, and (a, b) indicates the complex numberz = a + ib. Schellman's review [24] adopts cgs units. For the reasons discussed by McWeeny [41], in this work we prefer using atomic units for the numerical values of polarizabilities and we consider the SI system for the description of the electromagnetic field. For the sake of simplicity, we disregard solvent effects and adopt electromagnetic waves in vacuum. Therefore, the intensity of the electromagnetic wave required by Schellman's treatment (Equation (6) in Ref. [24]) is given by where the magnetic (B 0 ) and electric field (E 0 ) amplitudes are related one another by the relation B 0 = E 0 /c [42]. The electromagnetic energy absorbed per unit time as a consequence of the 0 → k transitions promoted by the electromagnetic perturbation at frequency ω is given by [24]: where ρ k (ω) is the density of states for the 0 → k transition, and it is normalized, i.e., dω ρ k (ω) = 1. By substituting Equation (1) into Equation (2), one obtains: To connect Equation (3) with the molecular polarizability tensors, it is instrumental to write the inner products among the transition dipoles (Equation (1)) as traces of outer products: We recall that the tensor (outer) product a ⊗ b of two vectors a and b has matrix elements c ij = a i b j . Therefore, the dot product a · b can be seen as For a molecule with eigenstates |k and transition energieshω k from the ground state |0 , the electric polarizability can be written as a sum over states, with opposite sign damping Γ in the resonance term (first term) and off-resonance term (second term) [35,[43][44][45]: By recognizing in the second line of Equation (5), the numerators as the outer products of the transition electric dipoles (namely, µ 0k ⊗ µ k0 and µ k0 ⊗ µ 0k ), one can compactly write the polarizability tensor as follows: By introducing the formal substitution (µ i → µ i , µ j → m j ) in the expression of the electric polarizability given by Equation (5), one can devise the expression of the complex electric dipole-magnetic dipole polarizability α em : With a similar reasoning, the following expression can be written for the magnetic dipole polarizability tensor α mm : Near resonance conditions with the 0 → k transition, Equations (6)- (8) can be used to relate the trace of molecular polarizabilities to the transition rate given by Equation (3) (see Appendix A for details): As well-known, we recall that the imaginary part of α ee relates to absorption phenomena, and the real part to dispersive ones. Because of the presence of the purely imaginary magnetic dipole transition moment in the expression of α em , the absorption and dispersive terms in α em are the real and the imaginary parts of α em , respectively (see also Equation (A7) in Appendix A). When neglecting vibronic effects, the 0 → k transitions introduced in Equations (6)-(8) are taken as vertical electronic transitions, and they can be conveniently computed by, e.g., TDDFT. This is the approach adopted by the polar code. In summary, by the calculation of frequency-dependent polarizabilities, Equation (9) can be used: (i) to simulate UV-Vis absorption spectra in the absence of vibronic effects: (ii) to simulate electronic CD (ECD) spectra [24]: Iω and (iii) to simulate the dissymmetry factor g, which is the ratio of CD (∆ ∝ w CD ) to ordinary absorption ( ∝ w abs ) [24]: To simulate IR and VCD spectra, one has to evaluate the molecular polarizabilities by considering the 0 → k transitions as vibrational. If one focuses on the vibrational transitions in the mid-IR range, within the double harmonic approximation, the sum over states is restricted to one-quantum transitions 0 a → 1 a over the 3N-6 vibrational modes of the system, based on mechanical harmonicity. Therefore, by making the formal change (0 → k) ⇔ (0 a → 1 a ) in Equations (6)- (8), one directly obtains the desired result: wherehω a is the vibrational quantum of the a-th normal mode). Moreover, the matrix elements of the electric and magnetic dipole over one-quantum vibrational transitions can be evaluated according to the following expressions [46,47], based on electrical harmonicity assumption: In Equation (14), L ia = ∂x i /∂q a = ∂ẋ i /∂q a is the transformation matrix relating the Cartesian nuclear displacements x i to the normal coordinate q a . Wilson's L matrix [48] is often named S in VCD literature, and individual indexes are adopted to label the atom λ and the ith Cartesian component of the nuclear displacement (i = 1, 2, 3). This double index notation is cumbersome when writing loops in programs. For this reason, we prefer to merely list the Cartesian nuclear displacements along a vector with 3N components: The real vectors ∂µ/∂x i = P i collect the Atomic Polar Tensors (APTs) [49][50][51]. The P ik matrix element represents the change of the kth Cartesian component of the electric dipole caused by a change in nuclear position along the x i coordinate (i = 1.3N). The real vectors M i = ∂m/∂ẋ i collect the Atomic Axial Tensors (AATs) [49][50][51]. The M ik matrix element represents the change in the kth Cartesian component of the magnetic dipole caused by a change in the nuclear velocityẋ i . By making use of Equation (14), one obtains the following expressions for the dipole and rotatory strengths of the 0 a → 1 a vibrational transition: Further details about the evaluation of Equation (13) are given in Appendix B. The relation of molecular polarizabilities with Lambert-Beer's law and the concept of cross section is given in Appendix C.

Polarizability Due to Electronic Transitions and UV-Vis Absorption
As an initial benchmark for the polar implementation, we report in Figure 1 the results of its application on TDDFT outputs of the benzene molecule computed for increasing values of the number of excited states (N). The correctness of polar implementation is checked against independent results obtained from frequency-dependent linear response theory (i.e., CPKS equations). The latter approach is implemented in Gaussian09 and computes the polarizability avoiding the sum over states. It can be noted that the convergence of polar results to the CPKS results is smooth over a wide range of photon energies. No appreciable differences are observed when the number of states in the sum of Equation (5) is large enough (in our case N = 1000). The static limit of the polarizability [Tr(α ee (0; 0))] is known to be proportional to the volume occupied by core and valence electrons, both for atoms [52] and molecules [53]. For this reason, to reach convergence with respect to CPKS results, one needs to include enough electronic excitations in the sum over states expansion of α ee , so to make sure that also contributions from core-excitations are present. As expected, convergence is faster when approaching resonance conditions (e.g., close to 7 eV in Figure 1a), for which single transitions dominate the sum over states expression of the polarizability.  In a further independent check, we have validated polar against the reported values of the complex electric dipole polarizability of 4,4'-diaminoazobenzene [54]. A direct numerical comparison is difficult in this case due to the different quantum chemistry code used here (Gaussian09) and in [54] (ADF). To approach the computational conditions of Haghdani et al. [54], we adopted for 4,4'-diaminoazobenzene the PBE functional and the cc-pVTZ basis set. Our results, reported in Figure 2, nicely compare with those reported in Figure 3b of Ref. [54]: the maximum value of the imaginary part of Tr(α ee )/3 is about 1000 bohr 3 and the maximum value of the real part of Tr(α ee )/3 is about 600 bohr 3 . We also note that the shapes and linewidths of both the real and imaginary parts of the isotropic electric polarizability in Figure 2 are very close to those reported in Figure 3b of Ref. [54]. Tr(α ee )) computed by polar from a TDDFT calculation at PBE/cc-pVTZ level: (a) the chemical structure of 4,4'-diaminoazobenzene and the graphical representation of the computed equilibrium structure, which is essentially planar, with the -NH 2 groups slightly pyramidalized; and (b) the results obtained by the polar code for two choices of N, demonstrating good convergence at N = 100 states in this photon energy range. To approach the conditions used in [54], the damping Γ for polar calculations is to 0.24 eV.

The Electric Dipole-Magnetic Dipole Polarizability and ECD
Turning to chiral molecules and chiroptical spectroscopy, the α em tensor becomes of relevance. To validate polar for this kind of calculations, we address in Figure 3 the case of hexahelicene. This is a representative molecule of the class of helicenes, which are π-conjugated systems consisting of ortho-fused benzene rings and have many possible applications [55]. They are characterized by inherent helical chirality, and possess interesting spectroscopic properties [56][57][58][59][60][61].  As shown in Figure 3, the trace of the α em polarizability tensor evaluated by polar on the basis of TDDFT calculations nicely match with the sign and relative magnitude of the rotatory strengths R 0k reported by Gaussian09 within the same TDDFT calculation. In comparison with hexahelicene, we have considered the hexamer model of a helically coiled graphene nanoribbon of recent synthesis [62] (HGNR-6), which is also a π-conjugated system with inherent chirality (see Figure 4a). Interestingly, if we compare hexahelicene with HGNR-6 (for the same P-helicity), we obtain the same sign of the rotatory strength for the optically more active low energy π → π * transition (computed at 391 nm in HGNR-6, see Figure 4b). Notably, being π-conjugation more extended in HGNR-6 than in hexahelicene, this causes a significant red-shift of the wavelength of such transition, and higher values of its rotatory strength. This is evident by comparing Figures 3b and 4b. Furthermore, following the assignment of this low energy π → π * transition in hexahelicene to a helical sense-responsive feature (H-type band) [60], and by observing that the associated rotatory strengths have the same sign in P-hexahelicene and P-HGNR-6, we infer that such low-energy transition (computed at 391 nm) is of H-type also in HGNR-6.  We have finally considered a third inherently chiral π-conjugated molecule, the higher fullerene C 84 , which belongs to D 2 point group symmetry [63] (see Figure 5). Similar to the case of another chiral fullerene, C 76 [12,[64][65][66][67][68], UV-Vis absorption and CD spectroscopies were among the techniques which were instrumental also in the discovery and characterization of C 84 . Here, we limit ourselves to C 84 , for which the assignment of the absolute configuration was carried out by means of DFT methods in reference [69].  [69] to the experimental CD spectrum represented in the bottom part of (c). The color scheme of (a,b) helps in distinguishing pentagons from hexagons in the fullerene. The experimental CD spectrum is taken from Ref. [69]. To simulate α em (−ω, ω) (Γ = 0.35 eV), we fed polar with data from a TDDFT calculation (B3LYP/6-31G(d,p) level, 400 states).

Vibrational Polarizabilities and IR/VCD Spectroscopies
To benchmark the vibrational contributions to molecular polarizabilities provided by polar, we begin by considering α ee,v in the static limit. Such a quantity is straightforwardly provided by the Gaussian code in standard vibrational frequency calculations, and for small molecules (CH 4 , CF 4 , and CCl 4 ) it is also available from published independent calculations [70] and experiments [71]. As shown in Table 1, the numerical data provided by polar have the same quality as those reported by Gaussian, and compare reasonably with other theoretical data and experimental results. Table 1. Vibrational electric polarizabilities of CH 4 , CF 4 , and CCl 4 (data are in atomic units [41]; B3LYP/aug-cc-pVTZ DFT calculations for Columns 1 and 2-present work). The minor discrepancy between Columns 1 and 2 is due to slight differences in the numerical values of the physical constants involved in the calculations (i.e., polar in Column 1 vs. Gaussian in Column 2). This is proved by the ratios of the values in Columns 1 and 2 being 1.0072, independently of the molecule.  Parallel to the case of the electronic polarizability, where the imaginary part of α ee (−ω; ω) relates to UV-Vis absorption spectroscopy and the real part of α em (−ω; ω) relates to circular dichroism (CD) spectroscopy, the calculation of the frequency-dependent vibrational polarizabilities α ee,v (−ω; ω) and α em,v (−ω; ω) paves the way to the simulation of IR and VCD spectra, respectively. It is known that VCD and IR (to a minor extent) are fairly sensitive to structural changes in flexible molecules [72,73] and to the equilibrium among conformers [74][75][76][77]. Even though it is possible to handle molecular flexibility in VCD (e.g., see Ref. [78] for a recent approach), in our validation of polar, we decided to consider rigid molecules, for the sake of simplicity. We selected four chiral molecules with unique conformations, which offers a more straightforward comparison with the experimental IR and VCD spectra. The application of polar to molecules with many conformations is a conceptually straightforward task, which could be considered in future versions of the code. The experimental VCD and IR spectra of (  [79], and the main features of the spectra where assigned with the help of DFT calculations. The satisfactory comparison between frequency-dependent polarizabilities and experimental IR and VCD spectra of the four molecules is presented in Figures 6 and 7. The dissymmetry factor g plays an important role in VCD spectroscopy [80,81]. The g factor analytically correlates with the concept of robustness [82,83], which was introduced in VCD by Nicu and Baerends [84]. This useful parameter characterizes the soundness of the numerical prediction of VCD intensities and signs, and its use has been extended by Góbi and Magyarfalvi through the dissymmetry factor [85]. For these reasons, we have decided to compare in Figure 8 the function g(ω) computed by polar with the corresponding experimental quantity measured for 1S-CAM, taken as a representative reference. Both the order of magnitude and the spectral shape of the computed g(ω) match reasonably well with the experimental results, which we take as a solid benchmark of the implementation of the polar code. Expt. IR Expt. VCD Figure 6. (a,c) The equilibrium structure of (1S)-FEN and (1S)-CAM, respectively; and (b,d) polar-simulated IR and VCD spectra of (1S)-FEN and (1S)-CAM (Γ = 10 cm −1 ), respectively, compared to their experimental counterparts [79]. The vertical black bars in the top panels are proportional to the rotatory strengths computed by Gaussian (R 0k ). To ease the comparison with the experimental data, the wavenumber axis of the simulated polarizabilities has been uniformly scaled by the empirical factor 0.98. DFT calculations were carried out at the B3LYP/aug-cc-pVTZ level.    To further validate the calculation of the transition magnetic dipole moments and the rotatory strengths performed by polar, we have computed the transition dipole moments and the g-ratios of methyloxirane ( Figure 9), a well-known small and rigid chiral molecule which was investigated in the past by three of us [82]. The above quantities have been compared with the corresponding g-ratios that can be straightforwardly determined from the dipole and rotatory strengths extracted (in cgs units) from the typical output of a VCD calculation made by Gaussian. We collect in Table 2 the satisfactory comparison between the g-ratios computed by polar and those determined from the Gaussian output. Figure 9. Representation of the equilibrium structure of (R)-methyloxirane computed by DFT at the B3LYP/aug-cc-pVQZ level of theory. Table 2. Comparison of the results from a VCD calculation carried out with Gaussian09 on (R)-methyloxirane at the B3LYP/aug-cc-pVQZ level, and the results obtained by a corresponding polar calculation based on the numerical data contained in the *.fchk file associated with the Gaussian calculation. Column (a) is the wavenumbers computed by Gaussian09 and polar are not distinguishable at the shown precision. The APTs and AATs stored in the *.fchk file are used by polar to compute the dipole and rotatory strengths (Columns (b) and (c)) in atomic units (see Equation (14) and Section 4 for details), from which we computed the dissymmetry ratios g 0k = 4R 0k /(cD 0k ) listed in Column (d)-please note the use of SI for the electromagnetic field, which implies the factor c = 137.035999074 in atomic units [41]. The dipole and rotatory strengths computed by Gaussian09 (Columns (e) and (f)) are used to compute the dissymmetry ratios g 0k = 4R 0k /D 0k (cgs system) reported in Column (g).

Introducing Intensity-Carrying Modes in VCD
Among the observations pointed out in [79], we highlight that "the intensities, both in absorption and in VCD, are generally larger for FEN than for MEFEN and for CAM than for MECAM". Indeed, the inspection of Figures 6a and 7a confirms that the polarizabilities of the carbonyl-containing molecules FEN and CAM display larger values than MEFEN and MECAM. To directly point out the role of the carbonyl group in causing the intensity increase in both VCD and IR calculations, we compute the Intensity Carrying Modes (ICM) of the four molecules. This useful concept was introduced by Torii et al. in the framework of IR and Raman spectroscopy [32,33], connected to vibrational polarizabilities (and hyperpolarizabilities) [33], and later extended to Raman Optical Activity (ROA) by Luber et al. [86]. However, to the best the authors' knowledge, the use of ICM in the context of VCD is unprecedented. Therefore, it is worth recalling the notation used for IR spectroscopy, and extend it to VCD.
ICMs are introduced as collective nuclear displacements along which the IR intensity (or Raman or ROA intensity) is maximized. This leads to an eigenvalue problem [33,86] for the M matrix, which, in the case of IR spectroscopy, can be expressed through the dipole derivatives with respect to the nuclear displacement Cartesian coordinates (x i ; i = 1...3N): For a given M matrix (either for IR, Raman, or ROA), the eigenvalue problem provides the eigenvalues λ k associated with the maximized intensities, and the related nuclear displacements coefficients c k that describe a selected kth ICM [33,86]. If one evaluates the dipole strength D 0 a 1 a in terms of APTs (Equation (14)), one can easily justify the expression of the M matrix provided in Equation (16): We mention that Torii's approach to ICM [33] makes use of internal coordinates instead of Cartesian nuclear displacements. This is just a matter of choice in the representation of the nuclear displacements, and it does not affect the conclusions reached in this Section. Equation (18) shows that the dipole strength of a given vibrational transition (0 a → 1 a ) is a quadratic form in the coefficients of the nuclear displacement of the associated normal mode (L ia ). Clearly, Torii's approach to ICM [33] seeks the eigenvalues of the M IR matrix (Equation (16)), which is exactly the matrix of the quadratic form associated to the dipole strength reported in Equation (18). This observation justifies a similar approach for IR and VCD, provided that a suitable M matrix is defined also for the latter case. Following the strategy described for IR, to define ICM in VCD one may proceed by looking for a quadratic form associated with the rotatory strength R 0 a 1 a . By making use of Equation (14), one can express the rotatory strength in terms of APTs and AATs as follows: The last identity in Equation (19) provides the expression of the quadratic form suitable for VCD, characterized by a matrix ∂µ ∂x i · ∂m ∂ẋ j . We note that this matrix is not symmetric. However, the result of the summation in Equation (19) cannot be distinguished from the symmetric version below, which is preferred for numerical convenience: with From the inspection of Equation (20), we conclude that ICMs suitable for VCD can be found by seeking the stationary values of the rotatory strength R 01 as a function of the nuclear displacement coefficients (subject to normalization of the vector of coefficients). Similar to the case of IR spectroscopy [33], this leads to the eigenvalue problem for the M matrix defined by Equation (20).
We implemented the calculation of ICM in polar, for both IR and VCD spectroscopy, considering the M matrices given by Equations (18) and (20). The eigenvalues resulting from the solution of Equation (17) are collected in Table 3 (for IR) and in Table 4 (for VCD). We observe that the number of non-zero eigenvalues of M IR is three (as previously found by Torii [33]), while it is six for M VCD . As pointed out by Luber et al. in discussing Raman and ROA ICMs [86], the rank of the M matrix indeed equals the number of independent components of the polarizability. Similarly, in VCD, we have three independent components of the electric and magnetic dipoles, which justifies the rank of M VCD being six. Furthermore, since the sum of the rotatory strengths is zero [47] (i.e., ∑ a R 0 a 1 a = 0), this justifies the fact that the sum of the eigenvalues is approximately zero, within the numerical accuracy affordable by the selected computational method.
The inspection of the data reported in Table 3 fully justifies the observation made in Ref. [79] that the IR intensities of FEN and CAM are overall stronger than MEFEN and MECAM (respectively). The eigenvalues of the IR-ICMs of the carbonyl-containing compounds FEN and CAM are quite similar, but they are remarkably larger than those of MEFEN and MECAM. Moreover, the inspection of the animations of the associated eigenmodes (Supplementary Materials) reveals that for λ 1 and λ 2 (which are the most IR-active ICMs of FEN and CAM) the carbonyl group is significantly involved: mode 1 has a large C=O stretching content, whereas mode 2 shows a significant C=O bending content. It is thus clear that the presence of carbonyl explains the stronger activity in the IR of FEN and CAM compared to MEFEN and MECAM (as expected).
Turning now our attention to VCD, the ICM eigenvalues reported in Table 4 reveal, similar to IR, the presence in FEN and CAM of two significant positive eigenvalues (λ 1 , λ 2 ), which pair up with their negative counterparts λ 6 , λ 5 . The inspection of the associated eigenmodes (Supplementary Materials) reveals also for VCD the involvement of stretching/bending carbonyl motions occurring in-phase (λ 1 , λ 2 ) or out-of-phase (λ 6 , λ 5 ) with respect to skeletal distortions of the rest of the molecule. Therefore, based on the analysis of VCD-ICMs of the four molecules, we can provide a clue for the IR and VCD strengths of FEN and CAM being stronger than MEFEN and MECAM, as observed in [79]. Indeed, the direct involvement of carbonyl motions in the stronger ICMs is responsible for the observed differences in both IR and VCD.  Table 4. VCD Intensity Carrying Modes (ICM)s of the four chiral compounds for which the corresponding simulated and experimental spectra are reported in Figures 6 and 7. The APTs and AATs required for the evaluation of ICMs have been computed at the B3LYP/aug-cc-pVTZ level of theory. All numerical data are given in the atomic units of the M VCD matrix.

Materials and Methods
All the polarizability calculations reported here were carried out through the polar code developed in this work. polar allows the numerical evaluation of the three polarizability tensors α ee , α em , α mm as a function of the photon frequency ω. Both electronic and vibrational polarizabilities can be computed. polar is written in the C++ language, and adopts complex numbers for better code readability.
Electronic polarizabilities are evaluated as sum over vertical electronic transitions (Equations (6)- (8)). The required input data are the excitation energieshω k of a selected set of excited states |k , and the associated electric (magnetic) transition dipoles µ 0k (m 0k ). The output file of an ordinary TDDFT calculation made by Gaussian code [5] can be parsed by polar to gather all information required to operate. By writing additional parsers, polar can be straightforwardly interfaced to other quantum chemistry packages.
The calculation of vibrational polarizabilities by Equation (13) requires the quantities computed by Gaussian in a standard VCD job. These are the Hessian of the energy vs. the Cartesian nuclear displacements, and the derivatives of the electric and magnetic dipole vs. the Cartesian nuclear displacements (the latter two quantities are also, respectively, named APTs and AATs, see for instance [49][50][51]). By using the Hessian matrix, polar solves the secular equation to obtain the vibrational frequencies (ω k ) and the normal modes, which are required to compute the dipole derivatives vs. normal coordinates (q k ) from APTs and AATs. The Hessian, APTs, and AATs are stored by Gaussian in the formatted checkpoint file (*.fchk) and are read by polar before addressing the evaluation of vibrational polarizabilities. For convenience, since Gaussian checkpoint and output files report computed data in atomic units, we adopt the same convention in polar. The quite useful and neat paper by McWeeny [41] can help the reader in converting all the quantities reported here from atomic to SI units. We merely remind here the following identities (where E 0 = e 2 /(κ 0 a 0 ) is the symbol for the Hartree, a 0 the symbol for the Bohr, and κ 0 = 4π 0 is the electric permittivity): The DFT and time-dependent density functional theory(TDDFT) calculations reported in this work have been carried out with Gaussian09 Rev. D.01 [5]. Because of the scattered nature of the literature reference data employed in this work to validate polar, we had to chose individual functionals and basis sets, as best suited to carry out the required comparison. Therefore, in the caption of each figure or table presenting computational data, we report the selected functional and basis set.

Conclusions
The polar code makes the simulation of CD spectra straightforward through the evaluation of the trace of the frequency-dependent electric dipole-magnetic dipole molecular polarizability, α em (ω). polar allows simulating both ECD and VCD spectra based on the quantum chemical outputs routinely produced by the Gaussian code [5]. These calculations allow for the validation of polar thanks to the available literature data.
As a remarkable byproduct, our theoretical analysis of the molecular polarizabilities also allows one to easily compute the dissymmetry factor g(ω), which compares well with the experimental counterpart of camphor, as reported in Figure 8. Finally, we have extended Torii's Intensity Carrying Modes [32,33] to the context of VCD spectroscopy. This provides an informative approach to the understanding of the origin of VCD intensities.
Moreover, the complex polarizabilities computed by polar are of relevance to spectroscopic applications in the presence of controlled distributions of electromagnetic fields that can be realized via nano-optical engineering in order to enhance the chiroptical response [87][88][89]. In this framework, the molecular calculations can serve both the purpose of introducing the chiral medium into the electromagnetic simulations via effective medium approaches and of evaluating the weight of the quadrupolar contributions in the CD spectra because of the presence of specific field gradients.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix B. Polarizability as a Sum over Vibrational States
By evaluating the electric dipole transition elements showing up in Equation (13) through the use of Equation (14), one obtains the vibrational contribution to the electric polarizability as follows:

Appendix C. Relation with Lambert-Beer law
For convenience, we recall here the connection between the transition rate w(ω) and the Lambert-Beer law by which absorption experiments are usually interpreted (for consistency with the rest of this work, we adopt SI units, which, when needed, consistently lead to atomic units [41]. For this reason, instead of the molecular absorbivity ε, we use here the molecular cross section σ (whose dimension is an area): A is the absorbance, l the sample length, and C represents the number of absorbing molecules per unit volume (number concentration). The differential form of the Lambert-Beer law, corresponding to the differential absorption dA caused by light traversing the differential path dl, is the following: Hence, the intensity loss due to absorption within the path dl is: We remark that −dI represents the energy loss per unit time over the cross section S of the light beam, and it is caused by the rate of energy absorbed by the molecule, w(ω), times the number of molecules dN = CdV within the volume dV = Sdl, divided by the cross section S of the light beam: