Deciphering divergent ferroelectric behaviour in hydrogen-bonded organic polymorphs

Abstract

Polymorphism is widely studied for tuning material properties, but is rarely exploited in organic ferroelectrics because of the requirement for polar point-group symmetry. Moreover, their strategic design demands a fundamental understanding of how molecular packing dictates polarization. In this study, we investigate two polymorphic polar forms, tetragonal and orthorhombic, of an organic enamine–imine compound using a combination of quantum crystallography and first-principles calculations. Through multipolar modeling of high-resolution X-ray and neutron diffraction data, complemented by dynamic theoretical structure factors, we confirm the presence of resonance-assisted hydrogen bonds along the proton tautomeric pathways in both forms. Our first-principles calculations indicate near-identical energy barriers to proton transfer in both cases, yet their macroscopic ferroelectric properties differ significantly. Sublattice polarization analysis demonstrates that the negligible net polarization of the orthorhombic form stems from an antiparallel arrangement of proton tautomeric pathways. These findings highlight how molecular packing and hydrogen-bond networks can be leveraged to tune ferroelectric performance, providing a blueprint for the design of next-generation organic electronic materials.

---

Introduction

Polymorphism,¹ the phenomenon of the existence of a compound in multiple crystalline forms, plays a crucial role in tuning the properties of materials, as polymorphic materials often display distinct physicochemical characteristics.^2–8 A specific property exhibited by one polymorphic form may not exist in another form. Polymorphism has thus been exploited to tune multifunctional properties, including ferroelectricity (FE),^5–10 a property that enables the switching of polarization under the influence of an applied electric field (E).¹¹ To demonstrate FE, a compound must crystallize in a polar point group, but crystallization of a material in a polar point group does not assure its ferroelectric characteristic, as the molecular packing arrangement plays a pivotal role in determining such properties. Understanding the structure–property relationship in materials offers deeper insights into their physical properties and helps understand their response to external stimuli, such as E.

For decades, the field of FE has been dominated by inorganic oxides, such as lead zirconate titanate and lead titanate, which contain toxic and heavy metals. Global concerns with toxicity and human health have shifted the focus to organic materials due to their potential for lower toxicity, environmental neutrality, and biodegradability.^12–14 Mechanisms such as order–disorder, proton transfer, displacive, and proton tautomerism have been established to understand the polarization reversal process in organic ferroelectrics.¹² Proton tautomerism mechanism (PTM)-based organic ferroelectrics exhibit a low coercive field (E_c) for efficient ferroelectric switching.¹³ Organic molecules containing keto–enol and enamine–imine as tautomeric fragments have been reported to exhibit FE via PTM.¹⁵ In their crystalline forms, a collective site-to-site proton transfer in a donor–acceptor assembly coupled with tautomerism in an intermolecular π-cooperative hydrogen bonding, also known as resonance-assisted hydrogen bonding (RAHB),¹⁶ takes place. Hence, probing RAHB is crucial for understanding PTMs in organic ferroelectrics.

The topological properties derived from multipole modeling¹⁷ of electron densities have been extensively studied to characterize non-covalent interactions.^18,19 This field has gained importance for extracting accurate features of intra- and intermolecular interactions, with wide interest in materials science.^20–23 Various intermolecular interactions have been quantified and classified for molecular systems using charge-density analysis,^18,24 and the electron density distribution has been investigated to gain insights into intermolecular RAHB.²⁵ Charge density analysis via multipole modeling of high-resolution X-ray diffraction data has been performed for the PTM-based single-component organic ferroelectric, croconic acid.²⁶ That study demonstrated the pivotal role of integral properties within the atomic basin in influencing the ferroelectric properties. An organic ferroelectric co-crystal of phenazine and chloranilic acid has also been investigated to understand the energetics and weak interactions in their ferroelectric and paraelectric phases.²⁷

Quantum mechanical calculations, coupled with modern crystallography (now called quantum crystallography),²⁸ have been synergistically exploited to provide a deeper understanding of non-covalent interactions. This approach has been employed to investigate the piezoelectric properties of perovskites²⁹ and organic molecules such as β-glycine³⁰ and croconic acid.³¹ The study on croconic acid provided insights into its detailed ferroelectric mechanism, including its electronic structure and electron density distribution in the hydrogen-bonding region. A comprehensive understanding of RAHB in hydrogen-bonded ferroelectrics requires precise locations and accurate anisotropic displacement parameters (ADPs) for the hydrogen atoms, especially for the proton responsible for generating spontaneous polarization (P_s). This can be achieved by detailed treatment of hydrogen atoms, either using neutron diffraction data or by quantum-crystallographic methods.³² A combination of high-resolution X-ray and neutron diffraction data provides accurate electron density features of hydrogen bonds, crucial for characterizing PTM-based organic ferroelectrics. Also, multipole modeling of electron density is often complemented with theoretical structure factors to circumvent systematic errors in diffraction data.^33–36

A profound understanding of PTM-based organic ferroelectrics has also been complemented by first-principles density functional theory (DFT) calculations. Such investigations aimed to estimate P_s, the minimum-energy path, and to understand the proton-transfer mechanism.³⁷ The computational studies on squaric acid and a cocrystal of the [H₅₅dmbp][Hca] salt provided a deeper understanding of the origin of antiferroelectricity by examining two sublattices with opposite polarization directions.³⁸ The study articulated the field-induced variation in the polarization of the hydrogen-bonded sublattices.

Recently, some of us have studied 2-(4-(trifluoromethyl)phenyl)-1H-phenanthro[9,10-d]imidazole (1),⁸ an enamine–imine-based single-component organic compound, to demonstrate its multifunctional properties. The compound exhibited trimorphism. One form crystallized in the tetragonal crystal system (1P_T), and the second form crystallized in an orthorhombic crystal system (1P_O), in the polar space groups P4₁ and Pca2₁, respectively. The third form crystallized in the centrosymmetric space group I4₁/a. The 1P_T form exhibited excellent FE (P_s = 4.6 µC cm⁻² and E_c = 5.8 kV cm⁻¹) and retained its ferroelectric phase well above room temperature (∼497 K); in contrast, 1P_O did not exhibit the expected FE. The P_s values estimated using the Vienna ab initio simulation package (VASP)³⁹ for the 1P_T and 1P_O forms are reported⁸ as 4.01 and 0.17 µC cm⁻², respectively (Scheme 1).

