Mechanochemical cocrystallisation in a simplified mechanical model: decoupling kinetics and mechanisms using THz-TDS

Abstract

Mechanical compaction plays a dual role in mechanochemical cocrystallisation, simultaneously influencing reaction rates and crystallisation pathways. Here, we introduce a simplified single-punch compaction model that allows systematic variation of force (5–40 kN) and dwell time under controlled conditions. Using terahertz time-domain spectroscopy (THz-TDS) combined with a dual kinetic fitting strategy (free-fit vs. fixed-n Avrami models), we decouple kinetic rate constants from mechanism changes and map force-dependent transitions in a TPMA–PE cocrystallisation system. Increasing compaction force from 5 kN to 40 kN reduced the fitted rate constant k_free from 0.2147 to 0.1195 h⁻ⁿ, while increasing the Avrami exponent n_free from 0.6409 to 1.2057. This suggests a force-driven transition from diffusion-limited or heterogeneous nucleation to more interface-controlled one-dimensional crystallisation. These findings provide new mechanistic insight into how mechanical energy inputs shape solid-state transformations, with implications for process optimisation in pharmaceutical manufacturing, especially in continuous production environments where compaction profiles vary dynamically.

---

1 Introduction

The chemical diversity inherent in drug molecules presents both opportunities and challenges within the pharmaceutical sector. These molecules manifest in diverse physical states, such as the crystalline form of theophylline and the amorphous state of amoxicillin.^1,2 Crystalline substances often exhibit enhanced stability, but their applicability for drug therapy is often limited by low aqueous solubility and chemical instability in clinical settings.³ The pervasive issues of solubility and stability in bioactive compounds constrain the therapeutic potential of drugs and complicate their manufacturing and storage, posing substantial hurdles in clinical and commercial domains.¹

Mechanochemical cocrystallisation has emerged as a sustainable and versatile approach to tailor the physicochemical properties of active pharmaceutical ingredients (APIs), complementing conventional strategies such as salt formation and polymorph control. Among the available analytical tools, terahertz time-domain spectroscopy (THz-TDS) offers distinct advantages for in situ monitoring of solid-state transformations, owing to its sensitivity to crystalline lattice vibrations and non-destructive tracking of phase evolution.⁴ While THz-TDS has been successfully applied to extract kinetic parameters in various mechanochemical systems, its potential for decoupling intrinsic reaction rates from mechanism changes under systematically varied mechanical conditions remains underexplored.

Most prior mechanochemical studies have used ball milling or uniaxial pressing without precise control over compaction force at fixed compaction times.⁵ This makes it challenging to isolate the role of force–time profiles in dictating crystallisation pathways, and the interplay between mechanical input, particle rearrangement, and nucleation–growth mechanisms has not been rigorously quantified. To address this gap, we employ a simplified single-punch compaction model capable of applying well-defined forces (5–40 kN), dwell and compaction times. By coupling this with a dual kinetic fitting approach (free-fit vs. fixed-n Avrami models) applied to THz-TDS data, we systematically map force–mechanism–rate relationships in the cocrystallisation of a TPMA–PE model system. This framework enables the decoupling of rate acceleration from mechanism transitions, providing mechanistic insights and actionable guidance for optimising pharmaceutical manufacturing, especially under dynamically varying compaction profiles in continuous production.

Beyond the mechanochemical framework outlined above, it is important to place this work within the broader context of pharmaceutical cocrystal technology. To address these challenges, alternative structural entities like cocrystals, salts and different polymorphic forms have been developed to enhance the bioavailability and stability of drug compounds.^6,7 Cocrystal technology, in particular, has emerged as a robust strategy to improve the physicochemical properties of active pharmaceutical ingredients (APIs).⁶ The choice of coformers for API cocrystals is influenced by factors including functional group types, physical forms, and molecular sizes.⁸ Guided by hydrogen bonding, cocrystals are a universal technique to manage molecular aggregation, affecting crystal stacking patterns.⁹ Hydrogen bonds are crucial in the intermolecular interactions during the cocrystal formation process,^2,10 facilitating stable cluster formation through these bonds and proton transfers.¹¹ In selecting subjects for cocrystal research, the focus of this study has been particularly on the theophylline and malonic acid cocrystal due to theophylline's widespread use as a bronchodilator and its clinical limitations related to solubility and bioavailability. Malonic acid, offering two binding sites for hydrogen bonding, serves as an ideal coformer for theophylline cocrystals, expected to significantly enhance these properties and thus improve clinical efficacy and patient compliance.

Pharmaceutical cocrystals represent an environmentally sustainable alternative to traditional solid forms, providing a method to reduce dosage requirements while enhancing the stability of active components. Mechanochemical milling, alongside techniques like solvent evaporation, is a promising method for the industrial production of cocrystals. Yet, the fundamental mechanisms driving these solid-state transformations under mechanical activation remain areas of active investigation.¹² Studies indicate that some milled samples continue to exhibit cocrystal growth at room temperature even after being removed from the container.¹³ Additionally, liquid-assisted milling generally produces faster crystallisation rates and higher purity compared to non-liquid methods.^14,15 The complexity of these processes is further highlighted by observations of product layer formation acting as barriers, the effects of extended milling leading to agglomeration, and the significant influence of temperature on resulting polymorphic forms.^13,16

Understanding these intricate solid-state reactions has been profoundly advanced by the development of real-time, in situ monitoring techniques. Pioneering studies utilising synchrotron X-ray diffraction provided the first direct visualisation of reaction pathways, kinetics, and transient intermediates within operating milling equipment.^17,18 These methodologies, alongside complementary techniques like Raman spectroscopy, continue to yield crucial mechanistic insights into mechanochemical processes.¹⁹ Within this context, THz-TDS offers a powerful, non-contact analytical approach operating between the microwave and infrared regions (typically 0.1 THz to 10 THz). This frequency range probes low-energy phenomena such as intermolecular interactions (e.g., hydrogen bonding dynamics) and lattice vibrations (phonons), making it particularly sensitive to the subtle structural changes accompanying cocrystallisation and polymorphism.²⁰ Consequently, THz-TDS has proven instrumental in monitoring cocrystal growth, differentiating polymorphs, and complementing computational studies,^{2,13,16,21–25} often surpassing traditional techniques like PXRD and FTIR in distinguishing between forms with similar structures.²⁶

