Vibrational effects on polarizability: insights from normal mode analysis

Before considering the findings obtained, it is important to ensure the reliability of our methodology. Assessing the suitability of the level is a crucial step when dealing with quantum chemistry simulations. This can be done by considering the estimates at different method levels. To this end, Table S1 presents the \(\alpha _{iso}\) and \(\alpha _{aniso}\) from the optimized geometry using three cardinal levels: the aug-cc-pVTZ, aug-cc-pVQZ, and aug-cc-pV5Z bases. Based on the calculations, one can note that the different bases yielded equivalent results. In particular, the isotropic component exhibits notable stability, with variations of \(\approx 1\times 10^{-2}\) Å\(^3\). The anisotropic scalar presented slightly more visible variation, although still below the NE dispersion. Therefore, the results suggest that divergences regarding basis choice do not affect the conclusions drawn based on the \(\omega\)b97X-D/aug-cc-pVTZ level of theory.

To demonstrate the convergence of the ensemble statistical metrics, we present in Figures S1 and S2 the evolution of, respectively, \(\overline{\alpha }{_{iso}}\) and \(\overline{\alpha }{_{aniso}}\) for increasing ensemble size (N). The means are normalized by the mean using 500 geometries. All compounds display qualitatively the same behavior. When the ensemble is small (\(\approx\) 100 geometries), the mean is unstable, exhibiting erratic shifts over minor additions of new geometries. However, these changes tend to gradually lose amplitude. As the ensemble reaches the \(N \approx\) 400 mark, all polarizability metrics start to display a plateau. The difference between the values with \(N>400\) and the \(N=500\) mark is significantly below 1%. Thus, further increases in ensemble size result in only minor adjustments to the mean, demonstrating convergence. Note that we had set the same N for all molecules, compounds that have different sizes and, consequently, distinct normal mode space dimensions. Because of that, \(N=500\) represents different levels of convergence for different compounds. However, it is important to realize that does not affect the proceeding analysis. All the cases exhibit significant convergence, where changes on the mean are below the 1% mark. The results in Figures S1 and S2 ensure the stability of the statistical estimates derived from the NE method.

Naturally, the application of the NE method should cause dispersion in the estimates of \(\alpha _{iso}\) and \(\alpha _{aniso}\). However, the mean of the distribution has to agree with the experimental observation for a reliable description. With this test in mind, we present in Figure 1 the experimental isotropic polarizability (\(\alpha ^{exp}_{iso}\)) as a function of the the nuclear ensemble means (\(\overline{\alpha }{_{iso}}\)) for CH\(_4\)⁴⁸, CO⁴⁹, HF⁴⁹, H\(_2\)O⁴⁹, NH\(_3\)⁴⁹, CF\(_4\)⁴⁸, CF\(_2\)Cl\(_2\)⁵⁰, and CF\(_3\)CH\(_2\)F⁵¹. A blue line with \(45^\circ\) inclination is placed to serve as a guide. The close the data is to this curve, the more accurate is the prediction. Note that the NE mean estimates provide a robust prediction. The calculated values show a close agreement with the reference values, validating the methodology.

Figure 2 presents the polarizability probability density distributions to better visualize this effect. The blue curves shown in Figures 2(a)-(c) refer to the \(\alpha _{iso}\) distributions of CF\(_2\)Cl\(_2\), CF\(_4\), and HF, respectively. The corresponding means and standard deviations are in the right corner of each plot. All distributions are mono-modal, but their center positions and dispersion vary significantly among the molecular species. Specifically, the mean isotropic polarizability of CF\(_2\)Cl\(_2\), CF\(_4\), and HF are 6.61, 2.56, and 0.82 Å\(^{3}\), respectively. For instance, the dispersion of CF\(_2\)Cl\(_2\) is 0.12 Å\(^{3}\), representing 1.8% of the mean. Next, CF\(_4\) has a variability with similar magnitude, but it corresponds to 6.7 % of the mean.