Scheme 1 Chemical structure of 2-(4-(trifluoromethyl)phenyl)-1H-phenanthro[9,10-d]imidazole (1).

Here, we probe the possible reasons for the contrasting ferroelectric characteristics of the two polar polymorphs, 1P_T and 1P_O. To confirm the RAHB characteristics of the proton tautomeric pathways in these two polymorphs, we report quantitative and qualitative analyses of electron densities derived from multipole modeling of high-resolution X-ray and neutron diffraction data, as well as dynamic theoretical structure factors generated from periodic quantum-mechanical calculations. We also investigate variations in the electron localization function (ELF) to characterize hydrogen-bonding interactions along the proton tautomeric pathway. Having established the RAHB characteristics from an electron-density perspective in both cases, we estimate the potential energy barrier for proton transfer along the tautomeric pathways. Finally, we calculate sublattice polarizations to define the mechanism behind the distinct ferroelectric characteristics.

Experimental section

High-resolution single-crystal X-ray diffraction (SCXRD) experiment

A Bruker APEX3 diffractometer equipped with D8 Venture IµS microfocus dual sources, a PHOTON II CPAD detector, and an Oxford cryogenic system was used to perform these experiments on single crystals of the 1P_T and 1P_O forms. Monochromatic Mo K_α radiation (λ = 0.71073 Å) was used to collect data at 100 K for 1P_T and 90 K for 1P_O. The crystal-to-detector distances were set at 50 and 55 mm for 1P_T and 1P_O, respectively. For 1P_T, 24 sets of runs with varying detector positions (2θ) and ω angles were collected using exposure times of 6 and 10 s for lower angle reflections (2θ = 20.69° and 35.69°) and 30, 80, and 130 s for higher angle reflections (2θ = 62.34°, −77.34°, and −92.34°). For 1P_O, 18 runs were collected using exposure times of 2, 80, and 95 s for 2θ = 0.94°, −76.86°, and −86.86°, respectively. The frame width was set to 0.5° throughout the data collection for both forms.

Single-crystal neutron diffraction experiment

A large block-shaped crystal of the 1P_T form was selected for the data collection. No suitable-sized crystal of the 1P_O form exists to date, precluding its single-crystal neutron study. The neutron experiment was performed on the KOALA-2 instrument at the Australian Nuclear Science and Technology Organization (ANSTO). The KOALA-2 instrument utilises a reactor-based neutron source, enabling the recording of neutron Laue diffraction images. For consistency with the high-resolution X-ray diffraction experiment, the neutron experiment was performed at 100 K using an Oxford cryosystems open-flow type nitrogen cryostat⁴⁰ with a nominal stability of 0.1 K. A total of 16 frames were collected with an exposure time of 267 minutes per frame. The data were reduced using the LaueG suite,⁴¹ and the structure was refined using the CRYSTALS.⁴² Due to the inherent limitations of the Laue method for accurately deriving unit cell parameters, those determined from X-ray diffraction experiments at 100 K were employed for all aspects of the neutron single-crystal study of the 1P_T form.

X-ray diffraction data analysis and structure refinement

The X-ray diffraction data reduction and integration were performed using the SAINT⁴³ program. The sorting, scaling, absorption corrections, and merging of reflections, except for the Friedel pairs, were performed using the SADABS program⁴³ implemented in the APEX5 suite. The structures were solved using SHELXT⁴⁴ and refined using Olex2.⁴⁵ The non-hydrogen atoms were located in successive difference Fourier maps and refined with ADPs. All hydrogen atoms were located and subsequently refined using a riding model, except those attached to nitrogen, which were located in difference Fourier maps. Subsequently, the Independent Atom Model (IAM)-based structures were considered for Hirshfeld atom refinement (HAR)⁴⁶ and for multipole refinement using the Olex2 and XD (Revision 2016.01, July 8, 2016) programs,⁴⁷ respectively. The hydrogen atom positions and their corresponding ADPs were estimated by the HAR method using the NoSphera2⁴⁸ program implemented in Olex2. The molecular wave function was generated using ORCA 5.0.4⁴⁹ at the BLYP/cc-pVTZ⁴⁶ level of theory, assuming vacuum as the medium, with a charge of zero and a multiplicity of 1, which was used as input to perform the refinement until convergence.

Multipole modeling of X-ray and neutron diffraction data

The multipole refinement was performed based on Hansen–Coppens multipole formalism.¹⁷ The structural information from the IAM model was imported to the XD program using the XDINI module. All X–H distances were restrained to the values obtained from neutron diffraction data. The ADPs for hydrogen atoms were derived by scaling the neutron-derived ADPs using the Blessing method,⁵⁰ as implemented in the SCALE-XN program. For the 1P_O form, in the absence of neutron diffraction data and being polymorphs, the multipole refinement was performed using the X–H distances and scaled ADPs of the 1P_T form. The XDLSM module was then used to refine the scale factor using all reflections. Next, the atomic positions and ADPs of non-hydrogen atoms were refined using high-order diffraction data (sin θ/λ ≥ 0.8 Å⁻¹). The subsequent refinements of multipoles up to hexadecapoles (l_max = 4) for fluorine atoms, up to octupoles (l_max = 3) for carbon and nitrogen atoms, and up to dipoles (l_max = 1) for hydrogen, and additionally the first component of quadrupole (l_max = 2) for the hydrogen atoms bonded to nitrogen atoms, were carried out in a stepwise manner. Twenty-three sets of kappas (κ) were assigned for the chemically different types of non-hydrogen atoms. For hydrogens, the κ values were fixed to 1.2. Finally, all these parameters were combined to perform refinement using all reflections. At each step, refinement continued until convergence was achieved. For modeling the anharmonic motions of fluorine atoms, the fourth-order Gram–Charlier (GC) coefficients were refined using high-order diffraction data.^35,51,52 Additionally, for the 1P_O form, the GC coefficients were refined for the C22A and C11B atoms. Furthermore, for both forms, κ and κ′ for hydrogen atoms and κ′ for fluorine atoms were transferred from the multipole model constructed from dynamic structure factors calculated using periodic mechanical quantum calculations, as described in the following section.