Despite these advances in monitoring and the exploration of cocrystal technology, a detailed understanding of the kinetic processes governing cocrystal formation and the precise influence of varying mechanochemical conditions (like controlled compaction forces) remains underdeveloped. This research gap hinders the development of predictive models and standardised approaches needed to fully leverage cocrystal advantages. In this study, we aim to delve deeper into the kinetics of compaction-assisted cocrystallisation using THz-TDS, employing a simplified mechanical model to isolate the effects of discrete impact events. THz-TDS, as a non-destructive technique sensitive to molecular interactions and structural phase transformations,²⁷ enables us to monitor the reaction progress in situ. The findings are expected to offer a foundational framework for understanding how mechanical forces modulate reaction kinetics and mechanisms in solid-state transformations, thereby facilitating the rational design of pharmaceutical formulations and manufacturing processes.

2 Methods

2.1 Materials

Theophylline (anhydrous, >99%) and malonic acid (>99%) were purchased commercially (Fluorochem, Hadfield, UK). High density polyethylene (PE, Induchem, Volketswil, Switzerland) powder was used as a diluent. The theophylline and malonic acid powders received were individually ground gently using a mortar and pestle to reduce particle size heterogeneity. The ground powders were then sieved using a mechanical sieve shaker (EVS1, Endecotts Ltd, London, UK) to isolate the particle size fraction between 100 μm to 250 μm. Sieved powders were stored under ambient laboratory conditions prior to mixing.

2.2 Sample preparation

2.2.1 Neat powder preparation and storage (sample 1). Sieved theophylline and malonic acid powders were weighed and thoroughly mixed in a 1 : 1 equimolar ratio using a vortex mixer (Whirlimixer, Fisher Scientific, Loughborough, UK) at the highest speed setting for 1 min. Approximately 1 g of this neat mixture (referred to as TPMA) was placed in an open 25 mL glass beaker. This beaker was stored inside a sealed desiccator maintained at 75% relative humidity (RH) using a saturated aqueous solution of sodium chloride at ambient laboratory temperature (approximately 20 °C), to eliminate humidity as a variable influencing cocrystallisation kinetics. This sample, stored as a loose powder, constituted sample 1 (see Table 1). It underwent spontaneous cocrystallisation under these humid conditions without any applied mechanical force beyond initial mixing, following a similar principle to previous studies.⁴ The reactants and the final cocrystal products produced by mechanochemical grinding described by the previous literature were analysed by powder X-ray diffraction (Cu Kα). The diffractograms were compared to a laboratory-prepared TPMA cocrystal reference and literature patterns. Within the instrumental resolution, all analysed pellets matched the target cocrystal; no additional reflections attributable to the starting components or alternative polymorphs were observed. At designated time points for analysis, approximately 40 mg of powder was carefully sampled from the beaker for THz-TDS measurement preparation (see section 2.4) (Fig. 1).

Table 1 Sample composition, storage form, and compaction conditions. TPMA refers to the 1 : 1 equimolar mixture of theophylline and malonic acid. Sample 1 was stored as powder and only pelletised (20 kN, 120 s) immediately before THz measurement

Sample ID	TPMA (% w/w)	PE (% w/w)	Storage form	Compaction during storage
1	100	0	Powder	None
2	20	80	Pellet	Hydraulic press (20 kN, 120 s)
3	20	80	Pellet	Simulator (5 kN, 0.6 s)
4	20	80	Pellet	Simulator (10 kN, 0.6 s)
5	20	80	Pellet	Simulator (15 kN, 0.6 s)
6	20	80	Pellet	Simulator (20 kN, 0.6 s)
7	20	80	Pellet	Simulator (30 kN, 0.6 s)
8	20	80	Pellet	Simulator (40 kN, 0.6 s)

---

Fig. 1 Experimental design illustrating the preparation and storage conditions for the different cocrystallisation methods. Method 1 involves storing the neat powder mixture (sample 1) under controlled humidity. Methods 2 and 3 involve preparing diluted pellets (samples 2–8) via different compaction techniques before storage under the same humidity.

2.2.2 Compaction-assisted cocrystallisation preparation (samples 2–8). For samples intended for compaction-assisted cocrystallisation, the TPMA mixture was first diluted with PE. Sieved theophylline and malonic acid were mixed in a 1 : 1 equimolar ratio. A 40 mg portion of this TPMA mixture was then combined with 160 mg of PE powder in a vial and homogenised using a vortex mixer to yield a 200 mg total sample mass with 80% w/w PE. This diluted mixture was then immediately compacted into pellets using one of two methods:

• Hydraulic press (sample 2, method 2): the 200 mg diluted mixture was compacted using a manual hydraulic press (15T, Specac Ltd., Kent, UK) equipped with a 13 mm diameter flat-faced die set. A force of 20 kN was applied and held for 120 s (2 min). It is not possible to record the force–time profile for the hydraulic press. Instead, compaction was carried out consistently with the same ramp/dwell characteristics for a 120 s (2 min) compaction process, comprising of a gradual ∼2 s ramp to target force, followed by constant-force dwell and subsequent rapid unloading by opening the hydraulic valve of the press.

• Compaction simulator (samples 3–8, method 3): the 200 mg diluted mixture was compacted using an HB50 compaction simulator (Huxley Bertram Engineering, Cambridge, UK) fitted with a 13 mm diameter flat-faced die set. Compaction was performed over a total of 0.6 s period using an HB50 simulator with a trapezoidal force profile, comprising a linear ramp-up (0.2 s) and ramp-down (0.2 s), and a dwell period (0.2 s) at peak force. Target peak forces ranging from 5 kN to 40 kN (as detailed in Table 1).

Immediately after compaction, pellets (samples 2–8) were transferred to the same 75% RH desiccator used for sample 1 and stored for the duration of the experiment (illustrated in Fig. 2) to allow compaction-assisted cocrystallisation to proceed. Pellets were periodically removed for THz-TDS analysis.

Fig. 2 Illustration of sample forms (neat powder in beaker, diluted pellets) stored in the 75% RH desiccator.

2.3 Pellet dimension monitoring

All pellets were prepared using a 13 mm diameter die set. The initial thickness of the compacted pellets (samples 2–8) typically ranged from 1.6 mm to 1.8 mm varying slightly with compaction force. During storage under 75% RH, the pellets exhibited some physical changes, including gradual expansion (increase in thickness) and potential surface friability. To account for this, the thickness of each pellet (samples 2–8) was carefully measured using a digital micrometre (±1 μm, 0 mm to 25 mm, Sealey, UK) immediately before each THz-TDS analysis. Care was taken during handling to minimise particle detachment. The measured thickness at each time point was used to calculate optical properties from the THz spectra. Reference pellets composed solely of 160 mg of PE were also prepared using the hydraulic press (20 kN, 120 s) and used for spectral referencing (section 2.4).

2.4 THz-TDS measurement and data processing

Terahertz time-domain spectroscopy (THz-TDS) measurements were performed using a TeraView Spectra 4000 spectrometer operating in transmission mode. The sample chamber was purged with dry nitrogen gas during measurements, minimising spectral interference from atmospheric water vapour.

