Vibrational effects on polarizability: insights from normal mode analysis

Abstract

The performance of sensing and optical devices is closely linked to the molecular polarizability response, a property significantly influenced by the polarizability tensor’s isotropic and anisotropic components. In this study, we investigate the impact of vibrational motion on polarizability. We analyze how intramolecular vibrations affect the polar properties of molecules commonly employed in refrigerants and gas detection technologies. Our computational estimates align well with experimental data. It is revealed that vibrational effects can introduce uncertainties of up to 6% in sensing devices and 50% in optical devices due to nuclear oscillations. Furthermore, normal mode analysis highlights specific vibrational modes that enhance molecular polarizability.

Introduction

Detection devices are paramount in nanotechnology, with applications in aerospace, environmental monitoring, and medical fields¹. Gas ionizing-based detection mechanisms are particularly appealing due to their superior selectivity, durability, and fast recovery times². These devices function based on electronic scattering triggered by gas-phase molecular compounds^3,4. In this context, the search for efficient and environmentally sustainable materials for gas detection remains a constant demand.

Nowadays, 1,1,1,2-tetrafluoroethane (CF$_3$CH$_2$F) and dichlorodifluoromethane (CF$_2$Cl$_2$) are particularly relevant. CF$_3$CH$_2$F is the most abundant hydrofluorocarbon (HFC) in the atmosphere⁵, being widely integrated in refrigerant, gas detection⁶and propellant devices. Its ionization properties have been extensively studied both experimentally and theoretically^7,8,9,10,11. Additionally, CF$_2$Cl$_2$, also known as CFC-12, is a chlorofluorocarbon (CFC) historically used in similar applications, whose electronic scattering dynamics have significant implications for environmental monitoring and gas detection technologies¹². Although CF$_2$Cl$_2$ and CF$_3$CH$_2$F have had their industrial use ceased or progressively reduced, respectively, due to international regulations, they were industry workhorses for many years. And the overall atmospheric lifetime estimate is around 102 years for CF$_2$Cl$_2$ and 14 years for CF$_3$CH$_2$F⁵. In addition, it is important to stress that, because of their key importance to High-Energy Physics (in particle detectors), they could not be easily replaced. As a result, both molecules were reported as radiator gases in such detectors many years after being no longer recommended as refrigerant gases^6,13. Hence, these substances are still relevant for gas-sensing development since a greater risk of leakage is expected in old equipment.

Polarizability measures the electronic response of a molecule to an external electric field and is closely related to ionization processes¹⁴. Specifically, previous work suggests that the maximum electron impact ionization cross section ($\sigma _{max}$) has a positive linear correlation with the molecular isotropic polarizability¹⁵. The determination of $\sigma _{max}$is fundamental for understanding the interaction between electrons and molecules¹⁶, with significant implications for various fields¹⁷. In particular, the parameter greatly influences the Townsend coefficient¹⁸, an indicator of the ionization efficiency that is fundamental for detection applications¹⁹. Moreover, the anisotropic component measures the material’s optical response^20,21, determining properties such as the Kerr coefficient, which represents the birefringence strength of the Kerr effect^22,23,24,25.

Given this multifaceted relevance, numerous computational works have focused on theoretically calculating polarizabilities^{15,26,27,28,29,30}. Frameworks based on density functional theory (DFT) are particularly promising, offering reliable and scalable descriptions. The established DFT protocol evaluates the polarizability of optimized geometries. However, scattering mechanisms might be strongly influenced by changes in geometry. This is particularly true for phenomena involving light-matter interactions. A clear example of the importance of nuclear vibration is the mechanism of exciton transfer in molecules, which is mediated by virtual photon scattering^31,32,33.

For this reason, methodologies that incorporate vibrational effects on polarizability estimates have practical and theoretical significance. The impact of nuclear oscillations on molecular properties has been a longstanding subject in quantum chemistry. One of the most prominent approaches is the zero-point vibrational average (ZPVA) correction. The idea consists of expanding any electronic property of interest along the normal mode coordinates via first and second-order derivatives taken at the equilibrium geometry^34,35,36. The resulting shifts capture the effect of nuclear oscillations. Since its development, the theory has been widely applied, offering corrections to various properties and phenomena, such as polarizability^36,37,38,39, spin-spin coupling⁴⁰, optical attributes^39,41and thermodynamic properties⁴².

