Fundamentally intertwined: anharmonic intermolecular interactions dictate both thermal expansion and terahertz lattice dynamics in molecular crystals

Abstract

We investigate the anisotropic thermal expansion behavior of a co-crystalline system composed of 4,4′-azopyridine and trimesic acid (TMA-azo). Using variable-temperature single-crystal X-ray diffraction (SC-XRD), low-frequency Raman spectroscopy, and terahertz time-domain spectroscopy (THz-TDS), we observe significant temperature-induced shifting and broadening of the vibrational absorption features, indicating changes in the intermolecular potential. Our findings reveal that thermal expansion is driven by anharmonic interactions and the potential energy topography, rather than increased molecular dynamics. Density functional theory (DFT) simulations support these results, highlighting significant softening of the potential energy surface (PES) with temperature. This comprehensive approach offers valuable insights into the relationship between structural dynamics and thermal properties, providing a robust framework for designing materials with tailored thermal expansion characteristics.

---

Solid-state dynamics are increasingly being found to be responsible for a wide range of material properties. For example, lattice dynamics play a direct and critical role in phenomena such as mass transport in porous materials,^1,2 charge carrier mobility in organic semiconductors,^3,4 and solid-state polymerization reactions.⁵ On the other hand, dynamics are indirectly related to properties such as polymorphism, as thermodynamics in solids are related to the energetic distribution of vibrational states – instead of the actual atomic motions themselves. Dynamic processes, and other illustrative examples, are often governed by crystalline phonon modes, which correspond to low-frequency (terahertz) vibrational motions.⁶ These phonon modes are directly related to both atomic-level interactions and bulk physical properties. Consequently, understanding the origins and behavior of solid-state dynamics is essential for the rational design of materials tailored for specific applications and processes.

One area where understanding solid-state dynamics is particularly crucial is in the control of thermal expansion behavior.^7,8 The design of materials with controllable thermal expansion properties is a long sought-after goal in materials science and crystal engineering, with applications in thermomechanical actuators, sensors, and aerospace engineering.⁹ The origins of thermally-induced structural changes lie in the nature of the intermolecular PES. Molecular crystals, in particular, often exhibit significant structural changes near ambient conditions, as they are often held together by weak, non-covalent forces.¹⁰ The potential energy curves of these weak interactions are often quite soft, which results in the generation of low-frequency (terahertz) vibrational dynamics – creating vibrational states that can be significantly populated by sub-ambient temperatures, e.g., a vibrational mode occurring at 1 THz (33.3 cm⁻¹) has a 7% population of the ν = 5 state at 100 K.

Despite the importance of the intermolecular PES, most studies on thermal expansion rely on static diffraction experiments.^11,12 While these methods provide valuable information on how structural modifications manifest with temperature, they do not directly offer insights into the various intermolecular potentials that drive thermal expansion. Instead, structural studies often attempt to deduce the role of dynamic effects indirectly, and recently there are numerous studies that have attributed thermal expansion phenomena to molecular dynamics.^13–17 However, increasing dynamics with rising temperature alone does not necessarily correlate with thermal expansion. Structural changes require asymmetric (anharmonic) interactions, as harmonic interactions will not produce any thermal expansion effects, nor change in vibrational frequency as a function of temperature.¹⁰ This distinction is crucial – while harmonic oscillators certainly exhibit larger-amplitude displacements with increasing temperature (and the associated increase in population of excited vibrational states), increased dynamics alone cannot account for thermal expansion. Rather, thermal expansion is derived from asymmetry in the intermolecular potential energy surface, which results in a change in the average position of the particles as a function of energy (or temperature).¹⁰