For sample 1 (stored as powder), approximately 40 mg of the TPMA mixture was sampled at each time point, mixed with 160 mg of PE powder (via vortex mixer), and immediately pressed into a 13 mm diameter pellet using the hydraulic press (20 kN, 120 s). This pellet was then measured. For samples 2–8 (stored as pellets), the pellets were removed from the desiccator at each time point, their thickness measured (section 2.3), and then measured by THz-TDS.

For each sample measurement, 100 waveforms were acquired and averaged. A reference measurement was taken using a pure PE pellet (160 mg, prepared as described above) under identical conditions. The time-domain sample (E_sam(t)) and reference (E_ref(t)) signals were Fourier transformed to obtain the frequency-domain spectra (E_sam(ω) and E_ref(ω)). The frequency-dependent absorption coefficient α(ω) and refractive index n(ω) of the sample material within the PE matrix were calculated using standard equations based on the complex transmission coefficient T(ω) = E_sam(ω)/E_ref(ω) and the measured sample thickness d at that time point.²⁸ The spectral analysis focused on the frequency range from 0.2 THz to 3.0 THz. Note that, for the calculation of percent cocrystal content, because the integrated peak intensity is normalised to that of a reference cocrystal sample, values exceeding 100% can occur when the monitored pellet exhibits a stronger 2.16 THz feature than the reference, for example due to higher long-range order or reduced scattering losses after compaction.

3 Results and discussion

3.1 Characterisation using THz-TDS

Terahertz spectroscopy is a useful tool in elucidating the intricate process of cocrystal formation between the API theophylline and the short-chain diacid coformer malonic acid. The potential for hydrogen bond formation between the carboxyl groups of malonic acid and either the amide nitrogen or carbonyl group of theophylline is a key driver of this interaction. Given that both theophylline and malonic acid are crystalline, their interaction to form a cocrystal hinges on the establishment of hydrogen bonds, which are particularly sensitive in the terahertz region (<3 THz).

As illustrated in the top panel of Fig. 3, distinct spectral features below 3 THz are observed for both pure compounds. The initial equimolar mixtures exhibit prominent peaks around 0.98 THz and 1.64 THz, attributed to theophylline, alongside a peak at 2.5 THz corresponding to malonic acid. These observations provide a clear spectral fingerprint of the individual components prior to cocrystallisation.

Fig. 3 Top panel: THz absorption spectra of TP (theophylline), MA (malonic acid), and the initial mixture (sample 1 at t = 0). Middle panel: Comparison of the THz spectra of the initial mixture and the TPMA cocrystal reference synthesised using mechanochemical grinding. Bottom panel: Comparison of the spectra of sample 1 after 263 h storage and the cocrystal reference.

The middle panel of Fig. 3 further explores this process by comparing the spectra of the initial mixture with that of the pure theophylline–malonic acid (TPMA) cocrystal reference. This reference was prepared by mechanochemical neat grinding using a milling cycle at 25 Hz for 20 min, and was subsequently stored in a desiccator at 75% RH for 3 days.

A distinct peak emerging around 2.16 THz serves as a definitive marker of cocrystal formation. This spectral feature allows for precise tracking of the cocrystallisation progress, showcasing the sensitivity of terahertz spectroscopy in monitoring such molecular interactions.

The final measurement of sample 1 (stored powder), as depicted in the bottom panel of Fig. 3, shows a remarkable overlap with the pure TPMA cocrystal reference spectrum after 263 h. This indicates that nearly 100% conversion was achieved primarily by storing the mixed powders under humid conditions. This result highlights the effectiveness of utilising a humid environment to facilitate spontaneous cocrystal formation in this system.

3.2 Evolution of absorption spectra over time

The evolution of the absorption spectra over time provides a detailed view of the cocrystal formation process. Fig. 4 (left panel) shows representative absorption spectra for sample 1 monitored from 0 h to 263 h, revealing significant changes indicative of the ongoing reaction. Over time, the peaks at 0.98 THz and 1.64 THz (theophylline) gradually diminish, while the new peak at 2.16 THz, characteristic of the TPMA cocrystal, emerges and grows. Values above 100% indicate that the sample's 2.16 THz peak intensity exceeds that of the reference cocrystal spectrum, and should be interpreted as relative rather than absolute percentages.

Fig. 4 Left panel: Evolution of absorption spectra for sample 1 showing reactant peaks diminishing and the cocrystal peak emerging over time. Right panel: Normalised peak heights at 0.98 THz (black triangles, reactants consumption proxy) and 2.16 THz (red triangles, cocrystal formation X_c) as a function of time for sample 1. Normalisation was performed using different reference points for the two characteristic peaks: the 0.98 THz reactant peak was normalised with respect to its values at t = 0 (100%) and t = 263 h (0%), while the 2.16 THz cocrystal peak was normalised relative to its values at t = 0 (0%) and the benchmark state from a fully converted theophylline–malonic acid cocrystal synthesised via the grinding method (100%). Error bars represent the standard deviation over two independent experiments.

Initial measurements for sample 1 were conducted at approximately 12 h intervals for the first 100 h, with subsequent measurements spaced every 24 h. The plot on the right side of Fig. 4 offers a quantitative comparison based on normalised peak heights, averaged over two independent experiments. The fractional conversion to cocrystal (X_c) was estimated by normalising the absorption coefficients of two characteristic peaks:

For the reactant peak at 0.98 THz, the absorption coefficient at 0 h was set to 100% and that at 263 h was set to 0%, representing complete reactant consumption.

For the cocrystal peak at 2.16 THz, the absorption coefficient at 0 h was taken as 0%, and the peak height of a reference sample of pure theophylline–malonic acid cocrystal synthesised via the neat grinding method was used as the 100% conversion benchmark. During the first 100 h, the decrease in the reactant peak occurs at a slightly faster rate than the emergence of the cocrystal peak. This suggests that while the reactants are being rapidly consumed or transformed from their initial crystalline state, the complete formation and growth of detectable cocrystal domains requires additional time. Mechanistically, the slightly faster decay of the 0.98 THz reactant marker relative to the build-up of the 2.16 THz cocrystal marker during the first ≈100 h suggests an early-stage reorganisation/disordering of reactant crystallites prior to the growth of long-range ordered cocrystal domains, which is consistent with the Avrami analysis in section 3.3.

As the process continues, the system gradually approaches an equilibrium or completion state, with near-complete conversion observed by the 150 h mark. The transition underscores the kinetic nature of the spontaneous cocrystallisation process under these conditions.

3.3 Kinetic analysis of spontaneous cocrystallisation (sample 1)