The results demonstrate that some molecules are considerably more prone to vibrational dispersion than others. That is an expected result. Photophysical attributes such as the absorption and emission spectra depend on intramolecular vibrations^33,43,44. Therefore, it is not surprising that nuclear oscillations might influence other collision mechanisms based on electronic transitions, such as the ionization process for gas sensing. Moreover, taking only the mean value may not represent the polar properties of the molecule, as some cases display considerable variation due to conformation changes. In other words, a more accurate description may benefit from accounting for the dispersion effect.

We recall that \(\alpha _{iso}\) is linearly proportional to \(\sigma _{max}\)¹⁵. Therefore, our results extend to the electron impact ionization process. Simply put, vibrational effects can induce variability in the ionization cross-section. This dispersion could represent an error in the detection range since sensing devices rely on this attribute. Importantly, some \(\alpha _{iso}\) distributions display some degree of skewness. This suggests that nuclear oscillations might lead to more profound changes in the physical description, beyond mere corrections to the mean value. Therefore, selecting molecules with high \(\alpha _{iso}\) might not be enough when developing these apparatuses. The effect of vibrations must be accounted for.

The anisotropic component exhibits even more pronounced variability. Figures 2(d)-(e), show the \(\alpha _{aniso}\) distributions for, respectively, CF\(_2\)Cl\(_2\), CF\(_4\), and HF. The anisotropic component follows an increasing trend similar to \(\alpha _{iso}\). That is, the value for CF\(_2\)Cl\(_2\) is higher than for CF\(_4\), which is higher than for HF. Specifically, the mean values of \(\alpha _{aniso}\) of CF\(_2\)Cl\(_2\), CF\(_4\), and HF are 2.10, 0.32, and 0.21 Å\(^{3}\). These values are systematically lower than the \(\alpha _{iso}\) counterparts. However, the dispersions of \(\alpha _{aniso}\) are within the same range: 0.18 (CF\(_2\)Cl\(_2\)), 0.12 (CF\(_4\)), and 0.08 (HF) Å\(^{3}\). Because of that, the overall vibrational effect becomes more relevant for \(\alpha _{aniso}\). For instance, the dispersion of the HF \(\alpha _{aniso}\) distribution represents 40.7 % of the mean. That is a considerable increase compared to the isotropic component, which showed a variation of 6.3 %. The results suggest that \(\alpha _{aniso}\) is more sensible to vibrational effects.

The distribution analysis provided important insights. Now, we extend the investigation to the remaining molecules. Table 1 presents the NE statistical results. Each row corresponds to a molecule. The second column gives the experimental values of \(\alpha _{iso}\) shown in Figure 1. The next column contains the isotropic polarizability mean (\(\overline{\alpha }{_{iso}}\)) and its dispersion \(\sigma _{iso}\). The subsequent column presents the vibrational-induced variability as a percentage of the mean value.

As indicated in Figure 2, the importance of vibrational effects varies depending on the molecule. CF\(_2\)Cl\(_2\), CF\(_4\), and CF\(_3\)CH\(_2\)F exhibit some of the lowest \(\alpha _{iso}\)vibrational sensitivities. Notably, these compounds are integrated into detecting applications⁶. This suggests that efficient detection devices rely on rigid molecules that maintain relatively consistent polar properties despite vibrational excitation.

The last two columns of Table 1 present the mean and proportional deviation of the anisotropic distributions. All molecules display significant vibrational sensitivity towards \(\alpha _{aniso}\). The case with the lowest variation is the CF\(_2\)Cl\(_2\), with a 8% dispersion. Alternatively, the most sensitive is H\(_2\)O, exhibiting almost 50% change. Traditional refrigerant agents such as CF\(_4\) and CF\(_3\)CH\(_2\)F have deviations up to 30%.

These high dispersions (when compared to the mean) show that the material’s optical response is sensitive to the nuclear oscillation. Specifically, we recall that the Kerr coefficient depends on \(\alpha _{aniso}\). Therefore, the variability of the anisotropic component extends to this optical attribute. As a result, intramolecular vibration might cause refractive index shifts, affecting the device performance^22,52.