Theoretical section

Molecular electrostatic surface potential and frontier molecular orbital calculations

The initial structure of 1 was optimized at the B3LYP-D3/6-31G(d,p) level of theory, followed by vibrational frequency calculations to confirm that the optimized geometry corresponds to a true minimum on the potential energy surface. Single-point calculations were subsequently performed on the dimeric crystal structures of the polymorphs 1P_T and 1P_O, as extracted from their respective crystal lattices. Electron density and electrostatic surface potential (ESP) distributions were derived from the computed wavefunctions, and the ESP maps were visualized using MoleCoolQt.⁵³ The corresponding frontier molecular orbital (FMO) analyses were carried out using ChemCraft.⁵⁴ The quantum chemical calculations were performed using Gaussian 16.⁵⁵

Dynamic structure factor and multipole modeling

For both polymorphs, periodic quantum-mechanical calculations were performed using the CRYSTAL23 code⁵⁶ based on the set of atomic coordinates obtained from the experimental multipole model. The wave function was generated using an all-electron Gaussian-type localized orbital 6-31G** basis set⁵⁷ and generalized gradient approximation of Perdew–Burke–Ernzerhof (PBE)⁵⁸ was used as the exchange–correlation functional in combination with Grimme's empirical dispersion correction (D3).⁵⁹ The energy convergence was achieved by setting the tolerance and level shifter to 10⁻⁷ and 0.3 Hartree, respectively. Along the reciprocal lattice vectors, the shrinking factor was set to 4. The periodic wave function and the experimental ADPs were then utilized to generate the dynamic structure factor for the set of experimentally observed reflections using the XFAC module. Subsequently, the multipole modeling was performed using F², refining the same set of multipoles as in the experimental multipole model.

Estimation of proton transfer barrier

The proton transfer barrier was estimated using first-principles DFT calculations. Here, we performed the climbing image nudged elastic band (CI-NEB)⁶⁰ calculations using VASP 5.4.4 with the projector augmented wave pseudopotentials⁶¹ for the electron–ion interactions. The wave function was expanded in a plane-wave basis using a kinetic-energy cutoff of 550 eV. We used the PBE exchange–correlation functional in combination with Grimme's D3 correction and a k-point mesh of 1 × 2 × 1 and 2 × 2 × 2 for 1P_O and 1P_T forms, respectively. The convergence criteria for the energy and residual Hellmann–Feynman forces were set to 10⁻⁵ eV and 10⁻³ eV Å⁻¹, respectively. The cell parameters and atomic coordinates of the experimental crystal structures of 1P_T and 1P_O were optimized at 0 K. The corresponding optimized structure is referred to as the “pristine structure”. Subsequently, the proton of the –NH group was interpolated linearly towards the N-atom of the adjacent molecule participating in the N–H⋯N bonding. The corresponding structure is referred to as the “modified structure”, which comprises N⋯H–N hydrogen bonds. For the modified structure, only the atomic coordinates were optimized using the same cell parameters throughout the NEB trajectory. Considering the collective proton transfer during polarization reversal in the polymorphic forms, a total of five images were used in the trajectory between the before- and after-proton-transfer configurations, ensuring a step width of 0.1–0.2 Å for tracing the proton transfer. The entire NEB trajectory is optimized by selectively allowing proton movement in the plane perpendicular to the hydrogen-bonding direction. The energy differences between the optimized structures before and after the proton transfer result in different barrier heights in the forward and backward directions. Therefore, the potential barrier,⁶² i.e., the kinetically resolved activation (KRA) barrier, was estimated upon averaging the potential barrier in the forward and backward directions and using the following relation:

ΔE^KRA = E^TS − ½(E^I + E^F)

where E^TS is the energy of the transition state and E^I and E^F are the energies before and after the proton transfer configuration, respectively.

Sublattice polarization calculation

For the 1P_O form, we calculated the contributions of the two sublattices to the net P_s using the modern theory of polarization based on the Berry phase method,^63,64 as implemented in the VASP 5.4.4 code, with the same level of calculation described above. The polarization calculation was performed using the strategy for estimating the polarization of PTM-based organic ferroelectrics.^5,8,13 However, in this case, the polarization was calculated layer by layer by removing the second layer from the unit cell. Additionally, the polarization of the two layers in the unit cell was also estimated at 0 K using the optimized geometry.

Results and discussion

Crystal structure analysis

Excellent quality crystals of the 1P_T form were grown via slow evaporation from a solution in m-xylene, and those of the 1P_O form were grown via vacuum sublimation. The superior quality of X-ray diffraction from the 1P_T crystal resulted in slightly higher resolution data compared to the 1P_O crystal (Table 1). The polymorphs 1P_T and 1P_O crystallize in the polar space groups P4₁ and Pca2₁, respectively. The neutron diffraction-based crystal data and refinement parameters of 1P_T are tabulated in Table S1. The N–H bond and N–H⋯N hydrogen bond distances from neutron diffraction data were compared with those from the HAR model and with the standard neutron distances⁶⁵ (Table S2). While the standard N–H neutron distances are comparable to those determined by neutron diffraction, the HAR model-based distances are slightly underestimated. Consequently, the N–H⋯N hydrogen bond distances from the HAR model are marginally elongated. Due to the limitations of X-ray diffraction for hydrogen-atom site determination, the proton tautomeric pathway for both polymorphic forms was further analyzed after resetting the N–H distances to the values obtained from the neutron diffraction experiment. The displacement ellipsoid plots based on the multipole model of the polymorphic forms are shown in Fig. S1. Both forms contain two crystallographically independent molecules in the asymmetric unit (Z′ = 2), labeled A and B.

Table 1 Crystal data and refinement parameters of the two polymorphs