The kinetics of cocrystal formation for sample 1 were analysed using the Avrami–Eroféev equation. This model is widely used to describe the dynamics of solid-state phase transformations by integrating both nucleation and growth mechanisms. The model assumes nucleation occurs uniformly throughout the sample volume and that crystal growth proceeds at a constant rate (G) from these nuclei, often simplifying the growth geometry (e.g., spheres, rods). The Avrami model implicitly accounts for the impingement and overlap of growing domains, although the explicit geometry and spatial distribution are simplified in this application.

Under these assumptions, the fraction transformed (X_c) at time t can be described by:

X_c(t) = 1 − exp(−ktⁿ) (1)

The rate constant k (h⁻ⁿ) reflects the overall speed of the transformation, encompassing both nucleation and growth rates, while the Avrami exponent n provides insights into the dimensionality of growth and the nature of the nucleation process. Integer values of n are typically associated with specific idealised mechanisms (e.g., n = 1 for one-dimensional growth with instantaneous nucleation, n = 3 for three-dimensional growth with constant nucleation rate, n = 4 for three-dimensional growth with instantaneous nucleation).

In the initial 100 h for sample 1, the observed discrepancy between the rates of reactant peak decrease and cocrystal peak increase (Fig. 4, right panel) may suggest complexity beyond the simplest Avrami assumptions, potentially involving intermediate phases or different sensitivities of the THz peaks to initial reactant disordering versus cocrystal lattice formation. However, for quantifying the overall process, the Avrami equation was applied to the fractional cocrystal conversion (X_c) estimated from the THz peak height at 2.16 THz.

Eqn (1) can be rearranged for linear fitting:

ln[−ln(1 − X_c(t))] = ln k + n ln t (2)

where a plot of ln[−ln(1 − X_c)] versus ln t should yield a straight line with slope n and intercept ln k according to eqn (2). As shown in Fig. 5, fitting the averaged data for sample 1 (average of two runs) yielded an Avrami exponent of n = 0.97 ± 0.19 and a rate constant of k = 0.147 h⁻ⁿ (units depend on n), with a coefficient of determination of R² = 0.938. The fitted value of n was obtained from the linear regression slope, and its uncertainty corresponds to the half-width of the 95% confidence interval.

Fig. 5 Avrami plot (ln[−ln(1 − X_c)] vs. ln t) for the spontaneous cocrystallisation of sample 1 (average data). The red dashed line represents the linear fit using eqn (2) (free fit). The black dashed line represents the fit constrained to n = 1.

The Avrami exponent n being close to 1 suggests that the rate-limiting step might approximate one-dimensional growth or perhaps surface-controlled growth. This directionality might be influenced by the inherent crystal habits of the reactants or, potentially, by experimental conditions such as a humidity gradient within the powder bed, where moisture penetration primarily from the top could lead to anisotropic growth. One-dimensional growth is often associated with processes where crystal growth is confined linearly, such as along a crystallographic axis or an interface.

Given that n values are often interpreted based on integer values representing idealised mechanisms, and the fitted n is very close to unity, we also performed the analysis constraining n = 1. This yielded k = 0.140 h⁻¹ with R² = 0.937 (Table 2). The close agreement in R² values suggests that assuming n = 1 is a reasonable simplification for this system under spontaneous conditions.

Table 2 Avrami parameters from fitting sample 1 kinetic data (average, N = 3)

Fit type	n	k (h^−n)	R ^2
Free fit	0.9727 ± 0.1883	0.147	0.938
Fixed (n = 1)	1	0.140	0.937

---

3.4 Compaction-assisted cocrystallisation: hydraulic press (sample 2)

Building on the understanding of spontaneous cocrystallisation, the investigation progressed to compaction-assisted cases. Initial attempts using pure TPMA were hindered by strong THz absorption (exceeding instrument dynamic range above ca. 1.5 THz) and rapid conversion (<24 h). To mitigate these issues, an 80% w/w PE matrix was used for dilution. Sample 2 involved compacting the 20% TPMA/80% PE mixture with a hydraulic press (20 kN, 120 s) and storing the resulting pellet at 75% RH.

Despite dilution, the kinetic profile of sample 2 remained faster than sample 1 (Fig. 6, left panel). However, the presence of PE, while moderating the rate, also reduced reactant contact area and likely hindered mass transport, leading to incomplete conversion, reaching only approximately 80% by the end of the experiment (ca. 150 h). This highlights a trade-off between moderating kinetics for measurement and achieving complete conversion, suggesting a need to optimise parameters like diluent concentration and compaction conditions. Fitting the sample 2 data yielded significantly different Avrami parameters compared to sample 1 (Fig. 6, right panel; Table 3), notably a lower n value (0.45).

Fig. 6 Comparison of cocrystallisation kinetics between sample 1 (powder, black) and sample 2 (hydraulic press pellet, red). Left panel: Fractional cocrystal formation (X_c) versus time. Right panel: Corresponding Avrami plots (ln[−ln(1 − X_c)] vs. ln t) with free-fit lines.

Table 3 Avrami parameters calculated from the free-fit model for all samples, including R² values of the linear fit

Sample	Compaction force^a (kN)	n	log k	k (h^−n)	R ^2
Notes:a Force applied during initial pellet preparation for samples 2–8. Sample 1 was stored as powder.b Units of k depend on the fitted value of n.c Sample 4 (10 kN) shows a modest deviation relative to adjacent forces; it is retained in the dataset (Table 3) but excluded from the empirical force–kfree trend fit in Fig. 8 for clarity. Inclusion/exclusion does not change the monotonic trends or mechanistic interpretation.
1	Powder	0.9727 ± 0.1883	−1.9190	0.1468	0.938
2	20 (hydraulic)	0.4454 ± 0.1043	−1.1969	0.3021	0.889
3	5 (simulator)	0.6409 ± 0.1825	−1.5384	0.2147	0.891
4	10 (simulator)^c	1.3004 ± 0.2005	−2.4359	0.0875	0.965
5	15 (simulator)	0.8984 ± 0.3361	−1.7135	0.1802	0.851
6	20 (simulator)	1.0264 ± 0.2468	−2.0071	0.1344	0.945
7	30 (simulator)	1.1196 ± 0.3033	−2.0517	0.1285	0.947
8	40 (simulator)	1.2057 ± 0.3668	−2.1244	0.1195	0.896

---

The faster initial kinetics observed in the compacted sample 2 suggest that increased particle contact and reduced inter-particle distance, induced by compaction, enhance the molecular interactions required for cocrystal formation. Although inert polyethylene (PE) particles constitute 80% of the mixture by mass, the compaction process significantly reduces the overall inter-particle distance within the pellet. This densification promotes spatial reorganisation of all components, thereby increasing the likelihood of contact between reactive particles despite the high PE content. As a result, the absolute number of effective contacts between reactive species per unit volume is enhanced, facilitating nucleation and growth. Under these specific conditions (hydraulic press, prolonged duration), this beneficial effect of increased particle contact appears to outweigh potential negative impacts associated with reduced porosity or limited molecular mobility within the denser pellet structure.