It is also instructive to assess how different the estimates are from the standard approach (based on the optimized geometry) and the ensemble average. Tight agreement between ensemble and optimized values is observed in most cases, as demonstrated by Table S1. However, a notable exception is CH\(_4\), where the optimized estimate yields \(\alpha _{aniso}^{pVTZ/pVQZ/pV5Z} = 0.00\) Å\(^3\) while the NE average gives \(\overline{\alpha }{_{aniso}} = 0.32\) Å\(^3\). This divergence is not entirely unexpected as optimized geometries often benefit from spatial symmetries that lower the total energy. As a consequence, these symmetry constraints might force selection rules. The excitation of normal modes can break such symmetries, thereby enabling transitions and revealing electronic cloud profiles that would otherwise remain forbidden. When this happens, there is a significant difference between the methods. This result evidences the importance of integrating nuclear oscillations when estimating certain electronic properties. In addition, the result validates the expectation that non-linear optical mechanisms, such as the Kerr effect, might have significant vibrational contributions³⁴.

The NE distributions suggest that certain normal modes might significantly influence the \(\alpha _{iso/aniso}\) means. Identifying these specific oscillators could reveal chemical substitutions that enhance their excitation, thereby increasing polarizabilities. In Figure 3, we illustrate the impact of all CF\(_4\) normal modes on the isotropic polarizability. The normal mode distortion heatmap of CF\(_4\) is shown in Figure 3(a), where the modes are ordered by ascending mass-weighted frequency. The colors indicate the normal mode displacement for that particular geometry. Regions in blue refer to configurations with contracted oscillators (negative displacement), while those in red expanded (positive displacement). The 500 geometries are ordered in ascending \(\alpha _{iso}\).

Most oscillators exhibit no apparent pattern when increasing \(\alpha _{iso}\). The exception is the 5th mode. This oscillator has a clear monotonic evolution, transitioning from a contraction state (blue region) to an expanded configuration as \(\alpha _{iso}\) grows. In other words, there is a positive correlation between this mode and the isotropic polarizability.

Figure 3(b) shows this particular mode in detail, presenting the oscillator distortion as a function of \(\alpha _{iso}\). As seen, the relationship is approximately linear. This suggests that exciting this vibrational mode can shift the mean estimates. The extreme \(\alpha _{iso}\) values are \(\approx\) 2.7 and 3.0 Å\(^{3}\). Therefore, controlling the 5th mode can induce variations of at least 0.3 Å\(^{3}\), about 10% of the mean.

Now, a question arises: what modifications can be made to CF\(_4\) to favor the excitation of the 5th mode? The inset in Figure 3(b) presents the force vectors of this particular oscillator. Their lengths are proportional to the dislocation. The mode consists of a breathing motion that moves the fluorine atoms away from the central carbon when positively excited. This result suggests that polarizability enhancement is possible through bond length stretching. With this in mind, we can propose potential changes. For instance, annexing carbon-based chains in the molecule could improve steric repulsion. Another change is to replace the atomic species with similar elements that have larger radii. Doping the molecules to alter the local charge net might also favor the mechanism by improving the Coulomb repulsion. While all these options are promising, they can trigger other effects we are not anticipating due to the complexity and nuances of molecular interactions. The modifications serve as guidelines to be tested through detailed investigations.

The same analysis is applicable to the other molecules. In most cases, the contribution of a single mode stands out compared to the others. Supplementary Figures S3, S4, S5, S6, S7, and S8 present, respectively, the effects of normal modes in \(\alpha _{iso}\) for CF\(_3\)CH\(_2\)F, CH\(_4\), CO, H\(_2\)O, HF, and NH\(_3\). The colors hold the same meaning as Figure 3. Increasing the separation of atomic species seems to increase \(\alpha _{iso}\) for all molecules. This can be achieved by either contracting or expanding specific normal modes. For instance, a positive displacement of the pertinent oscillator in H\(_2\)O, CO, and HF causes an overall separation of the atomic species. In turn, the atoms in CF\(_3\)CH\(_2\)F, CH\(_4\), and NH\(_3\) separate when the modes are contracted.