Polymorphs	1P_T	1P_O
Chemical formula	C22H13F3N2
Formula weight	362.357
Space group	P41	Pca21
a, b, c (Å)	14.4058(5), 14.4058(5), 15.8555(6)	18.6560(8), 10.2687(4), 17.9704(7)
Volume (Å^3)	3290.4(2)	3442.6(2)
Z; Z′	8; 2	8; 2
Resolution (dmin) (Å)	0.44	0.46
Density (g cm^−3)	1.463	1.398
µ (mm^−1)	0.111	0.106
Tmin, Tmax	0.973, 0.980	0.951, 0.991
Rint, Rmerge	0.0456, 0.0450	0.0543, 0.0541
Measured reflections	452 320	363 650
Unique reflections	40 552	37 127
Completeness (%)	99.9	99.8
Redundancy	11.13	9.55

---

Model	IAM
Reflection no. [I > 2σ(I)]	36 599	34 169
R(F^2), Rw(F^2)	0.0315, 0.0845	0.0346, 0.0906
Goodness of fit (S)	1.116	1.046
Δρmax, Δρmin (e Å^−3)	0.595, −0.347	0.690, −0.695

---

Model	HAR
Reflection no. [I > 2σ(I)]	36 599	34 169
R(F^2), Rw(F^2)	0.0219, 0.0416	0.0268, 0.0625
Goodness of fit (S)	1.1205	1.034
Δρmax, Δρmin (e Å^−3)	0.526, −0.330	0.571, −0.732

---

Model	Multipole
	Experimental	Theoretical	Experimental	Theoretical
Reflection no. (Nref) [I > 3σ(I)]	35 571	40 609	33 592	39 066
Nref/Nv	23.0	38.0	21.7	36.5
R(F^2), Rw(F^2)	0.0204, 0.0259	0.0055, 0.0021	0.0231, 0.0321	0.0050, 0.0018
Δρmax, Δρmin (e Å^−3)	0.148, −0.145	0.093, −0.128	0.202, −0.206	0.092, −0.102

---

Molecular packing analysis

The polymorphs adopt an ABAB type structure, and the A and B molecules form butterfly-like arrangements via asymmetric N_A–H_A⋯N_B and N_B–H_B⋯N_A hydrogen bonds along the c- and a-axes for 1P_T and 1P_O, respectively (Fig. S2). Furthermore, the packing analysis reveals that the difference in distance between the two successive N–H⋯N hydrogen bonds in 1P_T (0.047 Å) is shorter than that in 1P_O (0.161 Å) (Fig. 1). This suggests that the A-type molecules are displaced from the center position between the two B-type molecules in 1P_O by a greater amount than that in 1P_T. The molecular packing of both polymorphs was analysed and discussed in detail in an earlier report by some of us.⁸

Fig. 1 Depiction of N–H⋯N hydrogen bonds between the A and B type molecules in 1P_T (left) and 1P_O (right). Only H-atoms of the –NH groups participating in the RAHB are shown for clarity. The distances are in Å.

Scrutiny of multipole models

The final experimental and theoretical multipole models derived for both polymorphs were tested for accuracy based on the statistical analysis of the residual electron densities. The anticipated parabolic shape of the fractal dimension plot and the Gaussian distribution of the pixel population, plotted against the residual electron densities (Fig. S3 and S4), indicate the accuracy of the models. A maximum scale-factor variation of 5% from unity, as observed in the DRK plots (Fig. S5), suggests that these models are of high quality. The featureless residual electron densities (Table 1) confirm that the experimental (Fig. S6 and S7) and theoretical (Fig. S8 and S9) multipole models are accurate. It is to be noted that the slight negative electron densities at some of the atomic positions in the theoretical residual electron density maps are because the molecular computation in the CRYSTAL23 package is based on Gaussian-type orbitals centered on each atomic site, whereas Slater-type orbitals are used to define the core electron density in the multipole refinement.⁶⁶

Deformation electron density

The sharp electron density features of the covalent bonds and the lone pair of electrons of the N-atoms in both polymorphic forms are highlighted via two-dimensional (2D) and three-dimensional (3D) deformation electron density maps (Fig. S10–S15). However, the experimental 3D deformation electron densities of the –CF₃ groups of both forms appear slightly diffused due to the high anharmonic thermal motions of the F-atoms (Fig. S16). Nevertheless, they agree well with those plotted based on the theoretical multipole model, which display sharp electron density features (Fig. S17). The 3D maps along the proton tautomeric pathway in both forms reveal that the non-bonding electron density of sp² N-atoms exhibits different degrees of polarization towards the H-atoms of the adjacent molecule (Fig. S18). While in 1P_T N2A and N2B display a similar degree of polarization of their lone pairs (Fig. S19), in 1P_O the lone pair of N2B is more polarized than that of N2A (Fig. S20). The nature of the degrees of polarization in the two forms correlates well with the N–H⋯N distances shown in Fig. 1.

Topological properties and Laplacian

The topological analyses of electron densities derived from the experimental and theoretical multipole models were performed using Bader's quantum theory of atoms in molecules (QTAIM) approach.⁶⁷ The topological parameters, electron density (ρ) and Laplacian (∇²ρ) at the bond critical points (BCPs), for the covalent bonds, are listed in Table S3. Overall, the ρ values of the covalent bonds of A and B molecules of both forms are slightly overestimated in the experimental model than in the theoretical model. Grossly, the values of the topological parameters for 1P_T and 1P_O are comparable. The bond paths (BPs) and the associated BCPs along the covalent bonds, the N–H⋯N hydrogen bonds, and the ring critical points for both forms have been located and are shown in Fig. S21. The Laplacian maps, shown in Fig. S22–S25, highlight the charge concentration and depletion along the covalent bonds. For the two polymorphs, the bonds involved in the proton tautomeric pathway exhibit comparable topological parameters, and the values obtained from the experimental and theoretical models correlate well (Table S4 and Fig. 2). In both cases, the higher values of ∇²ρ for the N1–C2 bonds of both A and B molecules (values ranging from −18.270 to −23.753 e Å⁻⁵ for the A molecule and from −17.766 to −25.167 e Å⁻⁵ for the B molecule) as compared to those of the N1–C1 bonds (values ranging from −15.228 to −19.255 e Å⁻⁵ for the A molecule and from −14.166 to −21.316 e Å⁻⁵ for the B molecule), suggest a high degree of charge concentration in the conjugated bonds. The higher values of ellipticity in the N1–C2 bonds of both A and B molecules (values ranging from 0.14 to 0.30 for the A molecule and from 0.15 to 0.28 for the B molecule) than those of the N1–C1 bonds (values ranging from 0.05 to 0.16 for the A molecule and from 0.11 to 0.15 for the B molecule) suggest a higher covalency in N1–C2 bonds than in a typical N–C bond. The topological parameters are also in good accordance with the values reported for 2-nitroimidazole⁶⁸ and benzimidazole.⁶⁹ The above observations confirm the presence of resonance within the enamine–imine fragment and RAHB in both polymorphic forms. For the intermolecular hydrogen bonds, the positive values of ∇²ρ suggest their closed-shell nature. While the topological parameters of H1A⋯N2B and H1B⋯N2A hydrogen bonds are comparable in the case of the 1P_T form, they are significantly higher in the 1P_O form, which is in accordance with their distinct bond path length, R_ij (Table S4 and Fig. 2). Accordingly, a significant polarization in the valence shell charge concentration (VSCC) around H1 and N2 is observed in the Laplacian maps, plotted using both experimental and theoretical multipole models (Fig. S26 and S27).