Mechanistic interpretation. Relative to the spontaneous powder (sample 1; n ≈ 1.0, k ≈ 0.14 h⁻¹ to 0.15 h⁻¹ under fixed-n; Table 2), the hydraulic-press condition (20 kN, 120 s) exhibits a higher apparent rate (k_free ≈ 0.302 h⁻ⁿ) and a sub-unity value for n (≈0.45), consistent with densification-enhanced contact driving faster transformation while the mechanism reflects transport or interface limitations (Fig. 6; Table 3). Taken together, the faster yet sub-unity-n behaviour under prolonged compression resonates with PXRD/MCR-ALS observations on TP–MA where humidity and transient phases governed pathway selection during milling.²⁹ Here, by prescribing the force–time profile, we resolve why an apparent k_free increase can coexist with transport/interface limitations (low n) under densification-driven contact generation. To further elucidate how mechanical force influences the cocrystallisation mechanism across a wider compaction range, we systematically varied the peak force using a compaction simulator.

3.5 Compaction-assisted cocrystallisation: compaction simulator (samples 3–8)

The impact of varying compaction force was systematically examined using an HB50 compaction simulator, applying forces from 5 kN to 40 kN over a brief 0.6 s period to diluted TPMA/PE mixtures (samples 3–8). The resulting relative pellet thickness decreased non-linearly with increasing force (Fig. 7), indicating initial powder rearrangement and particle contact improvement at low forces (5 kN to 20 kN), followed by particle deformation/fracture and reduced compressibility at higher forces (20 kN to 40 kN), where thickness approached a plateau around 0.28 mm. This behaviour is consistent with typical powder compaction models like the Heckel equation³⁰ or exponential decay models. Fitting the thickness data to y = ae^bx + c yielded R² = 0.997, confirming the diminishing returns of compression at higher forces under these conditions.

Fig. 7 Relative pellet thickness (sample pellet thickness minus that of a pure PE reference pellet, both measured shortly after compaction) versus applied compaction force for samples 3 to 8 (200 mg total mass, 80% PE). Shaded regions distinguish approximate low-force (5 kN to 15 kN) and high-force (20 kN to 40 kN) regimes discussed in the text.

Kinetic analysis was performed on samples 3–8 using the free-fit Avrami model (Table 3). The 10 kN kinetic profile (sample 4) shows a modest deviation relative to neighbouring forces; however, it lies within the experimental scatter typical of in situ mechanochemical monitoring across THz/Raman/PXRD modalities.^31,32 Crucially, the global trends (n_free↑ with force, k_free↓ with force; k_fixed↑ under n = 1) and the mechanistic assignment are unchanged whether sample 4 is included in or excluded from the trend fit (Fig. 8; Table 4). For clarity of regression display, sample 4 is omitted from the empirical force–k_free fit in Fig. 8, but it is retained in the dataset summary (Table 3).

Fig. 8 Relationship between compaction force and the rate constant k derived from the free-fit Avrami model for samples 3 and 5–8. Sample 4 (10 kN) is omitted due to anomalous behaviour. The trend suggests decreasing kinetics at higher forces, potentially modellable by exponential decay (see text).

Table 4 Comparison of kinetic parameters derived from free-fit (n_free, k_free) and fixed-n (n_fixed = 1, k_fixed) Avrami model approaches for simulator-compacted samples. Data for 10 kN (sample 4) excluded. Units of k_free depend on n_free

Compaction force (kN)	n free	k free (h^−n)	R ^2free	n fixed	k fixed (h^−1)	R ^2fixed
5	0.6409	0.2147	0.8910	1	0.1165	0.6114
15	0.8984	0.1802	0.8510	1	0.1523	0.8400
20	1.0264	0.1344	0.9450	1	0.1402	0.9446
30	1.1196	0.1285	0.9470	1	0.1548	0.9366
40	1.2057	0.1195	0.8960	1	0.1680	0.8701

---

While the exact cause of this discrepancy remains unclear, possible factors include experimental variability in sample preparation or measurement, or a genuine transition in the cocrystallisation mechanism at this specific compaction force. Given its clear deviation from the general trend, sample 4 was considered an outlier and excluded from the analysis of the relationship between kinetic parameters and compaction force presented in Fig. 8.

Excluding the outlier sample 4, the rate constant k_free generally decreased as compaction force increased from 5 kN to 40 kN (Table 3, Fig. 8). Specifically, k_free decreased from 0.2147 h⁻ⁿ at 5 kN to 0.1195 h⁻ⁿ at 40 kN with slight variations observed among the intermediate samples.

To describe this trend, the k_free values of samples 3 and 5–8 were fitted to an exponential decay model of the form k = ae^bF + c, where F represents the applied compaction force. The fitting yielded parameters of a = 0.1539, b = −0.0597, and c = 0.1028, with a coefficient of determination of R² = 0.9316, indicating an excellent fit.

Mechanistic interpretation. Across 5 kN to 40 kN, the free-fit analysis shows n_free increasing from ∼0.64 to ∼1.21 (Fig. 9, left; Table 4), while k_free decreases from 0.2147 to 0.1195 h⁻ⁿ (Fig. 8; Table 4). Interpreted within the Avrami–Eroféev framework, the rise in n_free indicates a force-induced drift in the dominant mechanism (toward more interface-controlled growth), whereas the fall in k_free reflects the conflation of intrinsic rate with mechanism changes. To decouple these effects, a fixed-n analysis (n = 1) reveals that k_fixed increases from 0.1165 to 0.1680 h⁻¹ with force (Table 4), showing that compaction intrinsically accelerates the transformation even when the apparent (mechanism-convolved) k_free declines at high force. The short compaction (0.6 s) and higher densification at elevated forces plausibly reduce pore space and restrict mobility, consistent with pore-constrained growth,³³ thereby shifting the apparent mechanism while raising the intrinsic rate. Data at 10 kN (sample 4) exhibited atypical behaviour; as pre-specified, it is excluded from trend fitting but does not alter the mechanistic assignment or regression quality for the remaining forces (Fig. 9; Table 4).

Fig. 9 Comparison of fitted Avrami parameters for simulator-compacted samples (excluding 10 kN) under varying compaction forces. Left panel: Avrami exponent n obtained from free-fit analysis (n_free, black circles), with the fixed n = 1 assumption shown as a dashed red line. Right panel: Rate constant k values derived from free-fit (k_free, black inverted triangles) and fixed-n (k_fixed, red triangles) models.