To highlight this distinction, we have performed a comprehensive investigation into a recently discovered material that exhibits interesting anisotropic thermal expansion behavior that has been suggested to be due to molecular dynamics. Co-crystals consisting of a 2 : 1 ratio of 4,4′-azopyridine (azo) to trimesic acid (TMA) have been previously studied and the thermal expansion behavior characterized using variable-temperature SC-XRD experiments from 190–290 K.¹³ This material contains two symmetry-independent azo molecules (and one TMA molecule) in the crystalline unit cell that form hydrogen-bonded sheets. The crystal was reported to exhibit interesting and distinct dynamic behavior as a function of temperature on the basis of static structural studies. Specifically, the azo molecules have been reported to undergo a ‘pedal’-type motion, which manifests as dynamic disorder, with the disorder increasing as temperature is increased. In this work, we have crystallized TMA-azo following the previous reported procedure (see ESI†). Variable-temperature SC-XRD experiments were acquired on a Rigaku Synergy (see ESI†), and we have expanded the range of temperatures studied to cover the 100–190 K region, and we observe the formation of a fully-ordered structure below 170 K. Interestingly, the two symmetry-independent azo molecules exhibit different onset temperatures of disorder, with one molecule showing disorder beginning at 170 K, while the second molecule has an onset temperature of 230 K.

The SC-XRD experiments provide clear indication of dynamic motion within the TMA-azo co-crystal, but in order to understand the relationship between dynamics, structure, and the underlying PES, low-frequency vibrational spectroscopy experiments were performed. Low-frequency (terahertz) vibrational spectroscopy directly provides insight into the intermolecular interactions that are responsible for thermal expansion, as the weak non-covalent interactions that are most influenced by temperature result in vibrational modes that occur in the terahertz frequency range.^6,10 This is in contrast to mid-IR vibrations, which arise from stronger, often covalent, interactions, which produce higher-frequency vibrational modes. The terahertz region is ideal for the study of molecular crystals, as there are often many terahertz vibrational modes (also referred to as ‘phonons’ in crystalline materials) due to the large number of unique intermolecular interactions in such systems. Accordingly, terahertz motions are often complex, involving dynamics of entire molecules, or large portions of molecules, which efficiently map out a large region of the potential energy surface.

The experimental low-frequency Raman spectra are shown in Fig. 1, and the THz-TDS spectra are shown in Fig. 2 (experimental methods provided in the ESI†). Both datasets show a considerable amount of temperature-induced spectral changes, with nearly all of the absorption features shifting to lower-frequencies and significantly broadening as temperature is increased. In order to assign the experimental spectral features to particular vibrational modes, solid-state density functional theory (DFT) simulations were performed using the Crystal23 software package¹⁸ using the 6-311G(2d,2p) basis set,¹⁹ and the PBE-D3 functional.^20,21 Both the Raman and IR simulations are in excellent agreement with the experimental low-temperature spectra.

Fig. 1 Low-frequency Raman spectra for TMA-azo from 90 K to 300 K (offset by 0.05), with the DFT-simulated spectrum shown in black.

Fig. 2 THz-TDS spectra for TMA-azo from 80 K to 300 K (offset by 0.5), with the DFT-simulated spectrum shown in black.

Many of the vibrational mode types in the terahertz region involve torsional motions of the azo group, with a simultaneous displacement of the TMA carboxylic acid group. However, the region between ca. 100–150 cm⁻¹, in particular, contains vibrational modes that sample the coordinates most directly related to the dynamic disorder observed in the static SC-XRD experiments. Two illustrative examples are shown in Fig. 3 – an IR-active mode predicted at 121.09 cm⁻¹, and the Raman-active mode predicted at 125.69 cm⁻¹. In Fig. 3a, a mixed rotational and torsional mode involving the pyridine groups on the azo molecules is coupled to a torsion of the TMA carboxylic acid groups. In Fig. 3b, this Raman-active mode maps out a coordinate involving the antisymmetric ‘pedal’-type motion of the two nitrogens in the azo bridge, coupled with a simultaneous torsion of the pyridine rings and, again, a torsion of the TMA carboxylic acids.

Fig. 3 Eigenvector displacement (red arrows) for two vibrational modes, an IR-active mode (a) and Raman-active mode (b), that are related to the dynamic disorder present in TMA-azo.