Fig. 2 Representation of ρ (e Å⁻³) (values in blue), Laplacian (e Å⁻⁵) (values in red), and R_ij (Å) (values in black) at the BCPs for bonds involved in the tautomeric path for (a) 1P_T and (b) 1P_O based on experimental and theoretical (values in italics) models.

Electron localization function

The localization of electron density in bonding and non-bonding regions of a compound can qualitatively be evaluated using the ELF map. The highest value of ELF, i.e., 1, indicates perfect localization of electron density. The 3D-ELF, shown in Fig. 3, highlights regions of delocalization and regions with perfect electron density localization. The ELF features derived from the experimental model agree well with those from the theoretical model. The 2D-ELF maps reveal that the electron densities of the N2 atom in both polymorphs are delocalized (highlighted in red in Fig. S28–S31), indicating a high degree of polarization, and the degree of delocalization is higher in molecule B than in molecule A. Furthermore, the variation in ELF within the intermolecular hydrogen-bonding region is quantified (Fig. 4). The ELF profile displays three characteristic humps; the first and the third hump in the profile correspond to the nuclei of hydrogen and nitrogen atoms, respectively, whereas the middle hump corresponds to the intermolecular hydrogen bonding region. The ELF profile indicates that the electron density is localized more toward the nitrogen atoms. Between the polymorphs, the electron density distribution in the intermolecular hydrogen-bonding region shows minimal differences toward nitrogen atoms, whereas subtle differences are observed toward hydrogen atoms. For the H1A⋯N2B hydrogen bond, a slightly higher electron localization (near the BCPs at a distance of 0.57–0.62 Å from the hydrogen atom) is observed in the case of 1P_O than that of 1P_T (Fig. 4a). The higher electron localization in the H1A⋯N2A hydrogen bond in the 1P_O form correlates well with its high electron density at the BCP and a shorter BP (Table S4). The H1B⋯N2A hydrogen bond shows similar electron localization (near the BCPs at a distance of 0.65–0.69 Å from the hydrogen atom) in both polymorphs, except for the experimental multipole model of the 1P_O form, which is attributed to its lower electron density (Table S4). A similar ELF profile was observed for the O–H⋯O hydrogen bonds in the case of the well-known hydrogen-bonded ferroelectric, KDP.⁷⁰

Fig. 3 Experimental (left) and theoretical (right) 3D-ELFs along the proton tautomeric pathways of (a) 1P_T and (b) 1P_O. The color gradient bar with scale is shown in the inset.

Fig. 4 Variation of electron localization function for the (a) N1A–H1A⋯N2B and (b) N1B–H1B⋯N2A hydrogen bonds in 1P_T and 1P_O.

Estimation of potential energy barrier

The role of molecular packing in proton transfer, assuming collective proton transfer as the mechanism for polarization reversal, was investigated by calculating the energy barrier based on both the optimized “pristine structure” and the “modified structures” (Fig. S32). For both polymorphs, the lattice parameters, N–H bond lengths, N–H⋯N hydrogen-bond distances, and N–H⋯N angles of the optimized structure agree well with those of the crystal structure (Table S5). Furthermore, the N–H bond lengths and N–H⋯N hydrogen-bond distances between the pristine structure and the modified structures compare well. The N–H⋯N angles differ by 8.86° and 11.75° for the 1P_T form, and 1.22° and 13.08° for the 1P_O form, leading to a slight energy difference between the optimized and modified structures (1.069 and 0.801 eV for the 1P_T and 1P_O forms, respectively), as represented by the last data point in Fig. 5. However, the energy barriers for proton transfer for the 1P_T and 1P_O forms are estimated to be nearly equal, 0.407 and 0.406 eV/proton, respectively, indicating that the proton transfer is indeed possible in both polymorphs.

Fig. 5 Minimum energy pathway for the collective proton transfer mechanism for polarization reversal in 1P_T and 1P_O.

Estimation of sublattice polarization

Molecular packing analyses reveal a significant difference in the alignment of proton tautomeric pathways between the two polymorphic forms (Fig. 6). In the 1P_T form, the proton tautomeric pathways are aligned parallel to the polar c-axis (Fig. 6a). In the 1P_O form, there are two sublattices formed via N–H⋯N hydrogen bonds along the a-axis, as highlighted using pink shaded planes in Fig. 6b. The angle between the mean planes containing the N–H⋯N hydrogen bonds is 11.99°, indicating that the proton tautomeric pathways are almost antiparallel, which is likely to diminish the net P_s in this form. The consequence of the nearly antiparallel proton tautomeric pathways on the net P_s was probed by estimating the polarization components for each sublattice. The estimation of polarizations, i.e., P₁ and P₂ for the two sublattices, reveals P_1x = −P_2x = −1.88, P_1y = P_2y = 0, and P_1z = P_2z = −0.05 µC cm⁻² for the crystal geometry, and those of the optimized geometry are −1.75, 0, and 0.05. A slight difference in the polarization values between the crystal (|P₁| = |P₂| = 1.88 µC cm⁻²) and the optimized (|P₁| = |P₂| = 1.75 µC cm⁻²) geometries arises from the subtle differences between the neutron and optimized N–H distances. The significant value of P_x, i.e., polarization along the a-axis, coincides with the direction of proton tautomeric pathways. However, the equal amount of polarization but in opposite directions (P_1x = −P_2x) leads to the cancellation of polarization along the a-axis. In contrast, an equal, yet small, amount of polarization along the c-axis (polar axis) in the same direction yields a negligible net P_s in 1P_O (0.17 µC cm⁻²).⁸

