Patterned crystal growth and heat wave generation in hydrogels

The crystallization of metastable liquid phase change materials releases stored energy as latent heat upon nucleation and may therefore provide a triggerable means of activating downstream processes that respond to changes in temperature. In this work, we describe a strategy for controlling the fast, exothermic crystallization of sodium acetate from a metastable aqueous solution into trihydrate crystals within a polyacrylamide hydrogel whose polymerization state has been patterned using photomasks. A comprehensive experimental study of crystal shapes, crystal growth front velocities and evolving thermal profiles showed that rapid growth of long needle-like crystals through unpolymerized solutions produced peak temperatures of up to 45˚C, while slower-crystallizing polymerized solutions produced polycrystalline composites and peaked at 30˚C due to lower rates of heat release relative to dissipation in these regions. This temperature difference in the propagating heat waves, which we describe using a proposed analytical model, enables the use of this strategy to selectively activate thermoresponsive processes in predefined areas.

. Fitting parameters for additive adsorption model. 8 Table S4. Fitting parameters for pore penetration model.  Figure S6. Spatial temperature profiles of samples suspended in air. 12 Figure S7. Formulation of heat transfer problem.
13 Table S5. Fitting parameter values from global fits.
14 Table S6: Calculation of reference values from Table S5 using

S1. Supplementary video captions
Video S1. Crystallization of a metastable solution polymerized through a photomask. Sample was exposed to UV light in all areas outside the letter H, which was covered by a mask. Precursor solution composition: 2.8 M acrylamide, 26 mM N,N'-methylenebisacrylamide, 7.0 M sodium acetate, 2 mM a-ketoglutaric acid. This video also demonstrates the procedure for crystallizing samples suspended in air; a paper grid supports the sample. Ruler shown for scale; numbered increments are centimeters. Video is shown in real time.
Video S2. Video microscopy (10x) through crossed linear polarizers showing progression of crystal growth along interface between masked and unmasked regions. Video shown at 10x magnification and 0.2x speed with automatic exposure. First frame shows faint phase boundary between masked (left) and unmasked (right) regions before growth front reaches field of view. Crystals grow rapidly through masked region with large, needle-like domains, then slowly outward through unmasked region where additives have been photopolymerized, forming small crystallites. Higher magnification view of unmasked region is shown in Figure S1; electron micrographs of similar samples are shown in Figure S2. After some time has passed, a secondary cascade of small crystallites engulfs the large needles in the masked region; this corresponds to the growth of bright white areas visible in Figure 4A and

S2. Precursor solution compositions by weight and molarity
When all components are fully dissolved, each solution has a final volume of 10 mL.    Table S2. Approximate diameters of pores or crystalline domains in polymer solutions with a range of acrylamide concentrations. Extracted from SEM images such as those in Figure S2, plotted in Figure Table S2 fit with Equations S6 (red dotted line) and S8 (blue dashed line). Figure S4. Fits of various models described in this section to the data shown in Figure 2A.