The experimental features that correspond to these vibrations, around 120 cm⁻¹ and 125 cm⁻¹ in the Raman and THz-TDS spectra, respectively, exhibit significant shifting and broadening with increasing temperature. As an illustrative example, the temperature-dependence of the Raman absorption feature is analyzed in Fig. 4, as well as the total occupancy of the two minor disordered sites, which provides a holistic quantification of the amount of thermally-induced dynamics in this system. The onset temperature of the disorder from variable-temperature SC-XRD experiments for the two symmetry-independent co-crystallized azo groups are indicated by blue and red dashed lines, respectively. Interestingly, these data show similar trends between the various experimental measurements. The frequency of the vibrational mode shifts from 123 cm⁻¹ to 115 cm⁻¹ – a significant 6.5% change from 80 K to 300 K. However, the temperature-dependence of the vibrational frequency is not linear, and three distinct regions with differing behavior are observed, below 140 K, between 140 K and 220 K, and above 220 K. This is also observed in the peak broadening, which shows similar trends. Importantly, the frequency of the vibration and the peak width provide insight into related – but distinct – physical properties. The peak width increases nearly 150% between 80 K and 300 K – which likely arises from a combination of a decrease in vibrational lifetime,²² and a smearing of the transition frequencies due to excitations originating from many states. The decrease in vibrational frequency with increasing temperature implies a significant softening of the vibrational potential as temperature is increased.¹⁰ This in and of itself is a hallmark of vibrational anharmonicity, and is relatively common for low-frequency vibrational modes.

Fig. 4 The total disordered site occupancy from the SC-XRD experiments (left), the center frequency (middle) and full-width at half-maxima (right) for the Raman absorption feature around 120 cm⁻¹. The blue and red vertical dashed lines represent the onset temperature of dynamic disorder for the respective azobenzne molecules.

Vibrational anharmonicity of terahertz phonons implies that the intermolecular potential – the same potential related to thermal expansion effects – is correspondingly anharmonic. But this naturally begs the question – do dynamics induce thermal expansion, or are structural changes pre-defined as the multidimensional PES is a static quantity? The PES certainly varies with volume, but dynamics are not required to generate a volume-dependent PES – implying that instead of dynamic-induced expansion, both vibrational dynamics and structural changes with temperature are simply outcomes of the underlying nature of the PES.

Here, we argue that it is not appropriate to simply imply that increased dynamics result in robust thermal expansion behavior. To illustrate this, it is important to consider the harmonic oscillator model, where increasing occupation of excited vibrational states leads to increased dynamics. This is evident from harmonic oscillator wavefunctions and their associated probability densities, which become more delocalized as the vibrational quantum number increases. However, in this model, structural changes are non-existent, despite the increase in dynamics. While in reality, most materials are not perfect harmonic oscillators, this helps to place into focus that linking increased dynamics directly to thermal expansion is overly simplistic. Instead, it is crucial to consider the multi-dimensional PES, which provides a comprehensive understanding of both dynamic and thermal expansion phenomena without solely relying on the nature of the dynamics.

This can be illustrated by Fig. 5, where the vibrational potential for the Raman-active mode predicted at 125.69 cm⁻¹ is plotted as a function of volume (see ESI†). The vibrational potential undergoes a significant softening as temperature is increased, ultimately leading to an increasingly-shallow potential where large-amplitude vibrational dynamics would arise. These vibrational potentials have similar behavior to the observed experimental properties described above. The three smallest-volume potentials show a more harmonic behavior, with only a slight softening as volume is increased. The fourth potential, corresponding to a volume of 2645 ų, is markedly different, and the presence of a second local minima emerges. Finally, in the most expanded two volumes, the barrier between the global and local minima vanishes, and an extremely shallow potential is observed. This nicely correlates with the observed structural and dynamic behavior.

Fig. 5 Vibrational potentials for the vibrational mode assigned to the 120 cm⁻¹ absorption feature as a function of unit cell volume. The normal mode coordinate refers to the scalar applied to the eigenvector predicted by the unconstrained optimization and frequency calculations.