The progressive rise of n_free with force (Fig. 9) also mirrors molecular-level predictions that impact geometry (e.g. incident angle) controls mixing efficiency and lattice disruption in organic molecular crystals.³⁴ Our macroscale mapping thus provides experimental context for those simulations at tablet-relevant length and time scales.

Overall, these results indicate moderate reaction rates under high humidity. The decreasing trend in k_free with increasing compaction force (above 5 kN or 15 kN) suggests that while initial compaction may improve contact, higher forces under these short (0.6 s) simulator conditions may primarily limit molecular mobility or create less favourable reaction environments compared to lower forces.

Comparing across methods, the simulator-compacted samples (3, 5–8) generally reacted slower than the hydraulically pressed sample 2 (k_free ≈ 0.30 h⁻ⁿ). However, in the low-force regime (5–15 kN), these samples reacted faster than the undiluted spontaneous powder sample 1 (k_free ≈ 0.15 h⁻ⁿ). This highlights the complex influence of both compaction force magnitude and duration. Interestingly, the k_free value for sample 6 (20 kN simulator, k ≈ 0.13 h⁻ⁿ) is similar to that of the spontaneous powder (sample 1). This suggests a complex balance between factors enhanced by compaction (contact area, reduced distance) and factors potentially hindered (molecular mobility, moisture transport), with mobility possibly playing a crucial role under these specific diluted and humid conditions. It is acknowledged that humidity can, in some mechanochemical systems, facilitate transformations by enhancing surface mobility, acting as a plasticiser, or participating in proton transfer processes.⁴¹ In the present study, humidity was held constant at ∼75% RH to ensure that the observed variations in Avrami parameters could be attributed solely to changes in compaction force and dwell time. Systematic variation of humidity, particularly at high RH (e.g., ≥90%), will be explored in future work to investigate potential water-force synergies.

Role of PE as an inert diluent. The inclusion of 80% w/w polyethylene (PE) in the formulation could, in principle, dilute reactive contact points between TPMA and the coformer. Our results indicate that this dilution does not confound the mechanistic trends observed. Specifically, kinetic analysis reveals a consistent increase in the Avrami exponent n_free (from ∼0.64 to ∼1.21) and a concurrent decrease in the apparent rate constant k_free with increasing compaction force. These opposing shifts point to a force-induced mechanism transition towards more interface-controlled crystallisation, rather than diminished reactivity from dilution.³³ This interpretation is reinforced by fixed-n fits (n = 1) , which show a clear increase in k_fixed with force (from 0.1165 to 0.1680 h⁻¹), revealing intrinsic rate enhancement despite the falling k_free. Additional support comes from literature reports of similar high inert-matrix loadings (70–85%) showing kinetics dominated by compaction and mixing dynamics.^31,35 The presence of PE also ensured reproducible tablet geometry, mitigated excessive local heating, and provided a non-reactive medium for in situ THz-TDS without signal saturation. Preliminary experiments with pure TPMA tablets exhibited very rapid crystallisation (<40 h) and markedly different mechanical behaviour (e.g., higher elastic modulus, lack of plastic deformation), making them unsuitable as a direct baseline for comparison.

3.6 Decoupling growth dimensionality and kinetic effects in compaction-driven cocrystallisation

As emphasised in a recent critical review, non-integer or sub-unity values of n often reflect evolving mechanisms, transport limitations or time-dependent nucleation; hence cross-condition comparisons of k alone are unsafe.³³ Our paired free-/fixed-n analysis operationalises this caution by quantifying mechanism drift while exposing the intrinsic rate trend.

The analysis of solid-state kinetics using the Avrami–Eroféev model requires careful interpretation, as the fitted parameters, the exponent n and the rate constant k, are inherently linked and reflect the combined effects of nucleation, growth dimensionality, and the intrinsic reaction rate. In our initial free-fit analysis (Table 3), significant variations in n were observed across different compaction forces (ranging from approximately 0.6 to 1.2, excluding the 10 kN outlier). According to the theoretical framework,^33,36 such variations, particularly non-integer or sub-unity values, strongly suggest that the underlying crystallisation mechanism is not constant but evolves with the applied mechanical stress. These changes could involve shifts in nucleation behaviour (e.g., from instantaneous to time-dependent) or alterations in the dominant growth dimensionality influenced by diffusion limitations or structural constraints within the compacted matrix.

The variability in crystallisation mechanism complicates the direct interpretation of the rate constant k obtained from the free-fit Avrami model (k_free). As shown in Fig. 9 (right panel, black inverted triangles), k_free exhibits a generally decreasing trend with increasing compaction force, ranging from 0.2147 h⁻ⁿ at 5 kN to 0.1195 h⁻ⁿ at 40 kN. These values are taken from Table 4 which presents the kinetic parameters specifically derived from the free-fit model for simulator-compacted samples.

It is important to note that Table 4 differs from Table 3. While Table 3 summarises the general free-fit results across all samples, Table 4 focuses on samples analysed for direct comparison between free-fit and fixed-n (n = 1) models. The k_free values discussed here correspond to the free-fit results presented in Table 4 and plotted in Fig. 9.

The observed decrease in k_free with increasing compaction force raises the question of whether compaction intrinsically slows the reaction kinetics, or whether this apparent trend is predominantly influenced by concurrent changes in the crystallisation mechanism, as reflected in the variations of n_free.

To disentangle these intertwined effects and specifically probe how compaction influences the intrinsic reaction rate under a hypothetically constant mechanism, we performed a secondary kinetic analysis. We constrained the Avrami exponent to a fixed value, n = 1, across all compaction conditions (samples 3, 5–8) and re-fitted the data to determine the corresponding rate constant (k_fixed). We chose n = 1 as it represents a common reference point, often associated with phase-boundary controlled reactions or one-dimensional growth, and it falls within the range of n_free values observed, particularly at higher forces. This fixed-n approach allows us to isolate the apparent effect of compaction force on the rate constant, assuming the reaction proceeds via a consistent pathway.

The results from this constrained fit revealed a markedly different kinetic trend for k_fixed compared to k_free (Fig. 9, right panel, red triangles; Table 4). Instead of decreasing, k_fixed initially increased with compaction force, rising from 0.117 h⁻¹ at 5 kN to a peak of 0.152 h⁻¹ at 15 kN, before plateauing and slightly increasing again at the highest forces (0.168 h⁻¹ at 40 kN).