A. Additive adsorption mechanisms
One possible mechanism by which additives might inhibit crystal growth more strongly when polymerized than as monomers in solution involves binding to the surface of the growing crystal, blocking its progression. We have fit our sodium acetate trihydrate crystal growth front velocity data in the presence of unpolymerized acrylamide to three leading models of crystal growth kinetics in the presence of adsorbing additives ( Figure S4). These models differ slightly in their assumptions and functional forms as described below; an excellent overview of these models is provided in a monograph by Sangwal. 1 Each model expresses the average linear step growth velocity in the presence of additives (v) normalized by the velocity in the absence of impurities (v0) in terms of the additive concentration (cadd), the affinity of the additive for the surface (K), and an effectiveness factor (a). The normalized average step velocity v/v0 is assumed to equal the normalized average face growth velocity via the Burton-Cabrera-Frank screw dislocation mechanism. 2 First, the model developed by Cabrera and Vermilyea 3 (Equation S1) assumes that additive particles adsorb to the surface terrace and are immobile compared to the velocity of the step ledges. These particles prevent the progression of step growth when they are spaced close enough together to prevent a two-dimensional crystal nucleus of critical size from squeezing between them. A geometric mean of the step velocity is used. The effectiveness factor = , where r2D* is the critical radius of a two-dimensional nucleus, l is the average distance between adsorption sites on the surface, a is the diameter of a crystallizing molecule, gl is the linear edge free energy (related to the surface free energy g by gl = ga 2 ), kB is the Boltzmann constant, T is the temperature, and s is the relative supersaturation of the crystallizing species. In this model, 0 < a < ¥; additives with a > 0.5 can arrest crystal growth entirely at sufficient concentration. Using a Langmuir adsorption isotherm, the normalized velocity is expressed as Equation S1 .
Next, a similar model developed by Kubota and Mullin 4 instead uses an arithmetic mean of the step velocity; further analysis by Sangwal asserts that this model also applies in the case of mobile additives, which would adsorb preferentially to kinks in step ledges and block their further progression. 1 The effectiveness factor takes the same form as the Cabrera-Vermilyea model: = where x0 is the average distance between adsorption sites in the ledge (i.e. kinks).
Again, 0 < a < ¥; however, in this model, fully arresting growth requires that a > 1.
The earlier model developed by Bliznakov 5,6 and refined by Chernov 7 is also expressed using Equation S2, except that 0 < a < 1, implying that crystal growth impeded by this mechanism can never be fully arrested. This formulation eschews a consideration of critical nuclei; instead, = where vi is the step growth rate in the presence of impurities at full coverage.
Each of these models fits well to our data from uncured solutions of acrylamide ( Figure S4, Table  S3). The models fit less well to the polymerized data sets; however, the free radical polymerization process used in this work produces a statistical distribution of molecular weights rather than monodisperse polymers. It may therefore be inappropriate to use the initial monomer concentration as cadd in these models. However, polymerization can be reasonably expected to increase the binding affinity of an additive through cooperative effects, 8,9 which is reflected in the fit parameters in Table S3.