It is no surprise that CH\(_4\) and NH\(_3\) share a similar response since both have pyramidal-like symmetry. Additionally, CF\(_3\)CH\(_2\)F can be seen as the combination of two pyramidal blocks: a CF3 with a CFH2. A similar pattern can be observed when looking at the corresponding heatmaps of CO and HF. Both belong to the same point group symmetry. Moreover, excitation of their relevant modes separates the two atoms. This finding suggests that the most prominent normal modes are related to the molecule’s point group symmetry.

The extreme values obtained by the NE might be accessible through precise control of the normal modes that correlate with \(\alpha _{iso}\). Therefore, it is instructive to evaluate the potential gain (or loss) of \(\alpha _{iso}\) using this design approach. One way to estimate this is by calculating the maximum and minimum values of the isotropic polarizability from Supplementary Figures S3-S8. We report that the extreme variations in CF\(_3\)CH\(_2\)F, CH\(_4\), CO, H\(_2\)O, HF, and NH\(_3\) could represent almost 50% of the mean obtained from the unrestricted sampling of normal modes. Therefore, active control of normal mode activation might be a powerful tool for creating tailored responses for optical sensing applications.

Interestingly, the interplay of two modes influences CF\(_2\)Cl\(_2\) \(\alpha _{iso}\). Figure 4(a) presents the corresponding oscillator displacement heatmap. Note that the 4th and 5th modes exhibit monotonic behavior with respect to \(\alpha _{iso}\). A dashed rectangle highlights them. The fourth oscillator expands as the isotropic polarizability increases, while the fifth contracts.

Figure 4(b) depicts the mode distortions versus \(\alpha _{iso}\). Note that the 4th mode is compressed at low \(\alpha _{iso}\), while the 5th is expanded. The profiles are reversed when switching to high \(\alpha _{iso}\): the fourth oscillator expands, and the fifth shrinks.

For some molecules, it is not feasible to visually determine the highly correlated modes. This is evident in Supplementary Figure S9, which shows the normal displacement heatmaps of CF\(_3\)CH\(_2\)F (a), CF\(_4\) (b), and NH\(_3\) (c). These results suggest that the role of vibrational effects on \(\alpha _{aniso}\) is more complex, involving non-linear combinations of certain oscillators. Conversely, the normal mode of the diatomic configurations presents a monotonic relationship with \(\alpha _{aniso}\). Supplementary Figures S10 and S11 demonstrate this outcome for, respectively, CO and HF. This is expected, as there is no other oscillator in diatomic complexes to dilute the contribution to the geometry shift.

CH\(_4\) exhibits a behavior different from all molecules. As shown in Supplementary Figures S12(a) and S12(b), the increase of \(\alpha _{aniso}\) is not due to a monotonic change in the oscillator displacement. Instead, it appears that the anisotropic component is influenced by the magnitude of the oscillator distortion. Recall that the potential oscillator energy of an oscillator with displacement \(\Delta R\) is \((1/2)m\omega ^2\Delta R^2\). Therefore, \(\alpha _{aniso}\) is favored when the mode is excited, regardless of whether it is contracted or expanded. In other words, the oscillator energy affects the anisotropic component. Following this reasoning, a potential route to maximize this optical attribute is to raise the temperature of the thermal bath. That way, more thermal energy is transferred to the oscillator, raising its amplitude.

Nowadays, a substantial number of promising molecules are being developed for sensing and optical applications^{53,54,55,56,57}. In principle, the proposed protocol can be extended to them. Here, we opted for a more conservative set of molecules for two main reasons: 1) they have an extensive track record of experimental and theoretical data, facilitating interpretation and validation; 2) the proposed normal mode analysis is only feasible for small molecules. That said, the methodology could be evoked to novel compounds, provided they are small. Treating larger compounds might be challenging since the importance of individual modes would be diluted. In such cases, the difference made by a single mode on \(\alpha\) would be subtle. In other words, at first glance, the effect of one mode might not differ significantly from the others. Instead, a particular set of oscillators will be responsible for the polarizability variations. Because of this dilution, identification of the modes could not be viable through visual confirmation, as done in Figures 4 and 5.