Fig. 6 Molecular packing diagrams of (a) 1P_T and (b) 1P_O, highlighting the proton tautomeric pathways.

Molecular electrostatic surface potential and frontier molecular orbital analysis

Furthermore, to correlate molecular packing geometry with polarization, the structural descriptors, MESP and FMO, were analyzed. The MESP maps of 1 and its polymorphic forms (1P_T and 1P_O) reveal a highly anisotropic electrostatic distribution that governs the intermolecular recognition pattern in the crystal lattice (Fig. 7). The surfaces are predominantly characterized by green-yellow regions corresponding to moderately neutral ESP (V_S), while localized intense blue domains (V_S⁺_,avg ≈ 183.79 kJ mol⁻¹) highlight strongly electron-deficient sites associated with polarized hydrogen atoms acting as efficient hydrogen-bond donors. Notably, the persistence of these positive ESP regions (V_S⁺_,avg ≈ 183.79 and 183.41 kJ mol⁻¹) in the respective polymorphic forms indicates that molecular packing, i.e., supramolecular aggregation, perturbs the intrinsic electrophilic character of the donor sites to the least. In contrast, the localized red-orange regions representing negative ESP exhibit a slight attenuation from V_S⁻_,avg ≈ −52.51 kJ mol⁻¹ in the isolated molecular form to −47.28 kJ mol⁻¹ in the aggregated state, suggesting partial electrostatic compensation arising from intermolecular charge redistribution due to the molecular packing in the crystal lattice. The spatial complementarity between the electron-rich acceptor domains and electron-deficient donor regions provides direct evidence for directional donor–acceptor interactions that promote robust intermolecular hydrogen bonding and stabilize the supramolecular framework. Furthermore, the distinct redistribution of ESP in 1P_T and 1P_O highlights the subtle influence of polymorphism on local charge organization, indicating that variations in molecular packing modulate the extent of intermolecular electrostatic coupling without substantially altering the fundamental reactive surface characteristics of the molecular scaffold.

Fig. 7 MESP maps of the (a) optimized geometry of 1 and the crystal geometry of (b) 1P_T and (c) 1P_O, plotted on the 0.02 a.u. electron-density isosurface. The average maximum (V_S⁺_,avg) and average minimum (V_S⁻_,avg) ESP values are in kJ mol⁻¹.

From a frontier electronic perspective, both 1P_T and 1P_O exhibit well-defined donor–acceptor electronic architectures with significant π-electron delocalization across the molecular backbone (Fig. 8), which is highly favorable for polarization and domain switching. In 1P_T, the HOMO is predominantly localized over the electron-rich aromatic and heterocyclic donor of one of the two molecules in the asymmetric unit of the crystal lattice, whereas the LUMO is distributed toward the halogenated acceptor region of the other molecule (Fig. 8a), indicating pronounced intramolecular charge-transfer characteristics. Furthermore, a comparatively lower HOMO–LUMO energy gap of 307.11 kJ mol⁻¹ suggests its enhanced electronic softness, higher charge mobility, and easier polarization under an external E. Similarly, 1P_O demonstrates substantial orbital delocalization, however, with slightly denser population (Fig. 8b). Its higher energy gap (313.48 kJ mol⁻¹) indicates marginally reduced charge-transfer efficiency and greater kinetic stability relative to 1P_T. Nevertheless, the spatial separation between HOMO and LUMO densities in both systems supports effective electronic polarization and directional electron migration, which are essential parameters governing ferroelectric response and electronic switching behavior. Overall, the FMO analysis suggests that 1P_T may exhibit prominent ferroelectric properties and higher polarizability, owing to its narrower band gap and stronger donor–acceptor interactions. Compared to 1P_T, the 1P_O form exhibits enhanced (6.37 kJ mol⁻¹) electronic stability with moderate polarizability, leading to a distinct molecular packing in its crystal lattice (Fig. 6).

Fig. 8 FMO distributions of the HOMO and LUMO for the asymmetric unit (Z′ = 2) of (a) 1P_T and (b) 1P_O. The orbitals are visualized at an isosurface value of 0.02 a.u.

Conclusion

This study provides a comprehensive quantum crystallographic rationale for the divergent ferroelectric behavior observed in the tetragonal and orthorhombic polar polymorphs of an enamine–imine-based organic compound. By integrating multipole modeling of high-resolution X-ray and neutron diffraction and 3D-ELF analysis, we successfully mapped the electron density distributions and confirmed the presence of RAHB along the proton tautomeric pathways in both forms. Our findings present a compelling paradox. While first-principles calculations indicate nearly identical potential energy barriers for proton transfer in both polymorphs, their macroscopic polarizations differ significantly. The resolution to this discrepancy lies in the sublattice polarization analysis, supported by MESP and FMO analyses. We demonstrate that in the orthorhombic form, the proton tautomeric pathways are arranged in almost an antiparallel orientation, leading to a near-total cancellation of dipole moments internally. In contrast, the tetragonal form maintains a cooperative alignment necessary for robust FE. This study underscores that achieving a polar point group and low switching barriers are only the first steps in organic ferroelectric design. The spatial synchronization of molecular dipoles, dictated by crystal packing, remains the ultimate determinant of performance. By providing a deeper understanding of the interplay between electron density, hydrogen bonding, and symmetry, this work establishes a blueprint for leveraging polymorphism as a strategic tool for the rational design of next-generation organic electronic materials with tunable properties.

Author contributions

Conceptualization: P.M.; methodology: Y.G., S.D., B.M., A.J.E., V.K., and P.M.; investigation: Y.G., S.D., B.M., and P.M.; writing – original draft: Y.G. and P.M.; writing – review & editing: Y.G., B.M., and P.M.; funding acquisition: P.M.; resources: P.M.; supervision: P.M.