The Nuclear Ensemble (NE) method is an alternative approach that can be readily integrated into DFT workflows. Currently, the formalism is a central piece in obtaining absorption and fluorescence spectra beyond the Frank-Condon approximation^33,43,44,45. Based on the same principle, we have recently proposed a generalization of the reorganization energy, a crucial parameter for Marcus-type charge transfer⁴⁶. Moreover, Barbatti et al. built a machine learning-purpose database of multiple molecules through NE, where single-point properties such as gap, dipole magnitude, and isotropic polarizability were evaluated⁴⁷. Despite these advancements, it remains elusive how control over nuclear oscillations might influence the polarizability of important atmospheric pollutants such as CF$_2$Cl$_2$ and CF$_3$CH$_2$F.

In this work, we propose an integrated methodology combining DFT with the NE method to obtain vibrational-corrected polarizabilities of common molecules used in detection applications. Our results show good agreement with experimental data. Moreover, statistical analysis of vibrational effects indicates that optical and detection devices might have uncertainties up to, respectively, 50% and 6%.

Results

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.

Fig. 1

Experimental polarizabilities ($\alpha ^{exp}_{iso}$) versus the calculated mean ($\overline{\alpha }$) estimates of CH$_4$⁴⁸, CO⁴⁹, HF⁴⁹, H$_2$O⁴⁹, NH$_3$⁴⁹, CF$_4$⁴⁸, CF$_2$Cl$_2$⁵⁰, and CF$_3$CH$_2$F⁵¹. The error bars represent the Nuclear Ensemble dispersion.

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.

Fig. 2

Nuclear ensemble distributions of the polar properties. The blue curves in (a), (b), and (c) are the isotropic polarizability distributions of CF$_2$Cl$_2$, CF$_4$, and HF. The green histograms in (d), (e), and (f) correspond to the anisotropic distributions.

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.

Table 1 Nuclear ensemble statistical results for the isotropic and anisotropic polarizabilities. The column $\alpha _{iso}^{exp}$ corresponds to the experimental values of $\alpha _{iso}$. The third column contains the ensemble mean of $\alpha _{iso}$ with the respective dispersion. The next column shows the vibrational dispersion as a percentage of the mean value. The two last columns correspond to the same properties for the anisotropic component.

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.

Fig. 3

CF$_4$ normal modes effect on isotropic polarizability. In (a), the oscillator displacement heatmap is shown. (b) presents the 5th mode displacement as a function of $\alpha _{iso}$.

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.

Fig. 4

CF$_2$Cl$_2$ normal modes effect on isotropic polarizability. In (a), the oscillator displacement heatmap is shown. (b) presents the particular modes as functions of $\alpha _{iso}$. (c) and (d) depicts, respectively, the 4th and 5th modes of force vectors in the real space. The vector lengths are proportional to the displacement.

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.

Fig. 5

CFCl2 normal modes effect on anisotropic polarizability. In (a), the oscillator displacement heatmap is shown. (b) presents the particular modes as functions of $\alpha _{aniso}$. (c) and (d) depicts, respectively, the 4th and 0th modes of force vectors in the real space. The vector lengths are proportional to the displacement.

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.

Fig. 6

The maximum first-order isotropic derivatives with respect the nuclear coordinates of 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 indicate the most influential normal mode according to the visual analysis of the ensemble heatmaps.

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] Å.

Discussion

In conclusion, we estimated the polarizability of eight molecules relevant for refrigerant and detection applications. Our methodology combines DFT with the nuclear ensemble method, creating an integrated framework. The results enabled us to assess the impact of intramolecular vibrations on the isotropic and anisotropic polarizabilities. Statistical analysis of the calculated distributions showed good agreement with experimental estimates.

Moreover, vibrational effects cause significant dispersion in the measurement of polarizability properties, potentially leading to non-negligible uncertainties in the performance of optical and sensing applications. The result is particularly relevant for detection devices based on the maximum electron impact ionization cross section. Further analysis showed that traditional refrigerants are particularly robust against vibrational variability. For this reason, a possible approach to obtaining new (and efficient) molecules is to search for rigid compounds.