We recall that polarizability must be differentiable along normal mode changes. This attribute is preserved in our calculations. It is important to clarify that the oscillations seen in, for instance, Figures 4 and 5 are not the result of numerical instability. All SCF calculations were converged based on the standard threshold cutoff of the direct inversion in the iterative subspace (DIIS) metric, where the state is accepted when DISS \(< 1\times 10^{-8}\). The variations in \(\alpha\) emerge because the NE method excites simultaneously all modes to generate the out-of-equilibrium conformations. Because of that, even data points with arbitrary close oscillator displacements for a specific mode (such as in Figures 3 and 4) can have quite different \(\alpha _{iso}\) values depending on how different the remaining oscillators are excited.

In general, differentiability of \(\alpha\) is only assured when all but one mode remain inactive. To demonstrate this, we present in Figure S14 the effect on \(\alpha _{iso}\) when exciting a single mode for CH\(_4\) (a), CF\(_2\)Cl\(_2\) (b), CO (c), HF (d), CF\(_3\)CH\(_2\)F (e), H\(_2\)O (f), NH\(_3\) (g), and CF\(_4\) (h). The colors of each curve assign the mode that is being active. There are two types of response: linear and parabolic. Here, all the normal coordinates with the highest (isolated) change on \(\alpha\) have a linear dependence. Moreover, it appears that the normal coordinate associated with the highest slope is the one that influences \(\alpha _{iso}\) the most. For instance, CF\(_4\)’s most prominent mode, according to Figure 3, is the 5th, which is the same oscillator responsible for the highest slope in Figure S14(h). A similar pattern is observed for the other compounds.

This suggests that an alternative way of determining the most influential modes is by finding the one with the highest derivative. Although being a good first-approach assessment, this strategy may not be optimal in every case. The profiles of CF\(_2\)Cl\(_2\) and CF\(_3\)CH\(_2\)F illustrate a potentially overlooked feature. Here, there are two modes that have comparable slopes. Restricting to finding a single oscillator would inevitably leave one of the modes (with equivalent relevance) away from the analysis. Other limitations are more visible when considering the corresponding panel for \(\alpha _{aniso}\).

Figure S15 presents the effect of exciting one mode on the anisotropic component. While the linear and parabolic shapes remain, some different responses emerge compared to the isotropic case. The predominant modes from the molecules CH\(_4\) and CF\(_4\) exhibit a “V”-shaped curve. The presence of dominant even parity curves relates to the profile in Figure S12, where the oscillator displacement magnitude (regardless of its direction) was enough to induce the change.

Another interesting feature is the emergence of crossing events, as exemplified in Figure S15(g). Near the zero-displacement region, the 0th mode (black) produces a greater change in polarizability than the 5th mode (orange). However, the setting reverses for distortions > 0.05 Å. Beyond this threshold, the orange curve intersects the black one, eventually reaching the same variation as the 3rd mode (purple). In other words, the derivatives were not always constant, introducing a complex interplay that cannot be fully determined by the properties of the equilibrium configuration. This result shows that seeking the mode with the highest derivative can be a limiting approach in some cases, as the curves are not always linear. The mode importance might change over different displacement regimes. This is particularly relevant when considering vibrational averaged methods such as the ZPVA correction that relies on polarizability derivatives at the equilibrium geometry. Such methods may overlook crossing events, potentially electing modes as predominant that would be later dominated by other modes.

Finally, consideration of the standard formalism in light of these results is important. Both the ZPVA correction and the NE method aim to include vibrational corrections into molecular properties. These approaches share important similarities and may yield equivalent corrections under certain conditions. We recall that the ZPVA correction is relies on a Taylor expansion involving up to second-order \(\alpha\) derivatives evaluated at the equilibrium geometry. On the other hand, the NE method explicitly evaluates the polarizability at the distorted geometry. Therefore, the two methods are likely to converge when normal displacements are small.