Conflicts of interest

There are no conflicts to declare.

Data availability

Experimental and theoretical data supporting this article have been included as part of the supplementary information (SI). Supplementary information: crystal structure details for the neutron data, comparison of the hydrogen bond distances obtained from neutron experiment, average neutron distances, and HAR model, molecular packing and ORTEP diagrams for the polymorphic forms, DRK plots, residual density maps, 2D and 3D static deformation density maps, topological parameters, representation of BCPs and associated bond paths, Laplacian maps, 2D electron localization function maps, comparison of experimental and optimized lattice parameters and their hydrogen bond geometries, pristine and modified structure of the polymorphic forms, optimized coordinates of the compound. See DOI: https://doi.org/10.1039/d6ma00724d.

CCDC 2541141, 2541177 and 2541178 contain the supplementary crystallographic data for this paper.^71a–c

Acknowledgements

Shiv Nadar Institution of Eminence Deemed to be University is thanked for providing fellowships to Y. G., S. D., and B.G., computational resources through MAGUS-02, and overall research support. P. M. thanks the ANRF (Grant No. CRG/2023/008038) for financial support. We acknowledge the support of the Australian Centre for Neutron Scattering, ANSTO, and the Australian Government for the neutron research infrastructure used in this work via proposal number 18582.

References