B. Pore penetration mechanism
Another plausible mechanism by which polymeric additives may slow crystal growth requires considering the thermodynamics of crystallization in pores. When crystals are sufficiently small and have a high enough interfacial curvature (k), the interfacial energy between the crystal and liquid phases, g, contributes non-negligibly to the pressure inside the crystal pc according to pc = pL + gk, where pL is the pressure in the liquid. The normal melting point of a crystal, Tm, applies to infinitely large crystals with no curvature in contact with a liquid phase at its vapor pressure, pe. The actual melting point * / of a crystal with finite curvature in a liquid whose pressure may deviate from pe is expressed by Equation S3; this curvature-related difference is known as the Gibbs-Thomson effect and is reviewed along with other thermodynamic considerations of crystallization in porous media by Scherer 10,11 and more recently by Meldrum and O'Shaughnessy. 12 In this formula, vL and vc are the partial molar volumes of the liquid and crystalline phases and DSfus,v is the entropy of fusion of the crystal per unit volume.
To penetrate a pore of diameter x, a crystal must adopt a curvature = This reduction in the melting point reduces the degree of supercooling (∆ = * / − ) at room temperature, which is the driving force for crystallization in supercooled melts. This driving force is generally related to the crystal growth rate by a function of the form v = k(DT) g , where k is a growth rate constant and g is the "order" of the crystallization process (distinct from "order" in chemical kinetics; here g has no obvious fundamental meaning) where 1 ≤ g ≤ 2.5. 13 For constant T, constructing a growth rate function in this way results in a function of the form: The reduction of the driving force in porous media owing to the Gibbs-Thomson effect therefore lowers the crystal growth front velocity as the pore size is reduced.
Even in the absence of cross-linker, sufficiently concentrated "semi-dilute" polymer solutions have a characteristic mesh size. ("Semi-dilute," in the parlance of polymer physics, indicates that the polymer fraction exceeds a threshold required for coils to overlap, but is much less than 1). The mesh size x of semi-dilute linear polymer coils in good solvent generally scales with the polymer fraction F of the polymer according to ≅ Φ Fitting Equation S7 to the rate data where polymers were present while setting Tm -T to 38 K and allowing k, A, and g to vary yields the red curve shown in Figure S4 and the parameters shown in Table S4. The resulting curve is very similar to the results of fitting the Kubota-Mullin model to the polymerized data. However, the trend in crystal and pore sizes observed via scanning electron microscopy ( Figure  S2,3, Table S2) as a function of acrylamide concentration is poorly fit by Equation S6 (Figure S3, Tm -T to 38 K, b to 0.845 (derived from the empirical fit above) and allowing k, A, and g to vary yields the blue curve shown in Figure S4 and the parameters shown in Table S4. The quality of this fit is noticeably worse than that of Equation S7; other functions that empirically fit the size data as a function of acrylamide concentration in Figure S3 and Table S2 well perform similarly poorly when substituted into Equation S5 and fit to the rate data.
Note that in this section we have considered salt hydrates as supercooled melts. Similar formulations exist for crystal growth from supersaturated solutions (e.g. v = k(Dc) g ), where the driving gradient is the degree of supersaturation. Salt hydrates such as sodium acetate trihydrate, in which solvating water co-crystallizes with the solute to form solid crystals that melt at 58˚C at standard pressure, may be fairly considered in either way; in this case, it is simpler to consider them as melts.

C. Growth arrest and stochastic nucleation mechanism
The descriptions above present means by which the growth rates of single crystals may be suppressed in the presence of polymers. However, Figures 1C, S2 and S3 make clear that polymerized samples form highly polycrystalline composites with a crystallite size range that decreases with increasing polymer concentration. A mechanism that takes this phenomenon into account may explain the observed growth rate trends.
Invoking an argument first presented by Chernov, 7 Asenath-Smith et al. 15 have proposed that crystals growing into polymer fibers may experience an abrupt increase in resistance to mass transport to the interface. As sodium acetate crystals grow in a needle shape, the active growth face has a small area; small crystallites in particular may be fully occluded upon contact with the polymer mesh and halt their lengthwise growth completely. In such a case, the growth front would proceed only due to the continuous stochastic secondary nucleation of new crystallites.
We have simulated this scenario using MATLAB. To initialize the simulation, a predetermined polymer volume fraction is distributed across random voxels throughout a sample volume and initial crystal "seeds" are planted at one end. A loop of timesteps follows; during each timestep, a) each active crystal grows by a length unit, b) each crystal that has collided with a polymer, another crystal, or the sample volume boundary is rendered inactive, and c) each voxel occupied by a crystal has a set probability of nucleating a new crystal with a random orientation. A snapshot of a simulated volume containing polymers (red) and crystals (blue) is shown in Figure S5A. At each timestep, the position is recorded of the farthest plane from the starting point in which a threshold percentage of the voxels have been converted to crystal; thresholds of 10%, 1%, 0.1%, 0.01% are shown in black. By inspection, after a sufficient induction time has passed, the average rate of progression of these threshold planes are constant and equal to one another regardless of the threshold percentage ( Figure S5B). Induction times increase with increased threshold percentage; the velocities of higher threshold planes over time fluctuate less than in lower thresholds. In order to analyze front velocity trends, we tracked the velocity of the 1% threshold plane during the final half of growth simulations which ran for more than twice this induction time. We used a 100 x 100 x (variable length) sample volume; 25 initial seeds were evenly placed along the initial plane during initialization and randomly oriented at angles of more than 45˚ from the direction of front progression. We ran 3 simulations for each set of parameters. The two variables that affect the growth rate in simulated systems containing polymers relative to simulated systems without polymers are the polymer volume fraction and the probability of nucleation per timestep. The effects of these parameters on front velocity are shown in Figure  S5C. Growth front rates are depressed at increased polymer volume fraction; this dependence is reduced as nucleation rates are raised. Figure S5C is overlaid with experimental data from Figure  2A (assuming 100% polymerization and that the density of polyacrylamide in solution is the same as that of solid acrylamide monomers). While no one parameter set fits the experimental data perfectly, the overall trend of a monotonically decreasing front velocity with increasing polymer concentration is captured.