Another situation in which the ZPVA and the NE formalism may yield similar results is when the response at the equilibrium geometry is a global characteristic, meaning that the derivatives are constant along the normal displacement axis. In other words, if polarizability varies linearly with the oscillator coordinate, the two methods should produce similar effects. Note that in this case there is no restriction to normal displacement magnitude, as the response remains invariant.

However, different corrections may emerge when away from these scenarios. One example is when the molecule is subjected to a high-temperature bath. In this case, the sampling of stronger oscillator distortions is favored. In those circumstances, the evaluation of polarizability shift through a second-order Taylor series could underrepresent the corrections. Therefore, in a sense, while an excellent first approach, the ZPVA method could not grasp all vibrational-induced change.

Regarding the derivatives, Figure S14 illustrates some possible responses, exhibiting the polarizability change when exciting one mode per time. As can be seen, the isotropic component has predominantly linear curves, where the two methods might agree. Signs of potential divergences appear for certain molecules, such as CF\(_2\)Cl\(_2\) and CF\(_4\), which exhibit parabolic or asymmetric curves instead of a linear shape. The anisotropic case (shown in Figure S15) presents even more drastic responses. Here, there are “V”-shaped curves, where the derivatives at the equilibrium configuration are not defined and show different values depending on the normal coordinate sign.

Another case that the ZPVA approach might overlook is when crossing events occur, as illustrated in Figures S14(g) and S14(b). In these profiles, the dominant mode near the equilibrium configuration can be overshadowed by another oscillator at stronger deformation regions. This mechanism cannot be captured when considering only local changes in the optimized geometry. In the face of this, the NE could be more suitable for a general-purpose protocol that includes vibrational effects on electronic properties. Therefore, although aiming to account for similar effects, the ZPVA correction and the NE method are not entirely equivalent. Each has its own advantages and limitations.

To better showcase the evident correlation between the derivatives and the most relevant mode, we present in Figure 6 the maximum absolute derivative of each normal coordinate for CH\(_4\) (a), CF\(_2\)Cl\(_2\) (b), CO (c), HF (d), CF\(_3\)CH\(_2\)F (e), H\(_2\)O (f), NH\(_3\) (g), and CF\(_4\) (h). The points in orange refer to the normal modes identified as the most influential through the previous visual analysis. Note that the highest derivatives belong to the picked modes, demonstrating complete correlation. Therefore, finding the mode with the highest derivative could serve as an alternative way of determining the most important mode, as foreshadowed in Figures S14 and S15. This approach has a clear advantage compared to the previous analysis because it does not rely on visual confirmation. If the derivatives from the most important modes remain significantly higher than those of the remaining modes, they can be automatically filtered out, simplifying the selection process.

However, the approach must be applied with caution. The anisotropic component equivalent plot demonstrates this in Figure S16. In general, the correlation remains. However, the orange points are not always at the highest peaks. Take the case of CF\(_2\)Cl\(_2\) for instance. The panel in Figure 4 shows that the 4th and 5th have comparable influence, leading to similar polarizability variations. However, that relationship does not extend to the derivative plot shown in Figure S16(b). The 0th, 6th, and 7th modes have higher derivatives, but their profile in Figure 4(a) suggests a lower influence. We attribute this discrepancy to the complex interplay of the molecular interactions when multiple normal modes are excited simultaneously. While analyzing the isolated effect of each mode can provide an excellent initial assessment, neglecting the interplay of normal oscillations to generate the shifted geometry can lead to misleading conclusions.

Another point is that taking only the maximum absolute derivative involves certain approximations. In general, the curves shown in Figures S14 and S15 do not have a constant derivative. Only a parcel of them is linear. Therefore, reducing the curve to a single value may underrepresent the global effect of the normal mode on \(\alpha\), not capturing, for instance, crossing events. Moreover, evaluating the derivative at the origin has its own limitations. The derivative of “V”-shaped curves are not defined at the origin. Therefore, the maximum derivatives in Figure 6 and S16 may be taken from any point at the range [−0.1,0.1] Å.