Methods

Computational details

Fig. 7

Some of the optimized geometries from molecules studied in this work: CF$_2$Cl$_2$ (a), CF$_3$CH$_2$F (b), HF (c), CF$_4$ (d), CH$_4$ (e), and NH$_3$ (f).

We adopted Density Functional Theory to calculate the electronic structure of the molecules. Eight molecules were considered. Figure 7 displays some of them, namely CF$_2$Cl$_2$ (a), CF$_3$CH$_2$F (b), HF (c), CF$_4$ (d), CH$_4$ (e), and NH$_3$ (d). An integrated workflow combining DFT with the nuclear ensemble method is applied to determine the molecular polar properties with vibrational effects. Figure 8 summarizes the sequence of steps. First, the molecular geometry is optimized. Next, frequency analysis is performed around the equilibrium geometry, determining the normal mode frequencies. Then, the nuclear ensemble is applied using the vibrational spectrum. The result is the calculation of N near-equilibrium geometries. This set of molecular configurations constitutes the ensemble. The next step is to perform the calculation of polarizability for each conformation. Combining the ensemble results returns distributions. Statistical analysis of the data yields the molecule’s vibrationally averaged polar properties.

Fig. 8

Workflow to calculate the vibrational averaged polar properties through the nuclear ensemble method.

All steps involving electronic structure calculations, namely geometry optimization, frequency analysis, and polarizability calculations were performed using the Q-Chem software⁵⁸ version 6.1 with the level of theory $\omega$b97X-D/aug-cc-pVTZ. The reliability of the optimization was confirmed by ensuring the absence of imaginary frequencies. The nuclear ensemble method is implemented through the NEMO⁴⁴ python package version 1.1.0.

Polar properties

The molecular electronic cloud rearranges when exposed to a static electric field ($\vec {F}$). The electron-field interaction associated with this relaxation can be represented as a Maclaurin series^27,49,59

$$\begin{aligned} E(\vec {F}) = E_{0} -\sum _{i=1}^{3}\mu _i F_i + \frac{1}{2!} \sum _{i=1}^{3}\sum _{j=1}^{3}\alpha _{ij}F_i F_j +..., \end{aligned}$$

(1)

in which $E_{0}$ is the zero-field energy, $\mu _i$ is the i-th component of the dipole moment and:

$$\begin{aligned} \alpha _{ij} \equiv \frac{\partial ^2 E}{\partial F_i\partial F_j}\bigg |_{\vec {F} = 0} \end{aligned}$$

(2)

is the element (i, j) of the polarizability tensor, $\overset{\leftrightarrow }{\alpha }$.

Naturally, the tensor components depend on the choice of basis. Therefore, it is convenient to compute scalar invariants such as the isotropic polarizability:

$$\begin{aligned} \alpha _{iso} = \frac{1}{3}\text {Tr}\overset{\leftrightarrow }{\alpha } \end{aligned}$$

(3)

and the anisotropic polarizability²¹:

$$\begin{aligned} \alpha _{aniso} = \frac{1}{\sqrt{2}}\sqrt{(\alpha _{X}-\alpha _{Y})^2 + (\alpha _{Y}-\alpha _{Z})^2 + (\alpha _{Z}-\alpha _{X})^2}, \end{aligned}$$

(4)

where $\alpha _{X}$, $\alpha _{Y}$, and $\alpha _{Z}$ are the principal components (eigenvalues) of $\overset{\leftrightarrow }{\alpha }$.

The determination of $\alpha _{iso}$ and $\alpha _{aniso}$ characterizes the response of the molecular wave function in the presence of the external field. Specifically, an empirical relationship between $\alpha _{iso}$ and the ionization cross section ($\sigma _{max}$) has been reported¹⁵:

$$\begin{aligned} \sigma _{max} = (1.478) \alpha _{iso}. \end{aligned}$$

(5)