The divergence between the k_free and k_fixed trends is highly informative. Firstly, it suggests that compaction can intrinsically accelerate the reaction rate; the initial increase in k_fixed (up to 15 kN) indicates that bringing particles into closer contact and increasing local density likely enhances the fundamental reaction rate assuming the mechanism remains constant. Secondly, and crucially, it implies that mechanistic changes significantly impact the observed kinetics. The fact that the observed rate constant k_free decreases across the force range, contrary to the initial trend of k_fixed, demonstrates that changes in the reaction mechanism (reflected in the varying n_free) play a dominant role in determining the overall kinetics. At lower forces, factors leading to n_free < 1 (e.g., diffusion control, complex nucleation) appear to limit the overall rate even if the intrinsic rate is increasing. Conversely, at higher forces, the development of severe structural constraints likely hinders molecular mobility and alters growth pathways (contributing to the varying n_free), leading to the observed lower values of k_free despite the intrinsic rate (k_fixed) potentially plateauing or slightly increasing.

Therefore, this comparative modelling approach allows us to conclude that mechanical compaction is not merely a kinetic tuner but a powerful mechanistic modulator. The overall observed rate (k_free) results from a complex interplay between the intrinsic acceleration due to increased contact/density and the potentially counteracting effects of evolving mechanistic pathways and physical constraints imposed by the compacted structure.

Implications for manufacturing and process analytics. First, in continuous direct compression (CDC), peak force and, to some extent, dwell time are tunable levers. Our force–mechanism–rate map enables selecting conditions that target desired pathways (e.g. interface-controlled growth, higher n) rather than maximising apparent k_free alone, aligning with CDC work that stresses control of force history and in-line analytics.³⁷ Second, THz-TDS is a complementary PAT to Raman/PXRD with heightened sensitivity to low-frequency lattice and hydrogen-bonding modes that fingerprint co-crystal formation and subtle structural changes; recent demonstrations include resolving guest-water orientation in MOFs and discriminating closely related crystalline frameworks.^38–40 Together, these points support the opportunity for mechanism-aware scale-up through joint control of peak force, dwell and humidity, with THz-TDS enabling non-destructive tracking of both kinetics and mechanism in real time.

4 Conclusion

Building upon the kinetic trends observed in the previous sections, we now contextualise these findings within the broader literature and mechanistic models of solid-state transformation.

This study integrates in situ THz-TDS with a dual Avrami fitting framework (free-fit vs. fixed-n) under programmable single-punch compaction to explicitly map force–mechanism–rate relationships in a model TP–MA cocrystallisation. This combination goes beyond prior in/ex situ diffraction-led monitoring and generic TRIS overviews by separating intrinsic rate changes from mechanism drift via paired fits, thereby enabling defensible interpretation of k across varying n. Methodologically, it complements time-resolved in situ monitoring advances in mechanochemistry³² and recent PXRD/MCR-ALS kinetic extraction on TP–MA,²⁹ while our force-programmed single-punch approach addresses the role of impact geometry and loading history highlighted by MD/indentation studies³⁴ and caveats in Avrami-based assignments.³³

Our findings demonstrated that compaction significantly enhanced reaction rates by increasing contact points between particles, particularly when compared to loose powder systems. Extended compaction duration further accelerated the process, while higher compaction forces, under the same force profile, paradoxically slowed the reaction due to limited molecular mobility. The addition of polyethylene (PE) reduced the effective contact area, slowing the reaction and limiting cocrystal formation upon completion.

THz-TDS was effective in monitoring distinct absorption peak features driven by specific molecular vibrations and transitions, providing direct evidence of molecular interactions and structural changes. This non-invasive technique enables the real-time optimisation of cocrystallisation conditions, offering valuable guidance for the development of stable and effective pharmaceutical formulations.

Our results confirmed that spontaneous cocrystallisation can occur without external mechanical force, consistent with findings reported by Davies et al. (2022).⁴ However, in our system, approximately 100 h were required to achieve complete conversion, compared to 50 h in previous studies, attributed largely to the use of larger particle sizes (100 μm to 250 μm versus 50 μm to 70 μm). This further highlights the role of particle size and the potential for compaction-assisted methods to enhance cocrystallisation kinetics.

Analysis of the spectral evolution indicated complex initial stages, potentially involving rapid reactant transformation preceding the slower bulk nucleation and growth phase modelled by the Avrami–Eroféev equation. The overall rate constant k, determined via Avrami analysis and ranging from approximately 0.1 h⁻¹ to 0.3 h⁻¹ (depending on conditions and model fit), indicated moderate crystallisation kinetics at room temperature. The dimensionality exponent n, approximately 1 for the spontaneous case, suggested a predominantly one-dimensional growth mechanism might be operative. While mechanochemical reactions primarily accelerate cocrystallisation by increasing contact area, our findings suggest that moisture mobility was the dominant driving factor in facilitating molecular rearrangement under the humid conditions studied. The presence of PE affected the kinetics and final conversion but did not inhibit cocrystallisation itself. Moreover, the observed force–mechanism relationships are robust under the 80% PE dilution level, with mechanical input dominating over any potential concentration effects.

Importantly, this work demonstrates that mechanical compaction significantly modulates both the rate and the mechanism of spontaneous cocrystallisation. Increasing compaction force induces a transition in the kinetic behaviour, reflecting the critical role of physical constraints in directing solid-state transformation pathways. Through a dual kinetic fitting approach comparing free and fixed-n models, we could separate rate-driven and mechanism-driven effects, offering a valuable framework for analysing complex crystallisation behaviours under mechanical stress.

Across 5–40 kN, n_free increases from ∼0.64 to ∼1.21 while k_free decreases from 0.2147 to 0.1195 h⁻ⁿ, whereas under a fixed-n assumption (n = 1) k_fixed increases from 0.1165 to 0.1680 h⁻¹. This divergence clarifies that compaction concurrently accelerates the intrinsic transformation and shifts the apparent mechanism—precisely the rationale for pairing free- and fixed-n analyses.

These insights have important implications for designing mechanically assisted manufacturing processes, where precise control over crystallisation dynamics is essential. By adjusting factors such as diluent concentration, humidity levels, and compaction force/duration, future studies can systematically explore optimal conditions for achieving high yields, desired polymorphs, and complete conversion in cocrystal products, providing an intriguing avenue for further research.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study are openly available in the University of Cambridge's Apollo repository at https://doi.org/10.17863/CAM.119087.

Acknowledgements

The authors gratefully acknowledge the use of the terahertz time-domain spectroscopy (THz-TDS) facility at the Department of Chemical Engineering, University of Cambridge. We would also like to thank J. Lee, M. Ma, C. Leung, and other members of the Terahertz Applications Group for valuable discussions and support.

Notes and references