This indicates that it is worth considering the true nature of so-called “dynamic-induced thermal expansion”. We suggest that instead of relating thermal expansion directly to increased dynamics, it is more appropriate to consider the multi-dimensional PES – Fig. 5 shows how this is achieved – the vibrational potential energy can be defined at a continuous series of volumes (or lattice parameters, intermolecular distances, and so on). Thus, when quanta of thermal energy are absorbed by a solid-state system, the free energy landscape is altered, and a new section of the thermally-defined PES is accessed. As a result of this new potential, both the vibrational dynamics and the thermal expansion change in tandem – as both are intimately linked, not necessarily to each other, but to the underlying potential.

Overall, we suggest that the potential energy picture is the most robust and complete way to describe temperature-induced structural and dynamic changes in materials. By approaching thermal expansion and related effects, such as negative thermal expansion, in molecular crystals through this view point, a more complete understanding of the link between supramolecular design, lattice dynamics, and bulk physical properties can be obtained. Through a combination of static structural studies, dynamic spectroscopic investigations, and consideration of the PES, a complete picture emerges that is able to describe the thermal properties of crystalline solids.

NK, JLH, WBS, WWB, SZ, PS: investigation, formal analysis, writing review and editing; LC, TMK: supervision, validation, writing review and editing; MTR: conceptualization, formal analysis, investigation, supervision, writing review and editing.

MTR thanks the NSF, award numbers CHE-2055402, DMR-204648, CHE-1725028, and CHE-2050793.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Notes and references

1. S. Krause, et al. , Angew. Chem., Int. Ed., 2020, 59, 15325–15341.
2. R. A. Klein, et al. , J. Am. Chem. Soc., 2023, 145, 2195521965.
3. P. A. Banks, et al. , Chem. Commun., 2022, 58, 12803–12806.
4. T. Okamoto, et al. , J. Am. Chem. Soc., 2020, 142, 9083–9096.
5. K. Hema, et al. , Chem. Soc. Rev., 2021, 50, 4062–4099.
6. P. A. Banks, et al. , Nat. Rev. Chem., 2023, 7, 480–495.
7. D. Das, et al. , Nat. Mater., 2009, 9, 36–39.
8. H.-Y. Ko, et al. , Phys. Rev. Mater., 2018, 2, 055603.
9. J. Chen, et al. , Chem. Soc. Rev., 2015, 44, 3522–3567.
10. M. T. Ruggiero, et al. , Chem. Commun., 2017, 53, 3781–3784.
11. A. van der Lee and D. G. Dumitrescu, Chem. Sci., 2021, 12, 8537–8547.
12. Z. Liu, et al. , Chem. Commun., 2018, 54, 5164–5176.
13. N. Juneja, et al. , Chem. Mater., 2023, 35, 7292–7300.
14. S. Scherb, et al. , Commun. Mater., 2020, 1, 8.
15. S. Li, et al. , Chem. Mater., 2023, 35, 2421–2428.
16. N. Juneja, et al. , Angew. Chem., Int. Ed., 2022, 61, e202202708.
17. N. Juneja, et al. , Cryst. Growth Des., 2023, 23, 524–531.
18. A. Erba, et al. , J. Chem. Theory Comput., 2022, 19, 6891–6932.
19. R. Krishnan, et al. , J. Chem. Phys., 1980, 72, 650–654.
20. J. P. Perdew, et al. , Phys. Rev. Lett., 1996, 77, 3865–3868.
21. S. Grimme, et al. , J. Chem. Phys., 2010, 132, 154104.
22. A. Fendt, et al. , Chem. Phys., 1981, 57, 55–64.

---

Footnotes

† Electronic supplementary information (ESI) available: Experimental methods, peak fitting analysis, experimental spectra, CIFs. CCDC 2377598–2377602. For ESI and crystallographic data in CIF or other electronic format see DOI: https://doi.org/10.1039/d4cc03307h

‡ These authors contributed equally to this work.

---