D. Viscous mechanism
Finally, polymerizing the acrylamide in a solution drastically increases its viscosity (captured in the denominator of the parameter k in Equation S5). The various mass transfer processes involved in crystal growth (e.g. diffusion and orientation of solute ions) may be slower in a more viscous solution, slowing the velocity of the growth front. However, polymer solutions are complex fluids in which viscous interactions are dependent on both length and time scales; 16 untangling this contribution to crystal growth inhibition is non-straightforward. It may be that the processes involved in crystal growth take place on small enough characteristic length scales to be unaffected by the polymers' contribution to the viscous environment; a past analysis of ice growth in hydrocolloid suspensions discarded viscosity as the main driver of inhibition. 17

S5. Derivation and fitting of temperature profiles
The following derivation takes as its starting point the formulation that Hopper and Uhlmann (1973) solved for the one-dimensional temperature distribution around the moving crystal-solution interface during crystal growth at a constant velocity. 18 These authors found a steady state solution for the temperature profile within the liquid phase, but none exists for the temperature profile in the solid phase under the conditions they describe, which includes no heat transfer outside the system. This result is consistent with our observations of crystallizing systems suspended in air (Figure S6), in which temperatures at and behind the moving front rise over time. However, we observed that crystallizing systems on a cold plate held to a constant temperature quickly reach a temperature profile around the solid-liquid interface that remains stable over time (Figure 3A,B). When formulating the heat transfer problem under these conditions, we include a dissipation term representing conductive heat transfer to the cold plate (normal to the direction of crystal growth) that was absent in Hopper and Uhlmann's analysis; this addition gives rise to a steady state solution (Equations S17 and S18). In this analysis, we derive the one-dimensional temperature distributions in the solid and liquid phase (Ts(x), Tl(x)) in a crystallizing supercooled system on a cold plate held at temperature T¥ with a frame of reference centered around a solid-liquid boundary (x = 0) moving with constant velocity v ( Figure S7). The system is solid in the negative domain (-¥ < x < 0) and liquid in the positive domain (0 < x < ¥). Equations S10 and S11 quantify heat conduction, equations S12-S14 are boundary conditions, and equation S15 is a heat balance at the interface. In the equations below, a (m 2 s -1 ) is the thermal diffusivity, b (K) is the latent heat scaled by specific heat capacity, and g (s -1 ) is the thermal contact conductance of the interface between the system and the cold plate scaled by heat capacity; for simplicity, these parameters are assumed to be identical between the solid and liquid states.
Equations S10 and S11 each have solutions of the form shown in S16 where C1 and C2 are constants: After applying boundary conditions S12-S15, equations S17 and S18 emerge as solutions: Globally fitting these equations to the profiles at the final timepoints shown in Figure 3A,B with v fixed to the recorded front velocities (In Figure 3A, v = vunpoly = 0.896 mm s -1 ; in Figure 3B, v = vpoly = 0.131 mm s -1 ) and a, b, and g allowed to vary (but shared between data sets) yielded the dashed curves in the figures and fitting parameter values listed in Table S5, which displayed reasonable agreement with literature values (which were calculated based on parameters collected for sodium acetate trihydrate in somewhat different conditions, making some discrepancy unavoidable; see Table S6). Profiles for a variety of velocities with a, b, and g fixed to the fit values are shown in Figure 3C. The sharp peak present in the model fit is blunted in the experimental data for the unpolymerized sample; this may be explained by a rough (non-discrete in x) advancing front or the diffusion of this steep gradient within the thickness of the cover glass between the crystallizing solution and the infrared camera.  It is clear that varying the crystal growth front velocity alone -which can be achieved by varying the composition and polymerization state of the additives to the solution -can be responsible for large changes in the temperature profile of the system. The maximum temperature can be expressed as a function of crystal growth front velocity using Equation S19. This is plotted in Figure S8; it can be observed that lim (→?
/#-= ? + = 55˚C. Based on prior reports of crystallization in metastable sodium acetate solutions diluted beyond a 1:3 salt:water stoichiometric mixture, 21,22 we expect this temperature to approach but never reach sodium acetate trihydrate's melting point of 58 ˚C, which agrees with this result.