1. J. Bernstein, Oxford University Press Inc., New York, 2002.
2. S. Dutta, L. Negi and P. Munshi, Mater. Adv., 2024, 5, 7495–7515.
3. V. Jayant and M. Yusuf, Discov. Chem., 2024, 1, 52.
4. K. K. Jha, S. Dutta and P. Munshi, Cryst. Growth Des., 2018, 18, 1126–1135.
5. S. Dutta, A. Yadav, R. Boomishankar, A. Bala, V. Kumar, T. Chakraborty, S. Elizabeth and P. Munshi, Chem. Commun., 2019, 55, 9610–9613.
6. B. K. Saha, N. K. Nath and R. Thakuria, Chem. Rec., 2023, 23, e202200173.
7. P. Ambastha, V. Kushwaha, A. Magar, Y. Gupta, T. Parida, A. Kanjilal, R. Boomishankar and P. Munshi, NPG Asia Mater., 2026, 18, 4.
8. S. Dutta, Vikas, V. Thangavel and P. Munshi, ACS Appl. Electron. Mater., 2021, 3, 3633–3640.
9. H. Chung and Y. Diao, J. Mater. Chem. C, 2016, 4, 3915–3933.
10. C. M. Fernández-Posada, A. Castro, J.-M. Kiat, F. Porcher, O. Pena, M. Algueró and H. Amorín, Nat. Commun., 2016, 7, 12772.
11. M. E. Lines and A. M. Glass, Principles and applications of ferroelectrics and related materials, Oxford university press, 2001.
12. S. Horiuchi and Y. Tokura, Nat. Mater., 2008, 7, 357–366.
13. S. Horiuchi, Y. Tokunaga, G. Giovannetti, S. Picozzi, H. Itoh, R. Shimano, R. Kumai and Y. Tokura, Nature, 2010, 463, 789–792.
14. S. Horiuchi, K. Kobayashi, R. Kumai and S. Ishibashi, Nat. Commun., 2017, 8, 14426.
15. S. Horiuchi, S. Ishibashi and Y. Tokura, in Organic Ferroelectric Materials and Applications, ed. K. Asadi, Woodhead Publishing, 2022, pp. 47–84.
16. G. Gilli, F. Bellucci, V. Ferretti and V. Bertolasi, J. Am. Chem. Soc., 1989, 111, 1023–1028.
17. N. K. Hansen and P. Coppens, Acta Cryst. A, 1978, 34, 909–921.
18. P. R. Mallinson, G. T. Smith, C. C. Wilson, E. Grech and K. Wozniak, J. Am. Chem. Soc., 2003, 125, 4259–4270.
19. P. Munshi and T. N. Guru Row, Crystallogr. Rev., 2005, 11, 199–241.
20. P. Munshi, C. Jelsch, V. R. Hathwar and T. N. Guru Row, Cryst. Growth Des., 2010, 10, 1516–1526.
21. K. K. Jha, Y. Yadav, S. B. Srivastava, D. Chakraborty, P. Johari, S. P. Singh and P. Munshi, Cryst. Growth Des., 2019, 19, 3019–3029.
22. T. B. E. Grønbech, H. Kasai, J. Zhang, E. Nishibori and B. B. Iversen, Adv. Funct. Mater., 2024, 34, 2401703.
23. K. Tolborg and B. B. Iversen, Chem. – Eur. J., 2019, 25, 15010–15029.
24. P. Munshi and T. G. Row, CrystEngComm, 2005, 7, 608–611.
25. A. Nassour, M. Kubicki, J. Wright, T. Borowiak, G. Dutkiewicz, C. Lecomte and C. Jelsch, Acta Crystallogr., Sect. B: Struct. Sci., 2014, 70, 197–211.
26. V. V. Zhurov and A. A. Pinkerton, Z. Anorg. Allg. Chem., 2013, 639, 1969–1978.
27. G. Anil Kumar and V. R. Hathwar, Acta Crystallogr., Sect. B: Struct. Sci., 2023, 79, 450–461.
28. A. Genoni, L. Bučinský, N. Claiser, J. Contreras-García, B. Dittrich, P. M. Dominiak, E. Espinosa, C. Gatti, P. Giannozzi, J.-M. Gillet, D. Jayatilaka, P. Macchi, A. Ø. Madsen, L. Massa, C. F. Matta, K. M. Merz Jr, P. N. H. Nakashima, H. Ott, U. Ryde, K. Schwarz, M. Sierka and S. Grabowsky, Chem. – Eur. J., 2018, 24, 10881–10905.
29. J. Ciriaco-Pinheiro, R. Kondo, E. Aquino, M. Santos, M. Farias, S. Alves, F. Gil, A. Figueiredo, J. Lobato, R. Ferreira and O. Teu-Filho, inPiezoelectric Materials, ed. T. Ogawa, IntechOpen, London, 2016.
30. M. A. Khainovsky, E. V. Boldyreva and V. G. Tsirelson, Acta Crystallogr., Sect. B:Struct. Sci., Cryst. Eng. Mater., 2024, 80, 51–63.
31. Y. Cai, S. Luo, Z. Zhu and H. Gu, J. Chem. Phys., 2013, 139, 044702.
32. L. A. Malaspina, A. Genoni, D. Jayatilaka, M. J. Turner, K. Sugimoto, E. Nishibori and S. Grabowsky, J. Appl. Crystallogr., 2021, 54, 718–729.
33. B. Gruza, M. L. Chodkiewicz, J. Krzeszczakowska and P. M. Dominiak, Acta Cryst. A, 2020, 76, 92–109.
34. M. S. Pavan, S. Sarkar and T. N. G. Row, Acta Crystallogr., Sect. B:Struct. Sci., Cryst. Eng. Mater., 2017, 73, 626–633.
35. P. Munshi, E. Cameron, T. N. G. Row, J. D. Ferrara and T. S. Cameron, J. Phys. Chem. A, 2007, 111, 7888–7897.
36. Sakshi, Y. Gupta, A. Roy Choudhury and P. Munshi, Cryst. Growth Des., 2025, 25, 5395–5405.
37. S. Mukhopadhyay, M. J. Gutmann, M. Jiménez-Ruiz, D. B. Jochym, K. T. Wikfeldt, K. Refson and F. Fernandez-Alonso, Phys. Chem. Chem. Phys., 2017, 19, 32216–32225.
38. S. Horiuchi, R. Kumai and S. Ishibashi, Chem. Sci., 2018, 9, 425–432.
39. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
40. J. Cosier and A. M. Glazer, J. Appl. Crystallogr., 1986, 19, 105–107.
41. R. Piltz, J. Appl. Crystallogr., 2018, 51, 963–965.
42. P. W. Betteridge, J. R. Carruthers, R. I. Cooper, K. Prout and D. J. Watkin, J. Appl. Crystallogr., 2003, 36, 1487.
43. Bruker ( 2012), APEX5; SAINT; SADABS, Bruker AXS Inc., Madison, Winconsin, USA.
44. G. Sheldrick, Acta Cryst. A, 2015, 71, 3–8.
45. O. V. Dolomanov, L. J. Bourhis, R. J. Gildea, J. A. K. Howard and H. Puschmann, J. Appl. Crystallogr., 2009, 42, 339–341.
46. S. C. Capelli, H.-B. Burgi, B. Dittrich, S. Grabowsky and D. Jayatilaka, IUCrJ, 2014, 1, 361–379.
47. A. Volkov, P. Macchi, L. J. Farrugia, C. Gatti, P. Mallinson, T. Richter and T. Koritsanszky, University at Buffalo, State University of New York, NY, USA, 2006.
48. F. Kleemiss, O. V. Dolomanov, M. Bodensteiner, N. Peyerimhoff, L. Midgley, L. J. Bourhis, A. Genoni, L. A. Malaspina, D. Jayatilaka, J. L. Spencer, F. White, B. Grundkötter-Stock, S. Steinhauer, D. Lentz, H. Puschmann and S. Grabowsky, Chem. Sci., 2021, 12, 1675–1692.
49. F. Neese, WIREs Comput. Mol. Sci., 2022, 12, e1606.
50. R. H. Blessing, Acta Crystallogr., Sect. B:Struct. Sci., 1995, 51, 816–823.
51. A. O. Madsen, in Modern Charge-Density Analysis, ed. C. Gatti and P. Macchi, Springer, Netherlands, Dordrecht, 2012, pp. 133–163.
52. P. R. Mallinson, T. Koritsanszky, E. Elkaim, N. Li and P. Coppens, Acta Cryst. A, 1988, 44, 336–343.
53. C. B. Hübschle and B. Dittrich, Appl. Crystallogr., 2011, 44, 238–240.
54. G. A. Andrienko, See https://www.chemcraftprog.com, 2010.
55. M. E. Frisch, G. Trucks, H. B. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, G. Petersson and H. Nakatsuji, 2016.
56. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga and A. Cossard, J. Chem. Theory Comput., 2022, 19, 6891–6932.
57. A. Ivanauskaite, R. Lygaitis, S. Raisys, K. Kazlauskas, G. Kreiza, D. Volyniuk, D. Gudeika, S. Jursenas and J. V. Grazulevicius, Phys. Chem. Chem. Phys., 2017, 19, 16737–16748.
58. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
59. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
60. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
61. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953.
62. M. F. Hoedl, A. Chesnokov, D. Gryaznov, R. Merkle, E. A. Kotomin and J. Maier, J. Mater. Chem. A, 2023, 11, 6336–6348.
63. R. King-Smith and D. Vanderbilt, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 1651.
64. R. Resta, Rev. Mod. Phys., 1994, 66, 899.
65. F. H. Allen and I. J. Bruno, Acta Crystallogr., Sect. B:Struct. Sci., 2010, 66, 380–386.
66. M. Ahmed, A. Nassour, S. Noureen, C. Lecomte and C. Jelsch, Acta Crystallogr., Sect. B:Struct. Sci., Cryst. Eng. Mater., 2016, 72, 75–86.
67. R. F. Bader, Chem. Rev., 1991, 91, 893–928.
68. C. Kalaiarasi, M. S. Pavan and P. Kumaradhas, Acta Crystallogr., Sect. B:Struct. Sci., Cryst. Eng. Mater., 2016, 72, 775–786.
69. A. D. Stephen, R. Thomas, P. Srinivasan, V. Narayayanasamy and P. Kumaradhas, J. Mol. Struct., 2011, 989, 122–130.
70. B. Silvi and H. Ratajczak, Phys. Chem. Chem. Phys., 2016, 18, 27442–27449.
71. (a) CCDC 2541141: Experimental Crystal Structure Determination, 2026 DOI:10.5517/ccdc.csd.cc2r989n; (b) CCDC 2541177: Experimental Crystal Structure Determination, 2026 DOI:10.5517/ccdc.csd.cc2r99gv; (c) CCDC 2541178: Experimental Crystal Structure Determination, 2026 DOI:10.5517/ccdc.csd.cc2r99hw.

---