In turn, the anisotropic component is relevant for optical-related phenomena such as the Kerr effect^20,21, a field-induced birefringence mechanism. The scaling parameter that controls the amplitude of the refractive index shift, the Kerr coefficient, depends on $\alpha _{aniso}$^22,23,24,25. Therefore, anisotropic polarizability estimates the material’s performance for certain applications in photonics. The numerical calculations of $\alpha _{iso}$ and $\alpha _{aniso}$were performed via the 1st-order finite difference with analytic dipoles, as discussed in the Q-Chem documentation⁵⁸.

Nuclear ensemble method

The nuclear ensemble method is a procedure to sample molecular geometries from the excitation of normal mode oscillators. The nuclear potential near equilibrium can be treated as a harmonic interaction. In this formulation, the motion can be seen as a set of coupled harmonic oscillators. This complex picture can be replaced by $N_{osc}$ decoupled harmonic oscillators after a convenient basis transformation that diagonalizes the Hessian matrix. That is known as the normal mode representation. Exciting these decoupled oscillators corresponds to obtaining new geometries around the reference configuration.

Consider a quantum harmonic oscillator with frequency $\omega _i$ and reduced mass $m_i$. It can be shown that the probability of measuring this particular mode with a distortion $R_i$ in a thermal bath at temperature Tis⁴⁴:

$$\begin{aligned} \rho _i = \left( \frac{\omega _i m_i}{\pi \hbar } \tanh \left( \frac{\hbar \omega _i}{2k_bT}\right) \right) ^{1/2} \exp \left( -\frac{\omega _i m_i}{\hbar } R_{i}^2 \tanh \left( \frac{\hbar \omega _i}{2k_bT}\right) \right) , \end{aligned}$$

(6)

in which $k_b$ is the Boltzmann constant and $\hbar$ is the reduced Planck constant.

Since the oscillators are decoupled, the probability of obtaining, simultaneously, the first, the second, ..., the $N_{osc}$-th normal mode oscillators distorted by $R_1$, $R_2$, ..., $R_{N_{osc}}$ is

$$\begin{aligned} P(R_1,R_2,...,R_{N_{osc}}) = \prod _{i=1}^{N_{osc}} \rho _i. \end{aligned}$$

(7)

Let $(R_1,R_2,...,R_{N_{osc}})\equiv \vec {R}$ be an arbitrary set of normal mode coordinates and $(Q_1,Q_2,...,Q_{N_{osc}})\equiv \vec {Q}$the corresponding shift vector which, when added to the nuclear coordinates, returns a new geometry. The nuclear ensemble method consists of four steps⁴⁴. First, sample from Eq. 7 a set of normal mode oscillator deformations $\vec {R}$. Next, reverse the transformation to real space, obtaining the corresponding shift vector ($\vec {Q}$). Apply $\vec {Q}$ to the atomic species of the equilibrium conformation. At this point, a new geometry is calculated. Finally, repeat the first step until the ensemble size is satisfactory. It is important to emphasize that the NE method is a purely harmonic approach. There is no consideration of anharmonic effects.

Fig. 9

Illustration of the nuclear ensemble method. (a) Out-of-equilibrium geometries (blue dots) are sampled around the equilibrium geometry (red dots). The set of geometries constitutes an ensemble. (b) isotropic and anisotropic polarizabilities are evaluated for each member of the ensemble, generating distributions.

Figure 9 illustrates the use of the NE method. A plot of the potential energy surface as a function of the nuclear coordinates is shown in Figure 9(a). The bottom of the curve contains a red dot, representing the optimized geometry. Note that there are multiple blue points surrounding this conformation. These are the geometries obtained by the NE method. Then, the calculation of $\alpha _{iso}$ and $\alpha _{aniso}$ returns two distributions, as depicted in Figure 9(b). In the present study, we had set $N=500$^{33,44,45,60,61} and $T=300$ K.

In the context of the NE method, the ensemble elements represent the various spatial conformations induced by vibrational effects. However, most of the physical meaning within the procedure resides on statistical estimates based on the ensemble, rather than in the individual features of each geometry. Here, the mean polarizability and its dispersion represent the material’s polar response, incorporating corrections due to nuclear oscillations. This is because the NE formalism is a time-independent sampling procedure which does not account for transitions between conformations. Each geometry has the same weight over the statistical distribution.