1. A. V. Trask, W. D. S. Motherwell and W. Jones, Int. J. Pharm., 2006, 320, 114–123.
2. L. Yu, Adv. Drug Delivery Rev., 2001, 48, 27–42.
3. M. I. Jarvis, V. Krishnan and S. S. Mitragotri, Bioeng. Transl. Med., 2018, 4, 5–16.
4. R. D. Davies, N. J. Vigilante, A. D. Frederick, V. S. Mandala and M. A. Mehta, J. Chem. Crystallogr., 2022, 52, 479–484.
5. T. Friščić, I. Halasz and P. A. Julien, Chem. Soc. Rev., 2022, 51, 6132–6154.
6. S. Aitipamula and G. Bolla, Mol. Pharmaceutics, 2024, 21, 2153–2167.
7. H. Mei, N. Wang, D. Wu, Q. Rong, X. Bai, X. Huang, L. Zhou, T. Wang and H. Hao, Pharm. Res., 2024, 41, 142–153.
8. M. Singh, H. Barua, V. G. S. Jyothi, M. R. Dhondale, A. G. Nambiar, A. K. Agrawal, P. Kumar, N. R. Shastri and D. Kumar, Pharmaceutics, 2023, 15, 1269.
9. M. C. Etter and G. M. Frankenbach, Chem. Mater., 1989, 1, 10–12.
10. G. R. Desiraju, Acc. Chem. Res., 2002, 35, 565–573.
11. Z.-Q. Wang, Y. Liu, C. Wang, S. Jiang, Y. Feng, T. Huang and W. Huang, Atmos. Environ., 2021, 264, 118671.
12. K. Trzeciak, M. K. Dudek and M. J. Potrzebowski, Chem. – Eur. J., 2024, 30, e202402683.
13. J. Yang, S. Li, H. Zhao, B. Song, G. Zhang, J. Zhang, Y. Zhu and J. Han, J. Phys. Chem. A, 2014, 118, 10927–10933.
14. R.-A. Mitran, S. Ioniţă, D. Lincu, E. M. Soare, I. Atkinson, A. Rusu, J. Pandele-Cuşu, C. Iordache, I. Pongratz, M. M. Pop and V. Fruth, Materials, 2024, 17, 3096.
15. K. L. Nguyen, T. Friščić, G. M. Day, L. F. Gladden and W. Jones, Nat. Mater., 2007, 6, 206–209.
16. R. Yang, Y. Li, Q. Li, L. Chen, Y. Zhou, M. Zhang, T. Li and B. Qin, J. Mol. Struct., 2024, 1309, 138409.
17. T. Friščić, I. Halasz, P. J. Beldon, A. M. Belenguer, F. Adams, S. A. J. Kimber, V. Honkimäki and R. E. Dinnebier, Nat. Chem., 2013, 5, 66–73.
18. K. Užarević, I. Halasz and T. Friščić, J. Phys. Chem. Lett., 2015, 6, 3817–3826.
19. P. A. Julien and T. Friščić, Acc. Chem. Res., 2022, 55, 1295–1305.
20. H. Yada, M. Nagai and K. Tanaka, Chem. Phys. Lett., 2008, 464, 166–170.
21. B. Qin, J. Qiu, R. Yang, Y. Gan, H. Zhong, Q. Liao, S. Lu and Y. Li, IEEE Trans. Terahertz Sci. Technol., 2024, 14, 476–483.
22. Y. Bo, J. Fang, Z. Zhang, J. Xue, J. Liu, Z. Hong and Y. Du, Pharmaceutics, 2021, 13, 1273.
23. P. Wang, J. Zhao, Y. Zhang, Z. Zhu, L. Liu, H. Zhao, X. Yang, X. Yang, X. Sun and M. x. He, Int. J. Pharm., 2022, 620, 121759.
24. P. Wang, Y. Li, W. Han, Y. Yan, C. Zhang, Q. Qu, X. Zhang, L. Liu, X. Sun, X. Yang and M. He, Talanta, 2024, 278, 126489.
25. Y. Jing, J. Zhang, M. Wan, J. Xue, J. Liu, J. Qin, Z. Hong and Y. Du, IEEE Trans. Terahertz Sci. Technol., 2024, 14, 152–161.
26. E. P. Parrott, J. A. Zeitler, T. Friščić, M. Pepper, W. Jones, G. M. Day and L. F. Gladden, Cryst. Growth Des., 2009, 9, 1452–1460.
27. S. Huang, H. Deng, X. Wei and J. Zhang, Front. Bioeng. Biotechnol., 2023, 11, 1249767.
28. J. Lee, C. K. Leung, M. Ma, J. Ward-Berry, S. Santitewagun and J. A. Zeitler, J. Infrared, Millimeter, Terahertz Waves, 2023, 44, 795–813.
29. L. Macchietti, L. Casali, F. Emmerling, D. Braga and F. Grepioni, RSC Mechanochem., 2024, 1, 106–115.
30. R. W. Heckel, Trans. Metall. Soc. AIME, 1961, 221, 671–675.
31. K. L. Nguyen, T. Friščić, G. M. Day, L. F. Gladden and W. Jones, Nat. Mater., 2007, 6, 206–209.
32. A. A. L. Michalchuk and F. Emmerling, Angew. Chem., Int. Ed., 2022, 61, e202117270.
33. K. Shirzad and C. Viney, J. R. Soc. Interface, 2023, 20, 20230242.
34. M. Ferguson and T. Friščić, Phys. Chem. Chem. Phys., 2024, 26, 9940–9947.
35. T. Friščić, I. Halasz, P. J. Beldon, A. M. Belenguer, F. Adams, S. A. J. Kimber, V. Honkimäki and R. E. Dinnebier, Nat. Chem., 2013, 5, 66–73.
36. L. M. Vas, E. Slezák, K. Molnár and F. Ronkay, Thermochim. Acta, 2024, 741, 179836.
37. H. Bekaert, R. Dufour, I. Nopens, C. Vervaet, I. Nopens and T. De Beer, Int. J. Pharm.: X, 2021, 3, 100110.
38. S. A. Ajibade, L. Catalano, J. Kölbel, D. M. Mittleman and M. T. Ruggiero, J. Phys. Chem. Lett., 2024, 15, 5549–5555.
39. Y. Huang, J. Dai, S. Chang, M. Chen, X. Chen, M. Tan and F. Luo, Front. Chem., 2023, 11, 1215103.
40. C. Yuan, Y. Ji, Y. Bai, L. Li, X. Zhang, J. Gong, D. Tang, Y. Sun, C. Su and K. Zhang, et al. , CrystEngComm, 2024, 26, 3971–3983.
41. Y. Xiao, C. Wu, X. Hu, K. Chen, L. Qi, P. Cui, L. Zhou and Q. Yin, Cryst. Growth Des., 2023, 23(6), 4680–4700.

---