Swelling thermodynamics and phase transitions of polymer gels

TL;DR

This paper reviews the thermodynamics and phase transitions of polymer gels, highlighting how shear rigidity influences phase behavior and introduces the concept of extreme thermodynamics for large shape changes.

Contribution

It provides a comprehensive analysis of gel phase transitions, emphasizing the role of shear rigidity and introducing the novel concept of extreme thermodynamics for shape change.

Findings

Volume phase transition resembles mean-field binary mixture models

Shear rigidity causes failure of analogy at critical point

Phase separation can induce macroscopic shape change

Abstract

We present a pedagogical review of the swelling thermodynamics and phase transitions of polymer gels. In particular, we discuss how features of the volume phase transition of the gel's osmotic equilibrium is analogous to other transitions described by mean-field models of binary mixtures, and the failure of this analogy at the critical point due to shear rigidity. We then consider the phase transition at fixed volume, a relatively unexplored paradigm for polymer gels that results in a phase-separated equilibrium consisting of coexisting solvent-rich and solvent-poor regions of gel. Again, the gel's shear rigidity is found to have a profound effect on the phase transition, here resulting in macroscopic shape change at constant volume of the sample, exemplified by the tunable buckling of toroidal samples of polymer gel. By drawing analogies with extreme mechanics, where large shape changes are achieved via mechanical instabilities, we formulate the notion of extreme thermodynamics, where large shape changes are achieved via thermodynamic instabilities, i.e. phase transitions.

Figures

Equations

$$\mathrm{d}G$$

$$\mathrm{d}F$$

$$V\approx n_s{v}_s+n_m{v}_m,$$

$$V=V_0(\det \Lambda ).$$

$$F(T,\Lambda ,n_s;\lambda )=F(T,\Lambda ,n_s)+\lambda \left[V_0\det \Lambda -V(n_s)\right],$$

$$\varphi \equiv \frac{n_m{v}_m}{n_s{v}_s+n_m{v}_m}\approx 1-\frac{v_s{n}_s}{V},$$

$$F=\frac{v{n}_m}{\varphi }\mathcal{F}(T,\Lambda ,\varphi ),$$

$$F=F_{\mathrm{sol}}+\Delta F,$$

$$\Pi (T,\varphi ,\Lambda )=-{\left(\frac{\partial ({\varphi }^{-1}\Delta \mathcal{F})}{\partial ({\varphi }^{-1})}\right)}_{T,\Lambda },$$

$$Z=\int [{\mathrm{d}}^3{\mathbf{\sigma }}_i{]}_{i=1}^{n_s}\int [\mathcal{D}{\mathbf{R}}_j(s){]}_{j=1}^{n_{ch}}{e}^{-\beta E}\prod _{k=1}^{3n_c}\delta ({\mathbf{f}}_k^{\mathrm{network}}),$$

$${\mathbf{f}}_k^{\mathrm{network}}\equiv \sum _{i\ne j}^{n_{ch}}{a}_k^{ij}[{\mathbf{R}}_i(L_{ch})-{\mathbf{R}}_j(0)]=0,$$

$$\overline{\mathcal{E}}\approx \frac{v{k}_{B}T}{2}\left[{\chi }_{m-m}{\rho }_m^2+2{\chi }_{m-s}{\rho }_m{\rho }_s+{\chi }_{s-s}{\rho }_s^2\right],$$

$$\eqalign{\overline{\mathcal{E}}\approx\frac{k_{B}T}{2v}\big{[}&\chi_{m-m}\,\phi^{2}\\ &+2\chi_{m-s}\,\phi(1-\phi)+\chi_{s-s}\,(1-\phi)^{2}\big{]}\,.}$$

$$\Delta E=N{k}_{B}T\chi \varphi (1-\varphi ),$$

$${\Omega }_1(\mathbf{r})=\frac{V}{v}\prod _{i=1}^3\frac{N_i!}{\left(\frac{N_{i}a+r_i}{2a}\right)!\left(\frac{N_{i}a-r_i}{2a}\right)!},$$

$$\eqalign{{\rm ln}\,\Omega_{1}&\approx{\rm ln}\,\frac{V}{v}+\sum_{i}\bigg{[}N_{i}{\rm ln}\,N_{i}-\frac{N_{i}a+r_{i}}{2a}{\rm ln}\,\frac{N_{i}a+r_{i}}{2a}\\ &\mkern 124.0mu-\frac{N_{i}a-r_{i}}{2a}{\rm ln}\,\frac{N_{i}a-r_{i}}{2a}\bigg{]}\,.}$$

$$\ln {\Omega }_1\approx \ln \frac{V}{v}+N\ln 2-\frac{3|\mathbf{r}{|}^2}{2R_0^2},$$

$$S_1(\mathbf{r})=S_0-\frac{3k_B|\mathbf{r}{|}^2}{2R_0^2}+k_B\ln \frac{V}{v}$$

$$\Delta S_1(\mathbf{r})=-\frac{3k_B}{2R_0^2}\left({\Lambda }_{ki}{\Lambda }_{kj}-{\delta }_{ij}\right)r^i{r}^j+k_B\ln (V^{\prime }/V).$$

$$\eqalign{\Delta S_{n_{ch}}&=-\frac{3k_{B}}{2R_{0}^{2}}\left(\Lambda_{ki}\Lambda_{kj}-\delta_{ij}\right)\big{(}r_{1}^{i}r_{1}^{j}+\cdots\\ &\mkern 206.0mu\cdots+r_{n_{ch}}^{i}r_{n_{ch}}^{j}\big{)}\\ &\mkern 32.0mu+n_{ch}k_{B}{\rm ln}\,{\rm det}\,\Lambda\,,}$$

$$⟨r_i{r}_j⟩=\int {\mathrm{d}}^{3}r{P}_1(\mathbf{r})r_i{r}_j=\frac{R_0^2}{3}{\delta }_{ij},$$

$$\Delta S_{\mathrm{net}}=-\frac{1}{2}{n}_{ch}{k}_B\left[\mathrm{tr}{\Lambda }^T\Lambda -3-2\ln \det \Lambda \right].$$

$$\Delta S_{\mathrm{net}}=-\frac{1}{2}{n}_{ch}{k}_B\left[\mathrm{tr}{\Lambda }^T\Lambda -3-\ln \det \Lambda \right],$$

$$\Omega =\frac{N!}{n_A!n_B!},$$

$$S_{\mathrm{latt}}\approx k_B\left[N\ln N-n_A\ln n_A-n_B\ln n_B\right].$$

$$S_{\mathrm{latt}}\approx -N{k}_B\left[(1-\varphi )\ln (1-\varphi )+\varphi \ln \varphi \right],$$

$$S_{\mathrm{latt}}\approx -N{k}_B\left[(1-\varphi )\ln (1-\varphi )+\frac{\varphi }{\mathcal{N}}\ln \frac{\varphi }{\mathcal{N}}\right].$$

$$\eqalign{\Delta S_{{\rm mix}}&=V\mathcal{S}(\phi)-vn_{s}\mathcal{S}(\phi=0)-vn_{m}\mathcal{S}(\phi=1)\\ &\approx-Nk_{B}\left[(1-\phi)\,{\rm ln}\,(1-\phi)+\frac{\phi}{\mathcal{N}}\,{\rm ln}\,\phi\right]\;,}$$

$$\Delta S_{\mathrm{mix}}\approx -N{k}_B(1-\varphi )\ln (1-\varphi ),$$

$$\Delta F=\Delta F_{\mathrm{elastic}}+\Delta F_{\mathrm{mix}}+\lambda \left[\det \Lambda -\frac{{\varphi }_0}{\varphi }\right],$$

$$\eqalign{\Delta F_{\rm elastic}&=-T\Delta S_{\rm net}\\ &=\frac{1}{2}n_{ch}k_{B}T\left[{\rm tr}\,\Lambda^{T}\Lambda-3-{\rm ln}\,{\rm det}\,\Lambda\right]\,,}$$

Full text

††: Nanotechnology

\ioptwocol

1 Introduction

Within the realm of amorphous, rigid materials without crystalline symmetries, polymer gels possess an interesting duality, having a rubber-like elasticity whilst being able to undergo large volume changes due to mixing with a solvent. These materials are soft, being composed of large macromolecules whose interactions are often governed by thermal fluctuations. There are three essential ingredients: polymers, solvent, and cross-links.

Polymers are composed of many short segments (i.e., monomers) that are typically chemically bonded end-to-end, as shown in figure 1(a). Whilst there are energetically favorable bond angles between successive monomers, there often are several monomer-monomer bond conformations that are at comparable energies [1]. As the number $\mathcal{N}$ of monomers that constitute a polymer is typically on the order of $10^{3}-10^{4}$, there are many different mutually accessible polymer conformations, within a small energy window, that a polymer may be found in. Thus, polymers are said to have static flexibility [2]. Furthermore, there are modest energy barriers between bond angles, enabling frequent transitions that are driven by thermal fluctuations, causing the polymer to explore many different conformations over time. As such, polymers are also said to have dynamic flexibility.

Quite generally, as a consequence of such static and dynamic flexibility, correlations between monomer-monomer bond angles decay with distance along the backbone of the polymer. Beyond a certain persistence length, bond angles are barely correlated. Therefore, conformations of polymers that span many persistence lengths have the form of a “random walk,” an example of which is shown in figure 1(b). The radius of gyration $R_{g}$ specifies the characteristic size of the polymer, as also illustrated in figure 1(b); for large $\mathcal{N}$, $R_{g}$ scales with the number of monomers $\mathcal{N}$ as $\mathcal{N}^{\nu}$, where $\nu=1/2$ for a random-walk polymer, or “ideal chain,” for which the excluded volume interaction between polymer segments is neglected. More realistically, an isolated polymer does not intersect with itself, resulting in statistics of a self-avoiding (as opposed to ideal) random walk, for which $\nu\approx 3/5$ in three dimensions, reflecting the swelling of the polymer sequence. This is indeed the situation when the polymer is immersed in a “good” solvent, one in which the polymer is miscible, as opposed to the case of immiscibility in a “poor” solvent, where the polymer radius is decreased and a compact structure is formed due to the high energetic penalty for the mixing of the polymer and solvent. An intermediate case is the $\vartheta$-solvent, where the radius-decrease due to a mildly poor solvent counteracts the radius-increase due to the self-repulsion of the polymer, resulting in ideal-chain scaling, for which $\nu=1/2$ (see, e.g., [1]).

Now consider a solution of many such polymers. Focusing on the case where all polymers are composed of roughly the same number $\mathcal{N}$ monomers that are chemically identical, the solution can be brought to a polymer concentration for which individual polymer coils overlap spatially in equilibrium. In this case, application of a static stress induces steady state flow, since after a short-time elastic response, where the polymer coils deform, they are able to re-arrange in space continually. Rigidity results from the introduction of cross-links between neighboring polymers, since cross-linked molecules can no longer individually undergo substantial re-arrangement relative to cross-linked partners. If there is sufficient linking of different polymers then a container-spanning, percolating, network of linked polymer coils forms; this constitutes a gel [1]. In a gel, the cross-linked clusters of polymers are thus localized in space relative to one another, unable to explore the volume of their container via Brownian motion; such ergodicity breaking of the polymers due to the formation of a percolating polymer network is a hallmark of the onset of rigidity [3], at least in spaces of dimension $d>2$.

Despite the large variety of cross-links that can be formed, they are generally classified according to two categories: physical and chemical 111Note that we do not consider here polymer rings that may be topologically linked to form rigid “Olympic” gels [1, 4] or other knotted polymer networks [5, 6].. Examples of physical cross-links include entangled polymers as well as non-covalent bonds, such as ionic (i.e., electrostatic) bonds. Physical cross-links enable an elastic response over potentially extended timescales; however, they are not truly rigid, as they allow flow at sufficiently long times due to the reversible nature of the cross-linking process. Instead, we shall focus our attention on gels formed from chemical cross-links, which are typically induced by the introduction of small–when compared to the typical polymer size–cross-linking molecules that form covalent bonds between polymers, as depicted in figure 1(c). Chemical cross-links may be regarded as permanent so that the network topology of the gel is frozen in after cross-linking, much like the cross-links in rubber, resulting in a thermodynamically rigid material [7]. Moreover, unlike amorphous elastic materials, such as glasses, these gels are in a well-defined, equilibrium solid phase; unlike conventional solids, gels lack long-range order. Thus, gels are equilibrium amorphous solids [3, 8].

However, unlike rubber, polymer gels are typically cross-linked in the presence of a solvent, which permeates the polymer network of the gel. The magnitude of the osmotic pressure $\Pi$ due to the mixing of polymer with the solvent is set by the thermal energy scale $k_{B}T$, and is thus on the same scale as the entropic elastic stresses of the polymer network. As a result, both are important in determining the macroscopic equilibrium state of the gel. This is especially true for gels having low cross-link densities, which have the signature ability to undergo large macroscopic volume changes in response to varying solvent conditions. In the presence of a good solvent, the polymer network is well mixed with solvent molecules, and the gel incorporates a large volume of solvent and is said to be swollen; in the presence of a poor solvent, the polymer network is essentially segregated from the solvent molecules and is said to be deswollen. Suitably prepared gels have a remarkably large volume response, capable of swelling to an equilibrium volume on the order of $10^{3}$ times their deswollen volume by absorbing solvent [9]. Whilst a variety of different swollen volumes can be achieved by continuous changes in solubility (as induced, e.g., via temperature), certain gels can exhibit a discontinuous change in volume. For example, poly(N-isopropylacrylamide) (pNIPAM) gels in water gradually deswell under heating until $\sim\mkern-4.0mu32^{\circ}{\rm C}$. Beyond this, due to a change in solvent nature from good to poor, they abruptly expel most of their solvent [9]. In fact, there is a first-order phase transition between well-defined swollen and deswollen phases of gel. This phase transition is confirmed by varying the osmotic pressure, leading to a phase diagram having a first-order transition region separating the two phases, terminating at a critical point [10, 11].

In addition to swelling, polymer gels can undergo shape changes. Mechanical constraints, such as attachment to a stiff substrate, can frustrate homogeneous deswelling, resulting in inhomogeneous deswelling of the gel, which can lead to the formation of surface ripples [12, 13]. Gels that are subject to inhomogeneous swelling are of particular interest, as the resulting deformations typically cannot be realized in flat space [14], resulting in a variety of buckled shapes [15, 16], some of which mirror patterns found in nature [17, 18, 19]. This has led to origami-inspired [20] and biology-inspired [21] searches for ways to program certain shapes that are actuated upon swelling.

In this Topical Review, we first discuss in Section 2 the thermodynamic description of polymer gel swelling. We give a brief outline of the statistical mechanical treatment due to Flory and Rehner [22, 23] and mention some subtleties that arise in describing the rubber-like elasticity of the gel [7]. We then present a pedagogical review of phase transitions in Section 3, to develop an intuition for the volume phase transition and the critical behavior of gels in analogy with the Van der Waals theory of the liquid-vapor phase transition. Continuing with this analogy, we consider in Section 4 how a transition to phase coexistence between swollen and deswollen phases can be achieved by arresting the deswelling transition of a swollen gel. In Section 5, we show how the equilibrium phase-coexistent gel is characterized by a large deformation of the macroscopic gel shape that is distinct from the usual volume phase transition. Drawing on analogies with the extreme mechanics of shape-changing materials through programmed mechanical instability [24, 25], we propose that phase-coexistent gels provide a route to large deformation via thermodynamic instability. To provide an illustrative example of this extreme thermodynamics, we give a detailed description of the phase-coexistent equilibrium of gel toroids and the accompanying shape-buckling transition, which has been realized in experiments. Finally, in Section 6, we summarize some of the open problems in the field and highlight some of the gaps in the understanding of polymer gels that have yet to be filled.

2 Swelling thermodynamics

The thermodynamic description of polymer gels is based on that of non-ideal fluids, in which interactions between particles give rise to equations of state, such as the Van der Waals equation, that differ from the universal, ideal gas description. To begin, we consider the state functions that are required to describe the macroscopic state of the gel. Then we turn to the microscopic description and outline the Flory-Rehner [22, 23] mean-field theory which yields approximate equations of state of the gel that are analogous to the Van der Waals equation of state. Next, just as the Van der Waals equation predicts that a fluid expands with increasing temperature at constant pressure, we show how the Flory-Rehner equation of state predicts gel deswelling with increasing temperature. Finally, we examine the breakdown of thermodynamic stability predicted by the Flory-Rehner theory, highlighting analogies and differences with phase transitions in fluids.

2.1 State functions and thermodynamic potentials

Macroscopic materials, such as polymer gels, are composed of a vast number of microscopic degrees of freedom that are in continual flux, i.e., are thermally fluctuating. In a thermodynamic description of such materials, these many fluctuating microscopic degrees of freedom are averaged over time and space to yield state functions (see, e.g., [26]). For example, in a fluid consisting of $N$ identical particles occupying a fixed volume $V$, both $V$ and $N$ are state functions. Whilst the spacings between particles are not fixed, there are, on average, $\rho\equiv N/V$ particles per unit volume. If the fluid is isolated then its total energy $E$ is fixed, as are the total number $N$ of particles in the fluid as well as its total volume $V$. These three state functions are sufficient for characterizing the macroscopic equilibrium state of the system. In order to quantify what happens when the macroscopic degrees of freedom $(E,V,N)$ are changed, one employs a thermodynamic potential. The entropy $S(E,V,N)$ is one example of a thermodynamic potential, which has the fundamental property that if the energy, volume, or number constraints are relaxed, the equilibrium state that the system eventually attains corresponds to one of maximum entropy. Note that entropy is also a state function, corresponding to the number of microstates of the fluid that give rise to a fixed macrostate $(E,V,N)$. At times, it is useful to use the entropy as a state-characterizing function, exchanging it with the total energy $E$, which can then take the role of the thermodynamic potential, corresponding to a macroscopic description $(S,V,N)$.

It is often convenient to consider interactions between the fluid and a much larger “bath,” whose state is not affected by the presence of the fluid. If we imagine that the fluid is kept in a container that allows heat to flow between the fluid and the surrounding bath then the energy of the fluid and that of the bath are allowed to change. The total entropy $S=S_{\rm fluid}+S_{\rm bath}$ is maximized when the temperature $T$ of the fluid matches that of the bath, which characterizes a state of thermal equilibrium. This container can either maintain a fixed volume $V$ of the fluid or be flexible, in which case mechanical equilibrium is reached when the pressure $P$ of particles in the fluid is balanced by a similar pressure from the bath. Similarly, the container can either be impermeable, maintaining a constant number $N$ of particles, or it can be permeable, so that chemical equilibrium is reached when the chemical potential $\mu$ of the fluid matches that of the surrounding bath. In this way, the paired state functions $(S,T)$, $(V,P)$, and $(N,\mu)$ are considered conjugate to one another. Much like the density $\rho$ of the fluid, $(T,P,\mu)$ are intensive state functions that characterize material properties of the fluid, whereas $(S,V,N)$ are extensive state functions. Whilst there is freedom in choosing the three state functions that describe the macroscopic state of the fluid,222Recall, however, that due to the Gibbs-Duhem equation, which provides a link between the intensive parameters, at least one of the state functions must be extensive. let us consider the temperature $T$ and the number of particles $N$ as specified properties, i.e., constraints imposed on the fluid. There are two representations that can be considered: the Gibbs representation $(T,P,N)$ and the Helmholtz representation $(T,V,N)$, with corresponding thermodynamic potentials $G$, the Gibbs free energy, and $F$, the Helmholtz free energy. Changes in constraints lead to changes in the thermodynamic potential, described by a Gibbs equation for each representation, namely

[TABLE]

from which we see that the two potentials are related via a Legendre transform, resulting in the relation $G=F+PV$.

Now consider a sample of gel that is allowed to exchange solvent with its surroundings but contains a constant number of monomers (i.e., polymer segments). The gel is composed of $n_{s}$ solvent molecules, $n_{m}$ monomers, and $n_{c}$ cross-linking molecules. Thus, it is natural to assume that in the Gibbs representation the state of the gel is characterized by the state functions $(T,P,n_{s},n_{m},n_{c})$. However, we will assume that each molecule and monomer occupy volumes $v_{s}$ and $v_{m}$, respectively, so that the volume $V$ of the system is approximately given by

[TABLE]

i.e., we have neglected the very small contribution due to the cross-linking molecules since the number of cross-links is typically orders of magnitude smaller than the number of solvent molecules and monomers. Therefore, changing the pressure $P$ acts to change the volume per solvent molecule ($v_{s}$) and the volume per monomer ($v_{m}$). This, however, only happens at very high pressures and is not of significance in the situations of interest here. We will instead focus on the effect that mixing these two chemical species has on the macroscopic properties of polymer gels, and treat $v_{s}$ and $v_{m}$ as constants; for simplicity, we assume that they have the same value, namely $v_{s}\approx v_{m}\equiv v$. Furthermore, the number of monomers $n_{m}$ and the number of cross-links $n_{c}$ are imposed at the formation of the gel, and are also assumed constant. Thus, we are left with the state functions $(T,n_{s})$, where the volume $V$ is determined as a function of $n_{s}$ via equation (1b); this state characterization then amounts to a Helmholtz representation of the gel.

Unlike fluids, however, polymer gels possess a nonzero rigidity with respect to elastic deformations. Therefore, in addition to occupying a volume $V$, the gel is able to maintain a deformed shape indefinitely when subjected to stress. We therefore require additional state functions to account for this fact. One such deformation is the change in the three side-lengths $\{L_{1},L_{2},L_{3}\}$ of the box-shaped sample of gel shown in figure 2 to lengths $L_{i}^{\prime}=\Lambda_{i}L_{i}$. The deformation of the gel at constant volume is therefore set by the dimensionless ratios of length $\{\Lambda_{1},\Lambda_{2},\Lambda_{3}\}$ such that $\Lambda_{1}\Lambda_{2}\Lambda_{3}=1$. Note that specifying the side-lengths of a parallelepiped region of gel is but one example deformation that can be achieved. For general gel shapes, it is more appropriate to examine how the distance between any two points $\mathbf{r}$ and $\mathbf{r}+{\rm d}\mathbf{r}$ is altered upon deformation of the gel, which takes ${\rm d}\mathbf{r}$ to ${\rm d}\mathbf{R}$. For example, it is useful to imagine $\mathbf{r}$ and $\mathbf{r}+{\rm d}\mathbf{r}$ as two neighboring cross-links. Assuming affine deformations, for which the changes in lengths between representative points in the gel are independent of position, all lengths are transformed by a deformation matrix $\Lambda$, such that ${\rm d}R_{i}=\Lambda_{ij}{\rm d}r_{j}$, where we use and adopt Einstein’s summation convention over repeated indices. Deformed volume elements ${\rm d}^{3}R$ are related to the undeformed ones via ${\rm d}^{3}R=({\rm det}\,\Lambda){\rm d}^{3}r$, so ${\rm det}\,\Lambda$ is the ratio of the deformed volume to the undeformed volume. Thus, deformations that maintain the gel volume are characterized by ${\rm det}\,\Lambda=1$. Alternatively, we are free to choose a reference state, which we shall refer to as a reference configuration $\mathcal{R}$, where the volume of the gel is given by $V_{0}$, such that after a deformation of the gel, the volume of the deformed state, which we shall refer to as a target configuration $\mathcal{T}$, is given by

[TABLE]

Therefore, the determinant $({\rm det}\,\Lambda)$ may be expressed in terms of the amount of solvent $n_{s}$ in configuration $\mathcal{T}$ via equations (1b) and (1c). As cross-links undergo Brownian motion, some care has to be taken in relating macroscopic affine deformation to a corresponding microscopic deformation. However, it has been found [8] that average cross-link positions indeed undergo affine deformation.

In order to account for the effect of deformation on the equilibrium thermodynamics of the gel, it is necessary to introduce the deformation matrix $\Lambda$ as a state function. However, by equation 1c, the deformation matrix determines the volume of the gel. To account for this redundancy in state functions, we may express the Helmholtz free energy as

[TABLE]

where $\lambda$ a Lagrange multiplier accounting for the constraint associated to equation (1c). In this form, the Lagrange multiplier $\lambda$ is an additional state function and the constraint is an equation of state.

It is useful to define the polymer volume fraction $\phi$ via

[TABLE]

i.e., the fraction of the gel volume that is occupied by polymer; $\phi=1$ corresponds to a gel that is completely devoid of solvent, whereas $\phi=0$ is the limit of an infinitely dilute gel. Noting that since $F$ is a homogeneous first-order function in its extensive parameters [26], we can define $F(T,\Lambda,n_{s})=V\mathcal{F}(T,\Lambda,n_{s}/V)$, where $\mathcal{F}$ is a free energy density. Therefore, in terms of the polymer volume fraction $\phi$, we have

[TABLE]

where we have used the assumption $v_{s}\approx v_{m}\equiv v$.

Inasmuch as $P$ and $V$ are conjugate to each other for a fluid, for a gel there is a state function that is paired with the polymer volume fraction $\phi$. This is the osmotic pressure $\Pi$. If the gel is in equilibrium with a solvent bath, the chemical potential of the solvent in the gel, $\mu(T,P)$, equals the chemical potential of the solvent in the bath, $\mu_{0}(T,P)$, plus a contribution $\Delta\mu$ accounting for the presence of the polymer network. Then $\mu(T,P)=\mu_{0}(T,P)+\Delta\mu$. In addition, we may regard the boundary of the polymer network as a semipermeable membrane. Equilibrium then requires an additional pressure in order to maintain the imbalance in solvent concentration in and out of the gel, $\mu(T,P+\Pi)=\mu_{0}(T,P)$; this additional pressure $\Pi$ is, by definition, the osmotic pressure. It can be shown (see, e.g., [27]) that the osmotic pressure $\Pi$ is related to $\Delta\mu$ via $\Pi=-\Delta\mu/v$, where $v$ is the solvent particle volume.

The thermodynamics of the polymer gel is determined by how polymer mixes with solvent. Just as $\Delta\mu$ is the change in the chemical potential of the solvent due to the presence of the polymer network, we can decompose the total free energy $F$ as

[TABLE]

where $F_{\rm sol}$ is the part of the free energy due to solvent, without the effect of the polymer network. Therefore, $\Delta\mu=(\partial\Delta F/\partial n_{s})_{T,\Lambda}$. Then, using the relation between $n_{s}$ and the volume fraction $\phi$ as well as equation (1f), the osmotic pressure $\Pi$ is given by

[TABLE]

where it is evident that if $\Pi$ is analogous to a pressure then $1/\phi$ plays the role of volume. Note that, since the determinant of the deformation matrix $\Lambda$ depends on the polymer volume fraction via the volume-deformation relation (1c), it is important to enforce the Lagrange multiplier constraint in equation (1d) when computing the osmotic pressure. As $\Delta F$ is the part of the free energy that describes gel deformation, such as swelling, we shall refer to it as the deformation free energy.

2.2 Flory-Rehner equation of state

Just as we enumerated a set of macroscopic descriptors of the gel, let us consider some microscopic ones. The gel is a mixture of $n_{s}$ solvent molecules, $n_{m}$ monomers, and $n_{c}$ cross-links. Whilst the solvent molecules may have multiple internal degrees of freedom, e.g., rotational and vibrational, let us focus only on the center-of-mass degrees of freedom and treat them as point particles at positions $\{\bm{\sigma}_{i}\}$, each occupying a volume $v$, where $i$ runs from 1 to $n_{s}$. Rather than treating each of the $n_{m}$ monomers as individual particles, we will group them into polymers. For simplicity, assume that (i) we can ignore any “free-ends” or “loops” of polymers in the polymer network and consider only segments whose endpoints are cross-linked to other segments, and (ii) each of these segments, which we shall refer to as “chains,” are composed of $\mathcal{N}$ monomers. Much like the simple representation of solvent molecules, we opt for a simple representation of chains as one dimensional curves $\{\mathbf{R}_{j}(s)\}$, where $s$ is the arclength parameter, running from 0 to chain length $L_{ch}\approx v^{1/3}\mathcal{N}$, and $j$ runs from 1 to the number of chains $n_{ch}\equiv 2n_{c}\approx n_{m}/\mathcal{N}$.

In order to link these microscopic degrees of freedom to the macroscopic properties of the system, one approach is to fix temperature $T[\equiv 1/(k_{B}\beta)]$ and the number of particles of each species in the system, and determine the canonical partition function $Z$, given by

[TABLE]

where $E$ is the total potential energy of the system. The network topology is set by a collection of constraints on the $2n_{ch}=4n_{c}$ ends of the chains $\{\mathbf{R}_{i}(L_{ch}),\mathbf{R}_{i}(0)\}$. At each cross-link there are 4 ends that coincide; we may choose these ends such that 2 are at $s=0$ and 2 are at $s=L_{ch}$. However for each cross-link there are only three independent constraints; the fourth is automatically satisfied. For example, if a cross-link consists of the chains ends $\{\mathbf{R}_{1}(L_{ch}),\mathbf{R}_{2}(L_{ch}),\mathbf{R}_{3}(0),\mathbf{R}_{4}(0)\}$ then enforcing the constraints $\mathbf{R}_{1}(L_{ch})=\mathbf{R}_{3}(0)$, $\mathbf{R}_{1}(L_{ch})=\mathbf{R}_{4}(0)$, and $\mathbf{R}_{2}(L_{ch})=\mathbf{R}_{3}(0)$ automatically implies that the fourth constraint $\mathbf{R}_{2}(L_{ch})=\mathbf{R}_{4}(0)$ is satisfied. Therefore, there are $3n_{c}$ vectors that are constrained, yielding exactly $n_{c}$ independent vectors, describing the positions of cross-links in space. Thus, for each of the $n_{c}$ cross-links, there are 3 constraint equations that can be written as

[TABLE]

for $k=1\dots 3n_{c}$, where $a_{k}^{ij}$ is an adjacency matrix that is $1$ when the two polymer ends are joined by a cross-link and is [math] otherwise. These constraints are enforced by including a product of Dirac delta functions $\Pi_{k}\delta(\mathbf{f}^{\rm network}_{k})$ in the integrand of the partition function, ensuring that the only contributions to the sum over states are those where $\mathbf{f}^{\rm network}_{k}=0$ for all $k$. Note that these topological constraints pose a considerable technical difficulty in the evaluation of the partition function $Z$ and the free energy $F=-k_{B}T\,{\rm ln}\,Z$, due to the lack of a periodic structure. A mesoscopic representation of such a network with “quenched disorder” is shown in figure 1(c). However, for sufficiently large gels, there are many different mesoscopic network structures. Thus, instead of summing over polymer configurations with a certain fixed network topology, one can instead sample from a distribution of mesoscopic network structures, with the idea that they all appear somewhere in the gel; this is known as self-averaging. The partition function is then evaluated via the replica trick, where many copies or “replicas” of the gel are treated as new interacting degrees of freedom to be integrated over [7, 28, 3, 8].

However, instead of seeking a direct evaluation of the partition function $Z$ in equation (1i), we consider the classical construction of Flory and Rehner [22, 23], which amounts to a mean-field approximation. In particular, we seek a description of the deformation free energy $\Delta F=\Delta E-T\Delta S$, where $\Delta E$ is the change in the energy and $\Delta S$ is the change in the entropy due to the mixing of the solvent and polymer. The total potential energy $E$ describes the microscopic interaction energy and is approximated by the sum of three contributions: $U_{m-m}$ is the energy of monomer-monomer interactions, $U_{s-s}$ is the energy of solvent-solvent interactions, and $U_{m-s}$ is the energy of monomer-solvent interactions (see [1]), with each of these terms depending on particle positions. For example, $U_{m-m}$ contains an excluded-volume interaction $(v_{m}/2)\int{\rm d}s_{i}\int{\rm d}s_{j}\delta(\mathbf{R}_{i}(s_{i})-\mathbf{R}_{j}(s_{j}))$ between chains $i$ and $j$ (including self-interactions, corresponding to the case $i=j$). In the mean-field approximation, the interaction energy depends only on local densities of solvent and monomer, resulting in a simple form for the mean energy density $\overline{\mathcal{E}}$, namely

[TABLE]

where $\rho_{m}[\equiv n_{m}/V]$ and $\rho_{s}[\equiv n_{s}/V]$ are number-densities of monomer and solvent and $\{\chi_{m-m}$, $\chi_{m-s}$, $\chi_{s-s}\}$ are the various interaction strengths, relative to $k_{B}T$, associated to Van der Waals, excluded volume, and hydrophobic interactions [1]. Re-writing in terms of $\phi$, the monomer density is $\rho_{m}=\phi/v$ and the solvent density is $\rho_{s}=(1-\phi)/v$ so that

[TABLE]

Each particle, independent of identity, shares the same mean energy density; the total energy $\overline{E}$ is therefore $V\overline{\mathcal{E}}$. The quantity of interest is the change in energy $\Delta E$ due to mixing, which is $\overline{E}-\overline{E}_{s}-\overline{E}_{m}$, where $\overline{E}_{s}=n_{s}v\overline{\mathcal{E}}(\phi=0)$ is the the total energy of a fictitious system having the same number $n_{s}$ of solvent molecules but without any monomers, so that $\phi=0$; similarly, $\overline{E}_{m}=n_{m}v\overline{\mathcal{E}}(\phi=1)$ is a system of monomers alone. Introducing the total number $N=n_{s}+n_{m}$ of solvent molecules and monomers, the mixing energy $\Delta E$ is simply

[TABLE]

where $\chi[\equiv\chi_{m-s}-\chi_{m-m}/2-\chi_{s-s}/2]$ is the so-called Flory parameter [1]. If $\chi<0$ the interaction energy is minimized when $\phi=1/2$, corresponding to equal parts of solvent and polymer. Since this case occurs when $\chi_{m-s}<(\chi_{m-m}+\chi_{s-s})/2$, it describes a regime in which the energetic cost of monomer-solvent interactions is less the average cost of pure monomer-monomer and pure solvent-solvent interactions.

Deformations of the polymer network generally result in changes in the contact interactions between chains. At the mean-field level, these interactions are incorporated in the mixing energy $\Delta E$ through $\phi$ alone. For anisotropic deformations at fixed $\phi$, however, we will assume that interactions between polymers are somewhat less important and approximate the polymer network via a phantom chain model where the chain conformations are allowed to overlap one-another, leading to random-walk “ideal” polymers. We thus focus on single-chain deformations, implicitly assuming that this is the main contribution to the entropic cost of stretching the polymer network $\Delta S_{\rm net}$. Although some degree of realism is lost, the problem gains tractability whilst retaining the essential physics – the free energy cost of elastic deformations is simple to derive and has a form that reduces to the classical rubber elasticity model (see e.g., [29]) in the unswollen limit. With this assumption, the elastic free energy of the gel is approximated as proportional to the net conformational entropy change due to deforming $n_{ch}$ independent polymer chains. Thus, we require knowledge of (i) how deformations affect the conformational entropy of a single chain and (ii) how to determine the effect on an ensemble of many such chains. It should be noted, however, that this construction is limited in scope and fails to accurately capture the elastic free energy in the large shear-strain regime, where correlations between polymer fluctuations gain importance [30].

To begin, consider a single chain of length $L$ with one terminal end at position $\mathbf{r}$ and the other at the origin, as shown in figure 3(a). In order to determine the entropy of the chain, first note that within the phantom chain model, all random paths of fixed length $L$ have the same energy. Therefore, the entropy $S_{1}(\mathbf{r})$ of a single chain with end-to-end vector $\mathbf{r}$ is given by $S_{1}(\mathbf{r})=k_{B}{\rm ln}\,\Omega_{1}(\mathbf{r})$, where $\Omega_{1}(\mathbf{r})$ is the total number of microstates available to the chain. We can model chain conformations by considering a lattice model in which each lattice site has length $a$ and each component $r_{i}$ of the vector $\mathbf{r}$ can be expressed as $r_{i}/a=2n_{i}-N_{i}$, where $n_{i}$ is the displacement along the lattice along one axis and $N_{i}$ is the total number of steps taken along that axis. The number of microstates $\omega(r_{i})$ for this one dimensional random walk is given by $N_{i}!/(n_{i}!(N-n_{i})!)$. Then $\Omega_{1}(\mathbf{r})$ is simply the product $(V/v)\omega(r_{1})\omega(r_{2})\omega(r_{3})$, where $V/v$ is the number of possible locations for $\mathbf{r}=0$, and is given by

[TABLE]

which can be simplified in the limit of large $N_{i}$ by taking advantage of Stirling’s approximation, yielding

[TABLE]

Assuming that the polymer explores three-dimensional space isotropically, $N_{1}=N_{2}=N_{3}=N/3$, where $N$ is the total number of monomer units. Furthermore, we note that the polymer is most likely to be found with end-to-end distance $|\mathbf{r}|$ to be much smaller than its extreme maximum length $Na$. Expanding ${\rm ln}\,\Omega_{1}$ in powers of $|\mathbf{r}|/(Na)$, the leading order contribution is

[TABLE]

where $R_{0}=\sqrt{N}a$ [1]. The entropy $S_{1}(\mathbf{r})$ for a single chain is therefore given by

[TABLE]

where $S_{0}$ is a constant that depends on the number of monomers in each chain. Now consider a process that stretches the chain, resulting in a new end-to-end vector $\mathbf{r^{\prime}}$, and hence a new entropy $S^{\prime}(\mathbf{r^{\prime}})$. Writing the displaced vector as $r^{{}^{\prime}i}=\Lambda_{ij}r^{j}$, where $\Lambda_{ij}$ is a deformation matrix, the change in entropy $\Delta S_{1}(\mathbf{r})$ for a single chain is given by

[TABLE]

Next, consider a collection of $n_{ch}$ chains that are cross-linked to form a polymer network. The macrostate of these chains, prior to deformation, is specified by the collection of end-to-end vectors $\{\mathbf{r_{1}},\mathbf{r_{2}},\dots,\mathbf{r_{n_{ch}}}\}$, corresponding to vectors between cross-links as shown in figure 3(b). If this network undergoes affine deformation, as illustrated in figure 2, then all end-to-end vectors are transformed by the same deformation matrix $\Lambda$. Therefore, the change in entropy for this collection of chains, $\Delta S_{n_{ch}}$, is simply given by a sum over independent entropy contributions (1r), namely

[TABLE]

where we have assumed that all chains have the same length and thus the same value of $R_{0}$ and used the relation $V=V^{\prime}({\rm det}\,\Lambda)$ for affine deformations. This is the change in entropy for a given network, represented by the collection of end-to-end vectors. Since we restrict our attention to chemical gels, the network topology is set upon cross-linking, and the same collection of end-to-end vectors describes the polymer network for all processes. While each network has a distinct topology, we assume that (i) the polymer networks are sufficiently large that the same collection of end-to-end vectors is represented throughout every network, albeit in a possibly different arrangement (by the self-averaging property), and (ii) cross-linking occurs when a collection of polymers in solution are brought to a concentration where they overlap but not to the point where they would deform due to steric repulsion. Then the change in entropy $\Delta S_{\rm net}$ representative of a polymer network composed of $n_{ch}$ chains of fixed $R_{0}$ is found by averaging equation (1s) over the equilibrium values of $\mathbf{r_{1}}\dots\mathbf{r_{n_{ch}}}$ before deformation. In order to do this, note that the probability $P_{1}(\mathbf{r})$ that a single polymer will have an end-to-end vector $\mathbf{r}$ is given by $P_{1}(\mathbf{r})=\Omega_{1}(\mathbf{r})/Z_{1}$, where $Z_{1}=\int{\rm d}^{3}r\Omega_{1}(\mathbf{r})$. Note that in the phantom chain model, the energy is then equal to zero. Using the result that

[TABLE]

the change in entropy $\Delta S_{\rm net}\equiv\left<\Delta S_{n_{ch}}(\mathbf{r_{1}},\dots,\mathbf{r_{n_{ch}}})\right>$ due to deforming a polymer network is given by

[TABLE]

However, we still have not arrived at our final result for $\Delta S_{\rm net}$. Since chemically cross-linked chains share a common endpoint, some of the chain degrees of freedom must be eliminated [23, 31]. This will cause $\Delta S_{\rm net}$, as expressed in equation (1u), to decrease. To estimate this reduction, note that within the phantom chain assumption, the endpoint $\mathbf{r}$ of a chain is free to lie within any point in the volume $V$ of the gel, irrespective of the location of where the polymer is based. After deformation, $V\rightarrow({\rm det}\,\Lambda)V$ so the change in entropy due to the change in the volume of the gel that is accessible to the endpoint is given by $n_{ch}k_{B}{\rm ln}\,{\rm det}\,\Lambda$; this exactly cancels the last term in Eq. (1u). However, each cross-link between two chains, say chain $i$ and chain $i+1$, constrains the endpoint motion of the chains; as a result, there is a constraint function $f_{i}(\mathbf{r}_{i},\mathbf{r}_{i+1})=0$ for these two chains. Since there are $n_{c}=n_{ch}/2$ such constraints, there is an additional reduction of the total entropy by $(n_{ch}k_{B}/2){\rm ln}\,{\rm det}\,\Lambda$. Thus, the overall entropy change due to deformations of the polymer network is

[TABLE]

which is attributed to Flory and Wall [29]. Note $\Delta S_{\rm net}\propto n_{ch}$ and is independent of chain length.

We emphasize, however, that the argument for reduction in entropy due to cross-linking, as presented above, is somewhat flawed. In the seminal paper of Deam and Edwards [7], it was shown that this argument relies on the assumption that the cross-linked ends of the chains are free to explore the entire volume $V$ of the gel. However, the cross-linked ends of the chains, whilst able to undergo thermal motion, are localized to a much smaller volume $\omega$ when the gel is formed. Furthermore, whereas the volume $V$ of the gel depends on the affine deformation $\Lambda$, the volume $\omega$ of the localization is a much weaker function of $\Lambda$ owing to non-affine fluctuations of the cross-linked endpoints. Therefore, the ${\rm ln}\,{\rm det}\,\Lambda$ term in the Flory-Wall entropy (1v) is not completely justified. In addition, since the term can be re-written as ${\rm ln}(\phi_{0}/\phi)$, it only depends on the polymer volume fraction. Since the mixing of solvent and polymer result in a similar contribution to the total entropy, it is nevertheless difficult to assess the validity of the inclusion of this term in the Flory-Wall entropy.

There is an additional contribution to the entropy coming from the mixing of solvent and the polymer network: $\Delta S_{\rm mix}$. To estimate $\Delta S_{\rm mix}$, we model the space occupied by the gel by a lattice, as shown in figure 4, of $N$ sites, each occupied by either a solvent molecule or a monomer; because the system is densely filled, there are $n_{s}$ solvent molecules and $n_{m}=N-n_{s}$ monomers. By fixing the total number of monomer and solvent molecules, the entropy $S_{\rm latt}$ of arranging monomers and solvent into the lattice is given by $S_{\rm latt}=k_{B}{\rm ln}\,\Omega$, where $\Omega$ is the number of possible lattice arrangements. We must therefore count the number of ways that the lattice can be filled with solvent and monomers, where the monomers (i) are arranged into polymers that (ii) belong to a cross-linked network that spans space. We will start with the simple case of “free” monomers that are unassociated into larger polymer molecules, all able to explore space independently and recover the entropy of Bragg-Williams theory [32]. Subsequently, we will progressively introduce the necessary constraints by associating the monomers into polymers and then introducing the cross-linking constraints.

Consider a binary system, consisting of $n_{A}$ and $n_{B}$ particles of species ‘A’ and ‘B’, respectively; ‘A’ could, for example, represent solvent and ‘B’ could represent free monomers, as shown in figure 4(a). The total number of microstates of the lattice is given by

[TABLE]

so that the entropy, using Stirling’s approximation is

[TABLE]

Recalling that the volume fraction $\phi=n_{B}/N$,

[TABLE]

from which we find that the state of maximum entropy corresponds to $\phi=1/2$, which further corresponds to a mixed state composed equally of both species of particles.

We now associate monomers into polymer that are free to explore the entire space, whilst localizing individual monomers to much smaller volumes around the centers of mass of the polymers, as illustrated in figure 4(b). Let each polymer consist of $\mathcal{N}$ monomers such that $n_{p}=n_{m}/\mathcal{N}$ is the total number of polymers. Whilst the individual monomers have the adjacency condition, polymers are allowed full translational freedom on the lattice. Thus, the entropy of the localized monomers is negligible compared with the translational entropy of the polymers. The entropy $S_{\rm latt}$ is therefore dominated by the translational entropy of the solvent and the polymers such that

[TABLE]

To obtain the mixing entropy $\Delta S_{{\rm mix}}$, first define an entropy density $\mathcal{S}=S/V$. Following the definition of the mixing energy $\Delta E$, the mixing entropy is given by

[TABLE]

which corresponds to the Flory-Huggins result [33, 1] for polymer solutions.

Finally, we consider the case in which permanent cross-links are introduced, localizing polymers to small regions about the cross-link sites [see figure 4(c)]. In this case, the polymers have constraints that reach all the way to the sample boundary, resulting in rigidity. Therefore, the translational entropy of polymers is negligible compared with the entropy of the solvent. The result may be found by considering the limit of the Flory-Huggins theory for infinitely long polymer, i.e., taking $\mathcal{N}\rightarrow\infty$. The result is

[TABLE]

which is independent of network details [23, 1, 31].

The deformation free energy $\Delta F$ can finally be decomposed as

[TABLE]

where the elastic deformation free energy

[TABLE]

arises from the entropy change due to deformation of the polymer network, and the mixing free energy,

[TABLE]

is the net change in the free energy due to mixing polymer and solvent. In the last term of equation (1y), the constant $\phi_{0}$ corresponds to the volume fraction in the reference state of the gel, usually taken to be the volume fraction at which cross-linking is performed, or occasionally the volume fraction of a completely dry gel, namely $\phi_{0}=1$. Notice that both the elastic and mixing free-energies scale with the thermal energy $k_{B}T$—the only term that presents a non-linear scaling with temperature is the Flory parameter term since $\chi$ is a function of temperature $T$. It is therefore convenient to re-scale the total free energy $\Delta F$ by the thermal energy, i.e., $\Delta F/k_{B}T$, from which we find that the equilibrium state of polymer gels is determined by $\chi(T)$ alone.

A more careful and detailed look at the theory of gel elasticity confirms that the affine-deformation picture of classical rubber elasticity is inaccurate [8]. While the average cross-link positions in space undergo affine transformation under a homogeneous deformation of the gel at its boundaries, there are in fact large fluctuations in cross-link positions due to thermal motion as well as network inhomogeneities. In fact, these fluctuations are on the order of the mean cross-link spacing, which would melt ordinary solids, according to the Lindemann criterion, further highlighting the strangeness of these materials. Additionally, the separation of the total free energy into a contribution due to the network elasticity and a contribution due to solvent-polymer mixing ultimately fails due to these large fluctuations, which renormalize both contributions. Thus, while we will use the Flory-Rehner to illustrate the thermodynamics of polymer gels, it should be regarded as a semi-empirical model, that over-simplifies the true microscopic state of the gel.

2.3 Isotropic swelling

When allowed to equilibrate with a solvent bath, the amount of solvent in a polymer gel balances the osmotic pressure due to the thermal motion of the polymer network with the entropic cost of stretching this network. Changes in this equilibrium state can be brought about by changing the solvent quality, as characterized by the Flory parameter $\chi$. With the Flory-Rehner free energy (1y) in hand, let us determine the equilibrium volume fraction $\phi(T)$ in the case of an isotropic gel. Applying the volume constraint, viz., $\partial\Delta F/\partial\lambda=0$, we readily obtain that the deformation matrix $\Lambda$ is given by

[TABLE]

The free energy density $\Delta\mathcal{F}$ of the gel is therefore

[TABLE]

where $\nu_{0}\equiv n_{ch}/V_{0}$ is the density of chains in the reference state of the gel. The osmotic pressure follows from $\Delta F$ and is given by

[TABLE]

where we have taken the Flory parameter $\chi$ to be a function only of temperature $T$. Equation (1ac) is the Flory-Rehner equation of state relating the osmotic pressure $\Pi$ to the volume fraction $\phi$ and temperature (via $\chi$). Contours of variable volume fraction $\phi$ and Flory parameter $\chi$ at constant osmotic pressure $\Pi$ are shown in figure 5(a). For $\Pi=0$, the volume fraction $\phi$ increases with decreasing $\chi$, corresponding to a gel that is swollen (low $\phi$) for a good solvent ($\chi<0.5$) and that deswells as the solvent becomes poor. Positive values of osmotic pressure can be obtained through the addition of a solute to the surrounding solvent; equilibrium osmotic isobars for positive values of $\Pi$ are also shown in figure 5(a).

We can understand the behavior of the osmotic isobars in analogy with isobars from the Van der Waals equation of state

[TABLE]

which relates the pressure $P$ to the density $\rho$ of particles, as shown in figure 5(b). Note that for low density, these curves asymptotically approach their ideal gas form $\rho^{-1}\sim P^{-1}k_{B}T$. For sufficiently large positive pressure, the density $\rho$ decreases with increasing temperature, corresponding to an expanding gas. However, for low pressures and temperatures, the density is a multivalued function of temperature. In the case of the Van der Waals fluid, the emergence of the multivalued region is indicative of a loss of thermodynamic stability and the development of distinct liquid and gas phases. Therefore, we might expect that the Flory-Rehner theory of gels has a similar phase transition separating a distinct low $\phi$ swollen phase and a high $\phi$ deswollen phase. Such a phase transition indeed exists for gels, even though the Flory-Rehner equation requires a slight alteration to correctly capture it [34].

3 Phase transitions

While there is a useful analogy that may be drawn between the isotropic swelling of polymer gels, as modeled by the Flory-Rehner theory, and the thermal expansion of a fluid, as modeled by the Van der Waals equation of state, there is also a key difference. Examining the isobars in figure 5(b), one finds that for sufficiently low temperature and pressure $P$, the $\rho^{-1}=V/N$ versus $T$ plot is multi-valued. The value of pressure for which the well-defined, single-valued expansion curve becomes multi-valued is called the critical pressure. We highlight the situation in figure 6(a), which shows three different isobars: $P<P_{c}$, $P=P_{c}$, and $P>P_{c}$. For $P>P_{c}$, a fluid that is quasistatically heated from high-density (low $\rho^{-1}$) becomes lower-density (higher $\rho^{-1}$); this process is easily reversed upon cooling. For $P=P_{c}$, while the equilibrium heating and cooling paths remain the same, the change in density with temperature diverges at a certain critical temperature $T_{c}$. However, for $P<P_{c}$, a high-density fluid can be quasistatically heated so that the density traces the lower part of the curve shown in figure 6(a) until the curve folds back on itself. Heating beyond this transition temperature $T_{>}$ results in a discontinuous jump to a much lower density and the ensuing thermal expansion follows the upper branch of the isobar. Cooling the fluid from low density and high temperature, however, traces the upper branch, until the discontinuity is encountered at a lower transition temperature $T_{<}$. The low-density and high-density values of the fluid for $P<P_{c}$ distinguish separate fluid phases, which we recognize as the gas phase and the liquid phase, respectively. Thus, this appearance of (i) a discontinuous jump in fluid density that (ii) depends on the heating path is not a failure of the Van der Waals model but rather a successful description of a first-order phase transition.

Interestingly, experiments on certain polymer gels, notably pNIPAM, reveal similar discontinuous behavior in the equilibrium swelling curves [35, 36]. For low temperatures, when the polymer network is miscible in the solvent, the gel is swollen. Slowly increasing temperature increases the cost of polymer-solvent interaction, i.e., increases $\chi$, leading to gradual deswelling. Above a certain temperature, roughly $32^{\circ}{\rm C}$ for pNIPAM, the gel suddenly expels most of its solvent into the surrounding bath, reducing its volume by orders of magnitude, and becomes opaque. This discontinuity hints at a similar first-order phase transition of polymer gels and distinguishable swollen and deswollen phases. However, the osmotic pressure of the Flory-Rehner model, equation (1ac), does not exhibit the multi-valued behavior of the Van der Waals model. We will discuss the Erman-Flory extension of the Flory-Rehner model that allows for such a phase transition. However, we will first briefly discuss the theory of phase transitions and critical phenomena more broadly.

3.1 Preliminaries: general aspects of phase transitions

In order to understand phase transitions, let us first consider thermodynamic stability. While this discussion is generalizable, we will continue to use the example of the Van der Waals model of fluids. We plot the constitutive relation between pressure $P$ and inverse density $\rho^{-1}$ for fixed temperature in figure 6(b). Evidently, if we are able to fix the temperature $T$ of a fluid and specify the total number of particles $N$ and volume $V$, then equation (1ad) tells us the pressure of the fluid, assuming that the density $\rho$ of particles is uniform everywhere. Of course, at finite temperature, particles in the fluid undergo thermal fluctuations and the density $\rho$ varies in space and time. The microscopic length scale over which spatial variations in particle density can be resolved is the correlation length $\xi$. To connect to macroscopic physics, state functions like $\rho$ are found by coarse-graining, or averaging over many particles in a certain region of space, the size of which is set by a coarse-graining length scale $\ell$. By taking \ell{\;\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;}\xi, thermal fluctuations are averaged out and we may approximate the state functions of each coarse-grained region by their thermodynamic limit. For example, if we label a certain region by a position $\mathbf{x}$ so that the local coarse-grained density is $\rho(\mathbf{x})$ then the local pressure $P(\mathbf{x})$ can be approximated by the Van der Waals equation of state (1ad).

Now consider a disturbance to the fluid, such as a vibration or incident sound wave. The result of such an external influence is a spatial modulation in density, typically of a longer length scale than $\xi$. One coarse-grained region may have a slightly lower density of particles than its surroundings, which may have a slightly higher density of particles. Consulting figure 6(b), we find that for the higher temperature isotherms, pressure increases monotonically with density. Therefore, assuming that the fluid is maintained at fixed temperature, the higher density regions are at higher pressure and the lower density regions are at lower pressure. Subsequently, due to this pressure difference, particles will migrate from higher density to lower density in order to re-establish equilibrium. The higher density region and the lower density regions eventually settle to a uniform density, namely $N/V$. This resilience to perturbations is known as thermodynamic stability. Per this argument, the essence of Le Chatelier’s Principle, thermodynamic stability requires

[TABLE]

where $K$ is the bulk modulus, which is the inverse of the isothermal compressibility $\kappa_{T}$ [26]. As long as the temperature $T>T_{c}$, where $T_{c}$ is the critical temperature, the constitutive relation $P(\rho)$ obeys this stability requirement. However, at $T_{c}$, there is a critical pressure $P_{c}$ at which $\left(\partial P/\partial\rho\right)_{T}=0$. This critical point marks the loss of thermodynamic stability and the onset of different physics. The zero value of the bulk modulus $K$, or diverging compressibility $\kappa_{T}$, at the critical point is one example of the critical phenomena that one encounters.

Saving the discussion of critical phenomena for later, let us address the consequences of thermodynamic instability. If $T<T_{c}$ and $P<P_{c}$ then there are certain values of $\rho$ where $\left(\partial_{P}/\partial\rho\right)_{T}<0$ so the bulk modulus $K$ is negative. Higher density regions are at lower pressure than the lower density region, so there is mass flow away from low density regions. As a result, density fluctuations grow and the ultimate fate of such a fluid is phase-separation into regions of low density and regions of high density. However, density variations cannot grow forever: eventually, the high density and low density regions leave the unstable region of figure 6(b) and enter stable regions of positive compressibility. The difference between the higher and lower densities grows with the distance of the fluid from the critical point, characterized by a reduced temperature $t\equiv(T-T_{c})/T_{c}<0$; for large enough values of $|t|$, these two densities describe well-defined, distinguishable phases. Eventually, these two fluid phases attain a phase-coexistent equilibrium.

Instead of using the Van der Waals equation of state to determine the pressure at a given density, now consider it as a way to determine density at a given pressure. Much like the $\rho^{-1}(T)$ isobars in figure 5(b), the $\rho^{-1}(P)$ isotherms are multi-valued graphs for $T<T_{c}$. Thus, a high-density, liquid-phase fluid at $T<T_{c}$ undergoes a first-order phase transition to a low-density, gas-phase fluid for sufficiently low pressure. However, we have shown that there should also be cases where the fluid is in a phase-coexistent equilibrium between liquid and gas phases. In order for these phases to coexist in equilibrium, they must (i) have the same temperature $T$, (ii) the same pressure $P$, and (iii) the same chemical potential $\mu$; that is, they must be in thermal, mechanical, and chemical equilibrium. While $T$ and $P$ are specified, $\mu$ must be determined. We can take advantage of the Gibbs-Duhem relation $N\,{\rm d}\mu=-S\,{\rm d}T+V\,{\rm d}P$, which, at constant temperature, can be expressed as ${\rm d}\mu=\rho^{-1}\,{\rm d}P$ [26]. Therefore, chemical equilibrium is realized when

[TABLE]

where $\mathcal{P}_{1}$ is the path along the isotherm, shown in figure 7(a), that connects the point $(\rho^{-1}_{\ell},P^{*})$ to $(\rho^{-1}_{g},P^{*})$, where $\rho_{\ell}$ and $\rho_{g}$ are the respective densities of the liquid and gas phases. Joining these two points by a constant-pressure line $\mathcal{P}_{2}$, we can define a loop $\mathcal{P}=\mathcal{P}_{1}\cup\mathcal{P}_{2}$ as the union of these two paths; this is called a “Van der Waals loop.” Integrating equation (1af) by parts,

[TABLE]

we therefore find that coexistent phases are in equilibrium when the net area enclosed by the Van der Waals loop $\mathcal{P}$ is 0. The line $\mathcal{P}_{2}$ that joins coexisting densities $\rho_{\ell}$ and $\rho_{g}$ gives the equilibrium pressure and replaces the $S$-shaped curve $\mathcal{P}_{1}$. This “Maxwell construction” corrects multivalued isotherms in the Van der Waals equation of state [26].

Furthermore, in the two-dimensional $P\;vs.\;T$ phase diagram of fluids, we can identify a one-dimensional locus of coexistent equilibria consisting of the single pressure $P$ that yields coexistence for each isotherm $T$. This coexistence curve in the phase diagram terminates at the critical point $(T_{c},P_{c})$. A fluid may be brought around the critical point without passing through the coexistence curve via appropriate temperature and pressure change protocols. Whilst points on the coexistence curve describe coexistence between gas and liquid phases, points immediately to one side or the other are in the single phase region. Passage through the coexistence curve results in a first-order phase transition, a discontinuous jump between high-density and low-density fluids. Therefore, the only way to distinguish between gas and liquid phases of fluids is to pass through the coexistence curve. In fact, the ability for a system to support coexisting densities in equilibrium is the defining feature of distinct phases that are separated by a first-order phase transition.

3.2 The common tangent construction in phase-separating systems

The phase diagram 7(b) shows a phase-coexistent region for a particular set of temperatures and pressures. To land on this curve, one requires precise control over temperature and pressure, suggesting that coexistent phases are rarely realized. As it turns out, one can readily achieve phase coexistence at constant temperature, volume, and number of particles $(T,V,N)$.

Rather than working with the pressure-density equation of state, consider the Helmholtz free energy $F(T,V,N)$. Since $F$ is a thermodynamic potential, it is a homogeneous first-order function in $V$ and $N$. Therefore, we can write $F$ in terms of its density $\mathcal{F}$ as $F(T,V,N)=N\mathcal{F}(T,\rho^{-1})$, defined in terms of $N$. If the free energy density $\mathcal{F}$ describes a thermodynamically stable system then the requirement for positive isothermal compressibility is satisfied for $(\partial^{2}\mathcal{F}/\partial(\rho^{-1})^{2})_{T}>0$ and there is a single free energy minimum $(\rho^{-1})^{*}$ for fixed values of $T$, as illustrated in the inset in figure 8. However, for $T<T_{c}$, the free energy has a region of negative compressibility where $(\partial^{2}\mathcal{F}/\partial(\rho^{-1})^{2})_{T}<0$ and $\mathcal{F}$ is concave when plotted against $\rho^{-1}$. In this case, the free energy supports two local minima, separated by a local maximum, as shown in figure 8. Lacking a volume constraint, the system will seek to minimize the free energy so a local minimum that is not the global free energy minimum is considered meta-stable: eventually, given sufficient time, thermal fluctuations will drive the system to the global free energy minimum. For example, it is possible to “superheat” a homogeneous liquid-phase fluid above the transition temperature for phase coexistence, keeping the fluid in its liquid phase, a metastable equilibrium, for a prolonged period of time. To do this requires careful preparation, removing any possible nucleation sites for the gas phase from the liquid at, for example, small pockets of trapped gas. The introduction of a nucleation site, e.g., via disturbing the fluid, lowers the free energy barrier locally and allows a portion of the liquid to transition to the gas phase. Without the introduction of a nucleation site from external influence, the superheated liquid nevertheless has a finite, albeit much longer, lifespan as a significantly large density fluctuation, driven by thermal fluctuations, will eventually provide a suitable nucleation site. Since the gas-phase is of lower density than the liquid-phase, it occupies a larger volume than the same mass of liquid-phase fluid. However, if the total volume of the fluid is constrained to remain constant, then even though the free energy associated with the gas is lower than that of the liquid, not all of the liquid can freely transition to the gas-phase. Instead, a portion of the liquid can transition to the gas-phase, at the expense of increasing the density of the liquid phase, achieving a phase-coexistent equilibrium.

In order to determine the conditions for equilibrium phase coexistence at constant volume, recall that the equilibrium state minimizes the global free energy $F$. Let $\rho_{\ell}$ and $\rho_{g}$ be the densities of the liquid and gas phases, consisting of $N_{\ell}$ and $N_{g}$ particles, respectively. By conservation of mass, the total number of particles in the container remains unchanged: $N_{\ell}+N_{g}=N$. Therefore, we can define a fraction $f\equiv N_{g}/N$ of particles that are in the gas phase; by mass conservation, the fraction of particles that are in the liquid phase is $(1-f)$. Furthermore, the total volume of particles in the liquid and gas phases are given by $V_{\ell}=\rho_{\ell}^{-1}N_{\ell}$ and $V_{g}=\rho_{g}^{-1}N_{g}$. Volume conservation requires that $V_{\ell}+V_{g}=V$; dividing by $N$, this conservation may be expressed as $f\rho_{g}^{-1}+(1-f)\rho_{\ell}^{-1}=\rho^{-1}$, where $\rho\equiv N/V$ is the nominal density of the fluid, as if it were in a single phase. Therefore, the fraction $f$ of particles in the gas-phase is given in terms of the equilibrium densities by

[TABLE]

a result known as the Lever Rule.

The phase-coexistent equilibrium is determined by minimizing the total Helmholtz free energy $F$,

[TABLE]

which is the sum of free energy contributions from each phase, where $\mathcal{F}_{g}=F(\rho_{g}^{-1})/N_{g}$ and $\mathcal{F}_{\ell}=F(\rho_{\ell}^{-1})/N_{\ell}$, subject to the volume constraint, enforced by a Lagrange multiplier $\lambda$. The minimization condition ${\rm d}F=0$ requires that partial derivatives of the total free energy with respect to $\rho_{g}^{-1}$, $\rho_{\ell}^{-1}$, $f$, and $\lambda$ are all equal to 0. This results in the Lever Rule (1ah) along with three additional equilibrium equations, namely

[TABLE]

Consider an example free energy density $\mathcal{F}(\rho^{-1})$ at low enough temperature that it has two local minima, as plotted in figure 8. The first two equilibrium equations (1aja,1ajb) show that the slopes of the free energy density $\mathcal{F}(\rho^{-1})$ at $\rho_{g}^{-1}$ and $\rho_{\ell}^{-1}$ are the same. The third equation (1ajc), after taking into account that $\lambda$ is the slope of the lines that are tangent to the free energy at the equilibrium densities, shows that the two tangent lines overlap. Thus, the equilibrium densities satisfy a common tangent construction: the values $\rho_{g}^{-1}$ and $\rho_{\ell}^{-1}$ can be found graphically by drawing a straight line that is tangent to $\mathcal{F}(\rho^{-1})$ at two points. As long as the free energy has two stable equilibria that are separated by an unstable equilibrium, this construction yields unique values for the equilibrium densities, and thus the fraction $f$ via the Lever Rule (1ah). Importantly, appreciate how the equilibrium densities are not at the local minima of $\mathcal{F}$. Instead, $\rho_{g}^{-1}$ and $\rho_{\ell}^{-1}$ are close to those minima, as can be seen in figure 8.

There is a physical rationale behind the equilibrium equations. Note that the pressure $P=-(\partial F/\partial V)_{T,N}=-(\partial\mathcal{F}/\partial\rho^{-1})_{T}$ so that the first two equations (1aja,1ajb) yields the interpretation of $\lambda$, a generalized force that maintains the volume constraint, as the pressure $P$ of the two phases, which are in mechanical equilibrium. To interpret the third equation (1ajc), note that $\mathcal{F}-\lambda\rho^{-1}=(F-PV)/N=G/N$, where $G$ is the Gibbs free energy. Therefore, we have that $G_{g}/N_{g}=G_{\ell}/N_{\ell}$, a balance of the Gibbs free energy density for each phase. Noting that $G=\mu N$, this is also a balance between chemical potentials $\mu_{g}$ and $\mu_{\ell}$, so the two phases are in chemical equilibrium. This result details a robust way to achieve equilibrium phase coexistence: as long as the temperature $T$ is low enough for the system to be thermodynamically unstable for certain values of the state functions, then there is equilibrium between coexistent phases at constant $N$ and $V$.

3.3 A note on critical phenomena

While we will not linger on the rich subject of physics near a critical point, a discussion of phase transitions requires at least a cursory mention of critical phenomena. Until this point, we have focused on equilibrium thermodynamics in the vicinity of the coexistence curve. Crossing this coexistence curve results in a first-order phase transition, which allows us to distinguish phases, such as the gas and liquid phase of a fluid. Sitting on the coexistence curve, the system is not a single homogeneous phase but rather an admixture of two phases that are in equilibrium with each other. The critical point is the endpoint of this coexistence curve and thus marks the onset of distinction between the two phases; equivalently, coexistence breaks down as this end of the curve is reached and the two phases lose distinction.

To study thermodynamics close to a critical point, it is helpful to define an order parameter $\varphi$ that is zero when only a single phase exists (off of the coexistence curve) and is nonzero when multiple phases exist. For a fluid, a choice of order parameter is the reduced density, namely

[TABLE]

where $\rho_{c}$ is the average density of the fluid at the critical point. For the Van der Waals equation of state (1ad), the critical temperature is given by $T_{c}=8a\rho_{0}/27$ and the critical pressure is $a\rho_{0}^{2}/27$ so the critical density is $\rho_{c}=\rho_{0}/3$. For small negative values of the reduced temperature $t=(T-T_{c})/T_{c}$ along the critical isochore, where the fluid at the critical density $\rho_{c}$ is in unstable equilibrium, there are two new locally stable equilibria, one with $\rho>\rho_{c}$ and another with $\rho<\rho_{c}$; there is a positive and a negative value of $\varphi$. Defining a reduced pressure $p=(P-P_{c})/P_{c}$ and expanding the Van der Waals equation of state (1ad) yields a linear leading-order dependence of $p$ with $\varphi$ for small $\varphi$, that is $p\approx r\varphi$, where the coefficient $r\propto t$. Noting that this implies that $\partial p/\partial\varphi\propto t$, we find stability via the Le Chatelier principle for $t>0$ and instability for $t<0$. To recover the appearance of two new locally stable equilibria, there needs to be a dependence on $\varphi$ added to $p$ that yields positive bulk modulus for sufficiently large values of $|\varphi|$. The simplest addition that stabilizes the reduced equation of state is a cubic $\varphi^{3}$ dependence, yielding $p=r\varphi+u\varphi^{3}$, where $u>0$. We can integrate the pressure to find an approximate form of the free energy density $\mathcal{F}$ near the critical point, namely

[TABLE]

One immediate consequence of this approximation is that the equilibrium values of the new phases close to the critical point, i.e. small $|t|$, are given by

[TABLE]

which are symmetric about the unstable equilibrium $\varphi=0$. Note that this symmetry between the two phases holds only close to the critical point: further away, odd powers of $\varphi$ appear in the free energy. Still, close to the critical point, we find that the separation in density between the liquid and gas phases, $(\rho_{\ell}-\rho_{g})/\rho_{c}=2|\varphi^{*}|$, where $\varphi^{*}$ is a minimum of $\mathcal{F}$, scales with $|t|^{\beta}$ with $\beta=1/2$ along the critical isochore. The exponent $\beta$ is one example of a critical exponent. There are a variety of critical exponents for different thermodynamic quantities, such as the isothermal compressibility $\kappa_{T}\sim|t|^{-\gamma}$. Since $\kappa_{T}^{-1}\propto{\rm d}^{2}\mathcal{F}/{\rm d}\varphi^{2}$, we find that $\gamma=1$ along the critical isochore. These critical exponents are used to characterize the behavior of a system near a critical point.

The problem with the above analysis is that it does not yield good predictions for critical exponents that are measured in the lab. The measured value of $\beta$ is actually close to $1/3$, but does not seem to be a rational number [37]. Even though the Van der Waals equation of state works well for describing equilibrium physics of liquid and gas phases of fluids for much of the phase diagram, it seems to fail near the critical point. As it turns out, the root of the problem lies in the assumption that thermal fluctuations are not important. The key assumption was that thermal fluctuations, characterized by a correlation length $\xi$, are important only over small length-scales. To find the large length-scale physics, recall that in defining $\rho$ and other state functions as descriptions of the macroscopic state of a system, there was a coarse-graining length scale $\ell$ introduced over which the microscopic details, such as thermal fluctuations, were averaged over. This coarse-graining length $\ell$ was taken to be at least as large as the correlation length $\xi$. Since fluctuations are ignored in the resulting description of the thermodynamics, this is referred to as a mean-field theory. The predicted critical exponents are mean-field critical exponents, which, owing to the simple structure of mean-field theories, are always rational numbers.

In order to correct mean-field theory, we need to properly incorporate this coarse-graining length scale and investigate corrections due to thermal fluctuations in the order parameter $\varphi$. To do this, we can model a fluctuating order parameter near the critical point by a model free energy

[TABLE]

where the coefficient $c$ sets the energy cost of spatial variations in $\varphi$, and $d$ is the dimensionality of the space in which properties of the material vary. Note that this is completely phenomenological: the presence of the gradient term $|\nabla\varphi|^{2}$ simply provides a positive free energy cost for spatial variation. This term is necessary for the development of the coarse-grained model of the fluid as it describes a lower cutoff length for the wavelength of spatial fluctuations in the order parameter field $\varphi$. To see this, let $\lambda$ represent the wavelength of a fluctuation in $\varphi$. Then there is an energetic cost of this fluctuation that scales as $c|\varphi|^{2}/\lambda^{2}$ so that as the length scale $\lambda$ of the spatial variation in $\varphi$ decreases in size, the cost of this fluctuation grows.

To see how adjusting $c$ affects the length $\xi$ over which fluctuations $\delta\varphi(\mathbf{x})$ of $\varphi(\mathbf{x})$ about the equilibrium $\varphi^{*}\equiv\pm\sqrt{|r|/u}$ are correlated in space, we can determine the functional form of the fluctuation correlations, namely $\left<\delta\varphi(\mathbf{x})\delta\varphi(0)\right>$. To leading order in fluctuations $\delta\varphi$, the change $\delta F$ in the free energy is given by

[TABLE]

where $F_{0}=-r^{2}V/(4u)$ is the free energy corresponding to the homogeneous equilibrium $\varphi^{*}$. The probability that a particular fluctuation $\delta\varphi(\mathbf{x})$ is weighted by the Boltzmann factor $\exp(-\beta\delta F[\delta\varphi])$, where $\beta=(k_{B}T)^{-1}$. Therefore, the fluctuation correlations are determined by

[TABLE]

where $\int[{\rm d}\delta\varphi]$ represents a sum over all possible fluctuations $\delta\varphi$ of the order parameter field. Integrating the free energy fluctuation (1ajao) by parts to yield

[TABLE]

we recognize that the functional integrals in (1ajap) have a Gaussian form and are thus simple to evaluate. The fluctuation correlations are given by

[TABLE]

and thus satisfy the Green’s function equation [32]

[TABLE]

where $\delta(\mathbf{x})$ is the Dirac delta function in $d$ dimensions. Therefore, the position dependence of the fluctuation correlations is given by

[TABLE]

where $\xi\equiv\sqrt{c/(2|r|)}$ defines the correlation length between thermal fluctuations in the coarse-grained field $\varphi$ [37]. As long as $|r|>0$, the isothermal compressibility $\kappa_{T}\sim r^{-1}$ is finite, and we can always therefore coarse-grain to a length-scale $\ell$ larger than the correlation length $\xi$. However, as $r\rightarrow 0^{-}$ on the approach to the critical point, this length-scale is ill-defined because $\xi$ diverges! Therefore, fluctuations cannot be ignored and mean-field theory is destined to fail. Indeed, this is confirmed in experiment via the phenomenon of “critical opalescence” [38]. As an otherwise transparent fluid, such as water, approaches the critical point, it turns opaque. Whereas normally, the correlation length $\xi$ is shorter than the wavelength of visible light, on the approach to the critical point, it lengthens to the point that thermal fluctuations in the density of the fluid can scatter light. The color of the fluid is a milky white, revealing that all visible wavelengths are scattered, so that fluctuations exist at many length-scales concurrently. Furthermore, this opacity lingers even as $\xi$ increases in length closer to the critical point, confirming that fluctuations at visible wavelengths remain, even as $\xi$ moves into the infrared and beyond, eventually stopping at the macroscopic length-scale $L\sim V^{1/3}$ of the container. Essentially, near the critical point, the fluctuations become scale-free.

Interestingly, while the mean-field theory predicts one set of critical exponents, in reality, critical exponents can vary from system-to-system. However, there are certain, seemingly unrelated, systems that share sets of critical exponents. For example, fluids, ferromagnets, and binary alloys all have approximately the same critical exponents [37]. This commonality of critical exponents amongst diverse systems means that their critical behavior is similar, even though the microscopic physics at play is very different, a phenomenon known as universality. Universality amongst systems is due to symmetry rather than microscopic physics. The order parameter $\varphi$ that we introduced for fluids represents a density difference. It works just as well for ferromagnets, which are described by magnetic dipoles that either point up or down; here, positive values of $\varphi$ correspond to an average magnetic dipole moment that is up and negative represents an average that is down. Similarly, for binary mixtures consisting of species labeled $A$ and $B$, $\varphi$ represents the difference in densities of species $A$ and species $B$. Regardless of the underlying microscopic physics, the form of the free energy at the critical point is identical, and the result is identical critical exponents, even when fluctuations are accounted for. Systems represented by other types of order parameters typically lie in other universality classes. The universality class of fluids, ferromagnets, and binary alloys is the three-dimensional Ising model, owing to the discrete $+/-$ symmetry of the free energy $\mathcal{F}$, namely, $\mathcal{F}(-\varphi)=\mathcal{F}(+\varphi)$. If, for example, $\varphi$ was a complex order parameter instead of a real scalar and if $\mathcal{F}$ was invariant under continuous transformation of the form $e^{i\theta}\varphi$, the corresponding critical phenomena would fall into the XY model universality class. For example, many systems with a polar order parameter $\bm{\sigma}$ that exhibit continuous rotational symmetry, such as superfluids, certain superconductors, and hexatic liquid crystals, lie in the universality class described by the XY model. The predictive power of the dimensionality of space, the dimensionality of the order parameter, and the symmetries of the system allow useful models of the behavior near the critical point via Landau theory, where a simple free energy, such as (1ajan), is constructed based on these considerations alone [32, 37].

3.4 Swelling-deswelling phase transition in polymer gels

Whilst many polymer gels undergo continuous changes in their polymer volume fraction due to changes in solvent conditions, e.g., via changes in temperature, certain gels exhibit a seemingly discontinuous change in volume fraction $\phi$, jumping between a low-$\phi$ swollen state to a high-$\phi$ deswollen state. As we have illustrated with our discussion about fluids, a discontinuous change in density in response to changing other state functions, e.g., temperature and pressure, indicates a first-order phase transition between a low-density and a high-density phase. For fluids modeled by the Van der Waals equation of state, these are the gas and liquid phases, respectively. Furthermore, little distinction between gas and liquid phases can be seen microscopically—unlike crystalline phases, there is no broken symmetry that distinguishes the two phases. The only sure way to distinguish these two phases is passage through a coexistence curve in the phase diagram that either crosses through the discontinuous transition or ends in a state of equilibrium phase coexistence. Therefore, we are led to conclude that polymer gels can have distinguishable swollen and deswollen phases. However, unlike in the Van der Waals model of fluids, the Flory-Rehner model of polymer gels does not predict discontinuous change in $\phi$ for physically reasonable parameters. Furthermore, within the formulation of the Flory-Rehner model for the osmotic pressure that has been presented thus far, not all values of $\chi$ yield a corresponding equilibrium value of $\phi$ when $\Pi$ is negative. One way of achieving $\Pi<0$ is by applying a mechanical pressure to the boundary of the gel, leading to solvent flow out of the gel via “reverse osmosis.” Interestingly, negative osmotic pressure states are typically thermodynamically unstable, favoring de-mixing of a solution into pure solute and pure solvent [39], i.e., phase-separation. This lack of general applicability suggests that equation (1ac) is an incomplete equation of state.

The Flory-Rehner model describes a rather simple picture of polymer gels in which the osmotic pressure $\Pi$ is expressed as two separate contributions: $\Pi_{\rm mix}$, which is due to the thermal motion of polymers amongst solvent molecules, and $\Pi_{\rm el}$, which is due to the elasticity of the polymer network. For ionic gels, thermal motion of free counterions contribute $\Pi_{\rm ion}$ to the osmotic pressure. This addition, which can be simply approximated as an ideal gas of counterions within the gel, is enough to theoretically obtain a discontinuous transition in the context of the Flory-Rehner model [40]. Ionic gels are indeed known to exhibit a discontinuous transition between swollen and deswollen phases. However, some neutral gels, such as pNIPAM, can also undergo a discontinuous transition yet do not have another obvious osmotic pressure contribution that is not captured within the Flory-Rehner model (1ac). Hence, as we have already emphasized, the Flory-Rehner model should be regarded as semi-empirical. This, in part, is due to the important role of thermodynamic fluctuations as well as static inhomogeneities in the polymer network. In the presence of a poor solvent, it has been shown [8] that, beyond a straightforward renormalization of the elastic and osmotic contributions, the presence of network inhomogeneities can lead to phase-separation, either at high wavenumber (microphase separation) or at low wavenumber (macrophase separation).

It is possible to extend the Flory-Rehner model such that it describes a phase transition. This is accomplished by altering the mixing energy between polymer and solvent molecules, which is controlled by the Flory parameter $\chi$. This term describes a two-body mean-field interaction between polymer and solvent molecules. To see this, expand the mixing contribution $\Pi_{\rm mix}$ to the osmotic pressure in powers of the polymer volume fraction $\phi$:

[TABLE]

where the sum is the remainder of the power series expansion of ${\rm ln}(1-\phi)$. This expansion has the form of a virial expansion of the pressure $P$ of a fluid in terms of its density $\rho$, namely

[TABLE]

where the leading order term is the ideal gas contribution and the higher order terms are corrections due to interactions between particles, which become important with increasing density $\rho$ [41]. In particular, the virial coefficient $b_{2}$ captures the effect of two-body interactions. For the Van der Waals equation of state, $b_{2}=\rho_{0}^{-1}-a(k_{B}T)^{-1}$, which is positive for sufficiently high temperatures, meaning that two-body interactions contribute an additional pressure to the independent-particle ideal gas term. This is much like the low-$\chi$ regime of polymer gels, which corresponds to the swollen phase. However for low temperatures, the two-body terms contributes a negative pressure, which drives particles to condense to a liquid phase, much as the polymer gel deswells for high-$\chi$. Note that at $\chi=1/2$, the two-body term disappears in the Flory-Rehner model, describing a $\vartheta$-solvent; this is analogous to the Boyle temperature of the Van der Waals model $k_{B}T_{b}\equiv\rho_{0}/a$, for which $b_{2}=0$.

Whereas the virial expansion for the Van der Waals model has two parameters, $\rho_{0}$ and $a$, the virial expansion for the Flory-Rehner model only has one, $\chi$, which is taken to be independent of polymer volume fraction $\phi$. However, measurements of $\chi$ have shown a nonlinear dependence on $\phi$ [42, 43, 44, 45, 46, 47, 48]. In general, the Flory parameter is a function of $\phi$ and can be expanded as a power series, namely

[TABLE]

and yields a more general virial expansion for the mixing contribution,

[TABLE]

Using this expansion, Erman and Flory [34] have shown that a discontinuous transition as well as a critical point can be recovered by tuning $\chi_{1}$ and $\chi_{2}$, and ignoring all other terms, i.e., $\chi_{m>2}\equiv 0$. In particular, acceptable fits to experimental swelling data can be found by fixing $\chi_{2}>1/3$ and varying $\chi_{1}$ with solvent quality, that is, taking $\chi_{1}$ to be a function of temperature only. It is important to emphasize that this expansion is purely phenomenological and does not assign specific microscopic meaning to the values of $\chi_{m\geq 2}$ [49]. In this phenomenological model of polymer gels, there are now two independent parameters, $\chi_{1}$ and $\chi_{2}$, in the virial expansion of the osmotic pressure, much like the two parameters of the Van der Waals model.

Example osmotic isobars $(\Pi={\rm const.})$ are shown in figure 9(a) and curves of constant $\chi_{1}$, corresponding to isotherms, are shown in figure 9(b). Much like the analogous processes shown for the Van der Waals model in figure (6), there are continuous and discontinuous swelling processes, along with an identifiable critical point. This critical point occurs at $\chi_{1}=\chi_{c}$ and $\Pi=\Pi_{c}$, determined by the condition that the osmotic bulk modulus $K$ vanishes, i.e.,

[TABLE]

and that the osmotic pressure isotherm is flat at the critical point, i.e.,

[TABLE]

Note that these conditions are equivalent to requiring $(\partial\Pi/\partial\phi)_{\chi_{c}}=0$ and $(\partial^{2}\Pi/\partial\phi^{2})_{\chi_{c}}=0$. Indeed, measurements of the bulk modulus $K$ near the phase transition show a dramatic softening when compared with the shear modulus $\mu$ of the gel [49, 50]. The breakdown of thermodynamic stability can be rectified by the existence of equilibrium phase coexistence between swollen and deswollen regions at some constant value of $\Pi$ by the Maxwell construction, namely

[TABLE]

where $\mathcal{P}_{1}$ is portion of the $S$-shaped curve beginning and ending at the equilibrium value of $\Pi$ in the thermodynamically stable region, as shown in figure 9(b). Therefore, much like the phase diagram predicted by the Van der Waals model, there is a similar phase diagram for polymer gels described via the Erman-Flory model, with a coexistence curve with a terminal critical point, as shown in figure 10. The phase diagram of polymer gel swelling is similar to that of a fluid modeled by the Van der Waals equation of state. Note that the coexistence curve is to the right of the critical point, whereas the curve in figure 7(b) is to the left. This is because the gas-like low-density phase, the swollen phase, occurs at low values of $\chi_{1}$, whereas the liquid-like high-density phase, the deswollen phase, occurs at higher values of $\chi_{1}$. To understand the negative slope of the coexistence curve, we turn to the Clapeyron relation [26], which relates this slope to the discontinuity change in volume and entropy that occurs when crossing the curve. For fluids, the positive slope ${\rm d}P/{\rm d}T>0$ shown in figure 7(b) tells us that the increase in volume per particle that occurs when a liquid evaporates and becomes a gas accompanies a corresponding increase in entropy. Conversely, the negative slope ${\rm d}\Pi/{\rm d}\chi_{1}<0$ for the gel coexistence curve indicates that there is a decrease in entropy as the gel passes from the deswollen phase to the swollen phase. This is to be expected of elastomeric materials in general: an increase in volume stretches polymer chains, decreasing their configurational entropy. Indeed, if one stretches a rubber band, decrease in the entropy of the polymer chains necessitates a passage of heat from the band into its surroundings, making it momentarily feel warm.

It is instructive to examine the behavior of the free energy near the coexistence curve and the critical point. However, the Erman-Flory virial expansion has to first be incorporated into the mixing free energy $\Delta F_{\rm mix}$. The power series expansion of the Flory parameter $\chi$ (1ajaw) cannot be directly substituted into the mixing free energy as the calculated osmotic pressure $\Pi_{\rm mix}$ is inconsistent with that given in equation (1ajax). Instead, start with the Erman-Flory mixing osmotic pressure in (1ajax) and integrate the relation $\Pi_{\rm mix}=\phi^{2}(\partial(\Delta\mathcal{F}_{\rm mix}/\phi)/\partial\phi)$ to find the mixing free energy density $\Delta\mathcal{F}_{\rm mix}$, up to an integration constant. Since we require that the mixing free energy disappears when the gel is either purely polymer, $\phi=1$, or in the limit where it is infinitely dilute, $\phi\rightarrow 0$, the integration constant is fixed, yielding

[TABLE]

where the original form of the mixing free energy (1aa) is recovered if only $\chi_{1}$ is retained. With this alteration to the total free energy $\Delta F$, there are values of $\chi_{1}$ and $\chi_{2}$ that cause $\Delta F$ to be a non-convex function of the polymer volume fraction $\phi$. This has the implication that the osmotic equilibrium

[TABLE]

may be satisfied for multiple values of $\phi$. Transforming to an analogue of the Gibbs free energy $G=\Delta F+\Pi V$, where $V=v\phi^{-1}$, the osmotic equilibrium condition is

[TABLE]

so that coexistent equilibrium volume fractions correspond to local minimia of $G(\phi)$, as illustrated in figure 10. Plots of $G(\phi)$ are shown for sample values of $(\chi_{1},\Pi)$ on the phase diagram 10; note that the inverse polymer volume fraction $\phi^{-1}$ is plotted on a logarithmic scale, reflecting the large scale of the volume changes that occur. At the critical point $(\chi_{c},\Pi_{c})$, the free energy minimum is broad, reflecting vanishing curvature $(\partial^{2}G/\partial(\phi^{-1})^{2})_{\chi_{1},\Pi}=0$, i.e., diverging isothermal compressibility. Along the coexistence curve, two degenerate minima of $G$ emerge, corresponding to the two coexistent phases. As predicted from Landau theory, via the model free energy (1ajan), close to the critical point, these two minima emerge symmetrically from the equilibrium value of $\phi^{-1}$ at the critical point and their separation grows continuously as the parameters $(\chi_{1},\Pi)$ are tuned along the coexistence curve; this is reflective of a continuous phase transition at the critical point. Instead, if $(\chi_{1},\Pi)$ are tuned transverse to the coexistence curve, the degeneracy of the free energy minima is lifted and there is an absolute minimum either for low $\phi^{-1}$ (deswollen phase) or for high $\phi^{-1}$ (swollen phase). In this case, as the coexistence curve is reached, the free energy difference between the two minima goes to 0, marking the cross-over between the two phases. This cross-over is a discontinuous jump between values of $\phi^{-1}$ that mark the local minima of $G$, indicating a first-order phase transition. However, physically, the gel may not immediately switch to the new absolute minimum. For example, if the gel is brought from the swollen phase across the coexistence curve to the deswollen phase, then the swollen phase still has a local free energy minimum; it is metastable. There is a free energy barrier for the gel to leave this metastable equilibrium and attain the global free energy minimum. Given sufficient time, thermal fluctuations will cause the gel to surmount this barrier. However, practically, the gel remains in the swollen phase until $(\chi_{1},\Pi)$ is tuned sufficiently far from the coexistence curve so that the free energy barrier disappears. Therefore, the observed values $(\chi_{1}^{*},\Pi^{*})$ at which the transition occurs are not typically on the coexistence curve, but rather to the right of it. Similarly, for the reverse process of starting from the deswollen phase and passing to the swollen phase, the transition typically occurs to the left of the coexistence curve. These states that linger past the phase coexistence transition are the polymer gel analogies of the superheated liquid and supercooled gas states of a fluid. This explains and generalizes the hysteresis that is shown in the osmotic isobars of figure 9.

3.5 Fluctuations and criticality in polymer gels

So far, we have discussed how the swelling thermodynamics of polymer gels share features in common with fluids. In discussing isotropic gels that undergo homogeneous changes in state, we have ignored the rigidity of these materials. At the critical point and along the coexistence curve, however, the gel is inhomogeneous, either undergoing large thermal fluctuations or possessing two equilibrium volume fractions. Here, the rigidity of the gel proves to have a profound effect on the equilibrium thermodynamics: since the polymer network is maintained by permanent chemical cross-links, its connectivity should not change under deformation. Therefore, unless the gel is torn, the polymer network should remain contiguous even while supporting coexistent phases, as shown in figure 11. Spatial inhomogeneities in $\phi$ therefore result in anisotropic stretching of the polymer network. The phase diagram in figure 10 is thus a simplistic picture of swelling thermodynamics, starting from the placement of the critical point. The remainder of this Topical Review is concerned with the consequences of rigidity for the swelling phase behavior.

We begin with the fate of thermodynamic stability and how elastic effects alter the critical point of the gel. Consider again the picture presented by Le Chatelier’s principle, where spatial variation in density at positive compressibility results in a spatial variation in pressure, which causes mass-flow that corrects the spatial variation. Conversely, for negative compressibility, low-density regions have increased pressure compared with high-density regions, leading to runaway mass-flow that causes density variations to grow. In an isotropic polymer gel, this runaway mass-flow that occurs for negative compressibility leads to spatial variation in the polymer volume fraction $\phi$. In order to remain contiguous, the polymer network deforms inhomogeneously, which has a free energy cost that is not accounted for in the classical analysis presented so far. In fact, the gel is stable for small, negative values of the compressibility, as confirmed by light-scattering measurements [51]. It should be noted that materials with negative bulk modulus are generically unstable, since such a material will increase in volume under applied pressure. However, since this type of instability is predicated on the absence of inhomogeneity, it is difficult to observe in bulk materials [52, 53].

In order to concretely formulate this discussion, let us construct the coarse-grained free energy $F$ of an inhomogeneous polymer gel, given a reference state $\mathcal{R}$ with homogeneous polymer volume fraction $\phi_{0}$. This free energy is given by

[TABLE]

where the first term represents the free energy cost of spatial variations in $\phi$ and is related to the correlation length $\xi$ of fluctuations in the volume fraction $\phi$, i.e., $c\sim\xi^{2}$. The free energy density $\tilde{\mathcal{F}}$ is a density with respect to the reference state $\mathcal{R}$ rather than the current state of the gel; the two are related by $\tilde{\mathcal{F}}=(\phi_{0}/\phi)\mathcal{F}$. Fluctuations in the state of the gel take points $\mathbf{r}$ in the reference state to points $\mathbf{R}=\mathbf{r}+\delta\mathbf{r}$ in the current state of the gel. The homogeneous gel is stable if these small fluctuations always increase the total free energy. Therefore, stability is determined by the second variations of the free energy $\delta^{2}F$.

The fluctuations $\delta\mathbf{r}$ are represented by a displacement field $\mathbf{u}(\mathbf{r})\equiv\mathbf{R}(\mathbf{r})-\mathbf{r}$. The deformation matrix can therefore be written as

[TABLE]

from which the polymer volume fraction $\phi=\phi_{0}/({\rm det}\,\Lambda)$ can be approximated as

[TABLE]

so the variation of the volume fraction is given by $\delta\phi=-\phi_{0}\partial_{i}u_{i}$. Therefore, the second variation of the free energy is given by

[TABLE]

where $E_{ijkl}$ is the elasticity tensor, found by expanding $\tilde{\mathcal{F}}$ to second order in $\partial_{i}u_{j}$, and is thus a function of temperature $T$ and volume fraction $\phi_{0}$. As the reference configuration corresponds to a homogeneous, isotropic gel, the elasticity tensor has the form

[TABLE]

where $\mu$ is the shear modulus and $K$ is the bulk modulus. Note that the bulk modulus is given by $K=\phi_{0}^{2}(\partial^{2}\tilde{\mathcal{F}}/\partial\phi^{2})|_{\phi_{0}}$, which can adopt negative values when the mixing part of the free energy is near a local maximum with respect to $\phi_{0}$. If we take the approximation that the gel occupies all of space, then it is useful to use the Fourier representation of the displacement field,

[TABLE]

and the second variation in the free energy becomes:

[TABLE]

We can then decompose the displacement field $\mathbf{u}$ into a longitudinal part $u^{\ell}_{\mathbf{k}}\hat{\mathbf{k}}$ and a transverse part $\mathbf{u}^{t}_{\mathbf{k}}\times\hat{\mathbf{k}}$, which yields two independent contributions to $\delta^{2}F$, namely that due to longitudinal fluctuations,

[TABLE]

and that due to transverse fluctuations,

[TABLE]

There is a correlation length $\xi_{\ell}\equiv[c\phi_{0}^{2}/(K+4\mu/3)]^{1/2}$ associated with longitudinal fluctuations, whereas the transverse fluctuations do not have an associated length scale. Since the longitudinal fluctuations correspond to fluctuations in the polymer volume fraction $\phi$, the longitudinal fluctuation correlation length $\xi_{\ell}$ also describes the polymer volume fraction correlation length. Therefore, we can conclude that whereas the critical point for fluids is at vanishing bulk modulus $K=0$, the critical point for polymer gels is at vanishing longitudinal modulus $K+4\mu/3=0$ [54]. Note that for the classical rubber elasticity used in the Flory-Rehner model, the shear modulus $\mu$ is always positive. However, as we have shown, the osmotic bulk modulus $K=\phi\partial\Pi/\partial\phi$ can be negative. Whereas a negative bulk modulus marks the loss of stability for fluids and homogeneous elastic solids, polymer gels lose thermodynamic stability only when $K<-4\mu/3$. The stabilizing effect of the shear rigidity near the critical point is an effect of inhomogeneity produced by thermal fluctuations. It is remarkable that the critical behavior of polymer gels, amorphous solids primarily composed of liquid, should share some fundamental similarities with exotic solid state alloys that undergo structural phase transitions, such as ferroelastic and martensitic transformations, and are also capable of exhibiting negative elastic moduli [55].

As we have shown, there are features of the swelling behavior of gels that can be understood in analogy with the phase behavior of fluids. However, as we have alluded to, there are key differences due to effects of shear rigidity. While there is loss of stability for longitudinal fluctuations when $K+4\mu/3\leq 0$, it should be remembered that the transverse fluctuations retain stability since $\mu>0$ [56]. The fact that critical fluctuations in the gel are not seen at $K=0$ means that the gel cannot be simply modeled by the Landau theory used for fluids, summarized in equation (1ajan). Instead, as shown by Golubović and Lubensky in [57], there is an additional term of the form $\mu(\int{\rm d}^{3}x\phi)^{2}$ that must be added to the model free energy. Since this term is non-local, meaning that it cannot be folded into a free energy density, it yields a departure from the Landau theory of simple isotropic fluids. Furthermore, it is long-ranged, meaning that the equilibrium value of the field $\phi(\mathbf{x})$ depends on the value of $\phi$ at all other points in the gel. More concretely, we can see the effect of this long-range interaction when one considers equilibrium phase coexistence between swollen and deswollen phases of gels, where the gel shape plays a role in determining the equilibrium values of $\phi$. This is a notable departure from the behavior that we expect from fluids, where the density of the fluid at any point in its bulk is independent of the shape of the container.333The density of a fluid, as well as its fluctuations, do depend on shape when one considers finite samples and points near the boundary of the fluid [58]. However, for fluids, which have only local interactions, these boundary effects decay as points in the bulk are considered, and are therefore distinctly different than the shape-dependence that is seen in gels.

4 Arrested deswelling: an exotic way to phase-separate

Whilst the transition between swollen and deswollen gels is in many ways similar to the transition between gas and liquid phases, shear rigidity alters some key characteristics of the transition. As demonstrated in the previous section, the definition of the swelling critical point for gels is modified due to the stabilizing effect that shear rigidity has on thermal fluctuations. Shear rigidity also has profound consequences for the phase-coexistent equilibria of polymer gels, as the mechanical equilibrium condition requires balance of a generally anisotropic stress tensor. Thus, particularly in the case of macrophase separation, where macroscopic domains of swollen and deswollen gel give rise to such an anisotropic stress distribution, the result is a deformation of the gel at similarly long length scales. As the response of a solid to an applied stress distribution depends on the shape of the solid, we expect that the conditions and configurations of the phase coexistent equilibrium states depend on the shape of the gel.

Experiments on cylindrical samples of an ionic polymer gel have demonstrated phase coexistence at constant ambient osmotic pressure [49]. The mass of swollen compared to deswollen gel in a single sample is controlled by temperature, which is similar to coexistent gas and liquid phases. However, there is indeed a noticeable dependence on shape. If an unconstrained cylinder-shaped gel, initially in the swollen phase, is brought to the coexistence regime, the deswollen phase nucleates at the two ends, as shown in figure 12(a); similarly, a swollen phase will grow from the ends of a deswollen cylinder. However, if such a gel is instead stretched uniaxially, the new phase nucleates from the center, as shown in figure 12(b).

Demonstrating phase coexistence in neutral gels has proven elusive due to the relatively narrow temperature range of equilibrium phase coexistence, spanning less than $0.1^{\circ}{\rm C}$, at constant osmotic pressure $\Pi$ [49]. However, it has been shown that this temperature range can be broadened through the application of mechanical stress. For example, equilibrium phase coexistence has been demonstrated in cylindrical samples of neutral gels that are stretched via a mechanical constraint applied to the ends of the gel [59, 60, 61], as shown in figure 12 (b). This is another realization of the effect of shear rigidity: since the equation of state for the osmotic pressure depends on the deformation matrix $\Lambda$, osmotic isobars are affected by anisotropic deformation.

Phase coexistence has been carefully achieved at constant osmotic pressure. However, phase coexistence at constant volume remains unexplored, likely because it is potentially challenging to fabricate a volume-constraining material grafted onto the boundary of the gel. One alternative approach is to take advantage of the rich deswelling kinetics of polymer gels. Matsuo and Tanaka [9] studied the equilibration of 0.1 to 1 millimeter-radius spherical samples of neutral pNIPAM under heating and cooling through the first-order phase transition at zero osmotic pressure. In their experiments, they observed that spheres, initially swollen at low temperature, that are rapidly heated past the deswelling transition temperature at $\sim\mkern-4.0mu32^{\circ}$C have a two-step deswelling process. Immediately after the rapid heating process, there is some deswelling over the first few seconds, after which deswelling halts for tens of seconds. During this pause in deswelling, referred to as the “plateau period,” there is thin skin of deswollen-phase gel at the boundary of the sample that is effectively impermeable to the solvent. Towards the end of the plateau period, the deswollen skin appears thinner in some regions and thicker in others. These thicker regions form a network-like structure of edges and vertices, with the thinner regions as the faces, reminiscent in morphology to a foam, that spans the surface of the gel. The resulting inhomogeneous stress distribution about the skin causes the thinner-skinned regions to balloon outward, whereas the thicker-skinned regions are creased inward. As a result of the ballooning, the thinner-skinned regions experience a large extensional strain tangential to the surface of the gel, facilitating the passage of solvent out of the gel. This allows the gel to resume deswelling, eventually equilibrating and becoming a spherical deswollen gel. Similar experiments on cylindrical [62, 63, 64, 65] and toroidal [66] samples of pNIPAM gel reveal even richer phenomena. In addition to the eventual formation of similar balloon-shapes on the surface of the gel near the end of the plateau period, which lasts for minutes in the experiments on toroidal gels, there is also a dramatic shape change, where the torus adopts, in some cases, a saddle-like or “PringleTM”-like morphology.

The plateau period after rapid heating is a prolonged time during which the volume of the polymer gel is essentially unchanging. During this time, the interior of the gel is trapped in the swollen phase, even though it is at a temperature where the deswollen phase is an absolute minimum of the free energy. If we ignore the small initial solvent loss in forming the very thin deswollen skin, this interior is effectively under a constant-volume constraint and is in a state that is far from the free energy minimum. Much like a fluid in similar conditions, we expect that the gel will phase-separate, forming a high polymer volume fraction region at the expense of also forming a low polymer volume fraction region. Since the plateau period is a relatively long-lasting part of the deswelling kinetics of the polymer gel, it is reasonable to assume that the phase-separation will approach an equilibrium phase-coexistent state before solvent starts to leak out at an appreciable rate. Equilibrium phase coexistence in this situation is made possible because the thin deswollen skin is effectively impermeable, so the swollen interior is out of chemical equilibrium with the solvent bath, as illustrated in figure 13. Thus, the rapid heating accomplishes (i) a quench across the first-order phase transition, which takes the gel far from global equilibrium, (ii) the introduction of a volume constraint for the swollen interior, and (iii) a prolonged period during which the swollen interior can reach a phase separated state.

4.1 Revisiting the common tangent construction

To determine the phase coexistent equilibria for polymer gels of some given geometry, we need to minimize the total deformation free energy $\Delta F$ under the constraint that the total volume $V$ remains constant. For convenience, take the reference configuration $\mathcal{R}$ to be the swollen gel immediately before phase-separation, where the gel is at a homogeneous polymer volume fraction $\phi_{0}$. We adopt the following general form for the free energy density $\Delta\tilde{\mathcal{F}}$ of the gel in the reference state,

[TABLE]

where $\mu_{0}$ is the shear modulus corresponding to the reference state. Note that the first term is a general form for the elastic contribution, originating from classical theory of rubber elasticity [29], appearing in the Flory-Rehner model, as well as more sophisticated descriptions [8]. The free energy density $\Delta\hat{\mathcal{F}}(T,\phi)$ is the remaining part of the free energy density that only depends on the polymer volume fraction $\phi(\mathbf{r})$ at each point in space. Note that in expressing $\Delta\tilde{\mathcal{F}}$ in terms of the Flory-Rehner theory, $\Delta\hat{\mathcal{F}}$ contains the mixing part of the free energy, as well as the additional terms of the Flory-Wall network entropy describing translational degrees of freedom of the chains; the shear modulus is given by $\mu_{0}=n^{0}_{ch}k_{B}T$, where $n^{0}_{ch}$ is the chain density in the reference state of the gel. The total free energy $\Delta F$, given by

[TABLE]

incorporates a constant-volume constraint, enforced by the Lagrange multiplier $p$. Since we seek a description of macroscopic phase-separation, we do not include the free energy cost of microscopic variation, $c|\bm{\nabla}\phi|^{2}$, as in equation (1ajbe); this amounts to neglecting thermal fluctuations of the state functions, which is justified as long as the gel is sufficiently far from the critical point.

Now consider phase-separation that forms a solvent-poor region with polymer volume fraction $\phi_{p}$ at the cost of forming a solvent-rich region with polymer volume fraction $\phi_{r}$. Starting from a swollen phase of $0<\phi_{0}\ll 1$, we expect that a quench deep into the deswollen phase will result in \phi_{p}^{-1}\ll\phi_{0}^{-1}{\;\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}\;}\phi_{r}^{-1}. Furthermore, the solvent-poor and solvent-rich regions are distinct phases assuming that the gel is far enough from the critical point. Therefore, we expect that $\phi_{p}$ and $\phi_{r}$ deviate very little from their representative values in their respective regions. Ignoring these deviations, the volume conservation condition is given by

[TABLE]

where $f$ is the fraction of gel in the reference configuration that will be solvent-poor. We can solve for the fraction $f$, yielding

[TABLE]

which is identical to the Lever Rule of phase coexistence in fluids, equation (1ah).

Now consider a spherical gel, such as that in the experiment of Matsuo and Tanaka [9]. We use spherical coordinates $(r,\theta,\phi)$ in the reference configuration and $(R,\Theta,\Phi)$ in the target configuration, such that reference configuration points are $\mathbf{r}=r\hat{\mathbf{r}}$ and target configuration points are $\mathbf{R}=R\hat{\mathbf{R}}$. Due to the symmetry of the sphere, we expect that the phase-coexistent equilibrium maintains the spherical symmetry; broken symmetry may arise from instability of this equilibrium. Therefore, points $\mathbf{R}$ in the phase coexistent state depend only on the radial coordinate $r$ in the reference state; the deformed sphere is expressed entirely in terms of $R(r)$, where $R$ is the radial coordinate of the phase coexistent configuration. Consequently, deformations maintain conformal symmetry so that we can take $\Theta=\theta$ and $\Phi=\phi$ and $\hat{\mathbf{R}}=\hat{\mathbf{r}}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. The transpose of the deformation matrix $\Lambda$, as expanded in spherical coordinates, is given by

[TABLE]

with components,

[TABLE]

In the reference geometry, let $a$ be the outer radius of the sphere and let $b$ represent the radius of the interface between solvent-rich and solvent-poor regions, as in figure 14(a). Given that the sphere starts swollen and that in the experiments, there is a thin deswollen-phase skin on the boundary of the sphere, we will assume that the solvent-poor region grows from the deswollen-phase skin, inward. Therefore we will take the region $r<b$ to be solvent-rich and the region $b<r<a$ to be solvent-poor. In general, $R(r)$ is a continuous function since the gel must remain connected throughout the deformation. However, it is expected that there is a change in the behavior of ${\rm d}R/{\rm d}r$ at the interface $r=b$ because whereas the solvent-rich layer should increase in radius, implying that $R(b)>b$, the solvent-poor layer should become thinner after deformation, so $a-b>R(a)-R(b)$, or $R(b)-b>R(a)-a$. Therefore, we expect that ${\rm d}R/{\rm d}r$ is typically greater than 1 for $r<b$ and is typically less than 1 for $b\leq r<a$. If the solvent-poor region is taken to be much smaller than the solvent-rich region, i.e., $f\ll 1$, then $(a-b)/b\ll 1$, and $R(r)$ adopts a piecewise-linear functional form,

[TABLE]

where $\Lambda_{1}>1$, reflecting the further swelling of the swollen interior, and $0<\Lambda_{2}<1$, reflecting the deswelling of the shell region. The deformation matrix in the solvent rich region is given by $\Lambda_{ij}\approx\Lambda_{1}\delta_{ij}$; since ${\rm det}\,\Lambda=\phi_{0}/\phi$, we find that $\Lambda_{1}=(\phi_{0}/\phi_{r})^{1/3}$. In the solvent-poor region, we find that the deformation matrix is both inhomogeneous and anisotropic, with $\Lambda_{Rr}\approx\Lambda_{2}$ and $\Lambda_{\Theta\theta}=\Lambda_{\Phi\phi}\approx\Lambda_{1}$, with $r$-dependence appearing at higher order in $(a-b)/b$. Therefore, we have that $\Lambda_{1}^{2}\Lambda_{2}=\phi_{0}/\phi_{p}$, giving the result $\Lambda_{2}=(\phi_{0}\phi_{r}^{2})^{1/3}/\phi_{p}$. The total deformation free energy for the sphere $\Delta F_{\rm sphere}$, consisting of contributions from the swollen core and deswollen shell, is given by

[TABLE]

Minimizing $F_{\rm sphere}$ with respect to $\phi_{r}$, $\phi_{p}$, $p$, and $f$ yields three equilibrium equations, namely

[TABLE]

as well as the lever rule. However, the result cannot be cast as a common tangent construction, as evidenced by the dependence of $p$ on $f$ in equation (1ajbvb). Ultimately, this results from the anisotropic stress of the solvent-poor region, which comes about through the coherency strain: as the polymer network remains contiguous, the components of the deformation matrix that are tangent to the phase interface, namely $\Lambda_{\Theta\theta}$ and $\Lambda_{\Phi\phi}$, must be continuous through the interface. Therefore, the deformation matrix for the solvent-poor region shares two of its three components with the solvent-rich region; the third component, $\Lambda_{Rr}$, which is normal to the phase interface, is discontinuous and therefore takes on two independent values for the two regions. Note that this is alleviated for effectively one-dimensional gels: a gel that is constrained to undergo uniaxial deformation normal to the phase interface does not experience coherency strain at the interface and therefore the phase-coexistent equilibrium obeys the common tangent construction [67].

For a concrete solution, we use the Flory-Rehner model to provide a concrete form for $\hat{\mathcal{F}}$ and minimize the free energy in equation (1ajbu) numerically. As shown in figure 14(b), coexistence between separate swollen and deswollen phases happens above certain transition values of $\chi_{1}$, depending on values of $\phi_{0}$. For values of $\chi_{1}$ greater than the transition points, the equilibrium volume fraction bifurcates into a high-$\phi$ branch, corresponding to solvent-poor gel of volume fraction $\phi_{p}$, as well as a low-$\phi$ branch, corresponding to solvent-rich gel of volume fraction $\phi_{r}$. Whereas there is a discontinuous jump from the initial volume fraction $\phi_{0}$ to solvent-poor gel $\phi_{p}$, the solvent-rich volume faction decreases continuously from $\phi_{0}$ above the transition. The reason for this is attributed to the manner at which phase-separation occurs, namely a form of heterogeneous nucleation. Since the solvent-poor gel grows from the boundary of the gel, inward, the fraction $f$ of solvent-poor gel grows continuously, as shown in figure 14(c). As the growth of solvent-poor gel comes at the cost of diluting the solvent-rich region, the continuous growth of the solvent-poor shell results in continuous solvent addition to the core-region; therefore $\phi_{r}$ decreases continuously as $f$ and $\phi_{p}$ grow. This continuous growth of the solvent-poor region is consistent with predictions seen elsewhere, in the context of solvent-poor phase growing around a hole [68]. Note that the transition value of $\chi_{1}$ shifts to higher values as $\phi_{0}$ is set to higher values. This is because for hight values of $\phi_{0}$, there is less interaction between solvent and polymer and therefore the total energetic cost of maintaining the homogeneous phase is less; this energetic cost increases for lower values of $\phi_{0}$.

5 Large shape change via thermodynamics

The characteristic swelling behavior of polymer gels makes them practical for many applications. Their ability to absorb and retain solvent, enables them to remove unwanted liquids, such as water, making them attractive for cleaning applications [69]. At the same time, the ability to expel a liquid given certain stimuli has led to drug-delivery applications [70, 71, 72]. Furthermore, as we have outlined, the equilibrium swelling thermodynamics of polymer gels is well modeled by the Flory-Rehner model, given certain fitting parameters. However, as we have discussed in the previous section, there are interesting consequences that emerge when a polymer gel is brought to a state in which having a single, homogeneous polymer volume fraction $\phi$ corresponds to a thermodynamically unstable situation, leading to a phase-coexistent equilibrium. Moreover, due to the mechanism of skin formation after a rapid quench from the swollen phase to the deswollen phase, this situation is attainable if the volume phase transition is approached rapidly. Thus, while the phase-coexistent state remains relatively unexplored, understanding it is, potentially, of considerable practical importance as it dramatically alters the expected behavior of a gel. While in many cases, the deswelling-arrested, phase-coexistent state of a polymer gel may be an unwanted, troublesome feature of the polymer gel thermodynamics that must be avoided, it is possible that phase coexistence and the accompanying shape change can be useful, and seen as a potential design feature. This is not without precedent, as exemplified by extreme mechanics, in which mechanical instability, traditionally regarded as a nuisance to be avoided, is harnessed in the design of novel material properties and responses to achieve shape change that would otherwise be unattainable.

5.1 Preliminaries: Extreme mechanics

Unlike many other rigid materials, the immutable network topology of polymer gels and other elastomeric materials enable them to undergo large elastic deformations. As swollen polymer gels are primarily composed of solvent by weight, they are very soft materials that are comparable to soft biological systems. Indeed, changes in the swelling state of a polymer gel can be used to mimic the growth of soft tissues. This combination of elasticity and swelling enables experiments on growth-induced instabilities as seen in nature. In particular, swelling of polymer gels or other materials that are tethered to a rigid substrate exhibit a wide range of surface patterns, ranging from wrinkles, resembling fingerprints, to folds and creases, much like the sulci of brain tissue [17, 19, 18, 13, 12, 73, 74]. Interestingly, the appearance of these patterns are seemingly chosen at random, given the underlying symmetry of the substrate. Such spontaneous symmetry breaking of gels and other elastic materials under some sort of mechanical constraint is typically the result of mechanical instability [75, 76]. Furthermore, the result of such pattern formation often yields new, and perhaps a priori unexpected mechanical response [77, 78, 79, 80, 25].

Perhaps the most famous example of mechanical instability, due to its analytical tractability and appearance in structural engineering and nature, is Euler buckling. For concreteness, we will consider a three-dimensional material that is described by an isotropic linear elastic energy

[TABLE]

where $E$ is the Young’s modulus, $\nu$ is the Poisson ratio, and $u_{ij}\equiv(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}\mathbf{u}\cdot\partial_{j}\mathbf{u})/2$ is the finite symmetric strain tensor, corresponding to the displacement field $\mathbf{u}(\mathbf{r})$. The nonlinear terms in the strain tensor, which are usually neglected, are retained here in order to properly compute the second variation of the elastic energy when determining stability. If this material is formed into a slender column of length $L$ and circular cross-section radius $a\ll L$ that is under a compressive stress $\sigma_{zz}=T/(\pi a^{2})$, as illustrated in figure 15(a), then the total energy is given by

[TABLE]

where $\Delta L$ is the a change in length. Assuming that the rod is symmetric about the $z$-axis, the deformation will also be axisymmetric so that the strain will have longitudinal $u_{zz}$ and transverse $u_{\perp\perp}$ parts that, in the slender rod approximation, can be approximated as constants, i.e., $u_{zz}=\Delta L/L$ and $u_{\perp\perp}=\Delta a/a$. The total energy can therefore be written as

[TABLE]

The new, stressed equilibrium of the rod is found by minimizing the total energy with respect to $u_{\perp\perp}$ and $u_{zz}$, yielding equilibrium strains

[TABLE]

which shows that for compression $(T>0)$, the rod decreases its length and for most elastic materials $(\nu>0)$, it increases its width proportionally. This shortening of the length of the rod details the expected change in mechanical equilibrium of the rod from an unstressed configuration to a stressed configuration. However, it is not guaranteed that the stressed equilibrium is stable to deformations of the rod that break its symmetry 444We will assume, however, that the unstressed equilibrium is stable. For polymer gels, as we have shown, this is a nontrivial assumption that depends on temperature and composition.. Indeed, it has been long known that if one overloads a column, the column may bend, as illustrated in figure 15(b). To probe the stability of the rod with respect to deflections that bend the rod, we require a description of the elasticity of bending deformations. It can be shown (see, e.g., [81]) that if the centerline of the rod is deflected from a straight configuration to a curved configuration, where it acquires a non-zero curvature $\kappa(z)$, the associated elastic energy cost is

[TABLE]

where $B$ is the bending modulus of the rod, a composite of the Young’s modulus $E$ and the second moment of area of the cross-section of the rod; for rods of circular cross-section, $B=(\pi/4)Ea^{4}$. Stability of the straight configuration against transverse deflections $\mathbf{u}_{\perp}$ of the rod is ensured when the second variation of the total energy, $\delta^{2}E$, with respect to these deflections is positive. Expanding the energy $E$ about the stressed equilibrium, that is, where $\Delta L\approx-LT/(E\pi a^{2})+\delta(\Delta L)+\delta^{2}(\Delta L)$, the second variation of the energy is given by

[TABLE]

which can then be expressed in terms of the transverse deflections $\mathbf{u}_{\perp}$. This is done by noting that the shape of the rod is given by $\bm{\gamma}(z)=z\hat{\mathbf{z}}+\mathbf{u}_{\perp}(z)$ so that the curvature change $\delta\kappa$ after deformation is given by

[TABLE]

and the second length variation $\delta^{2}L$ after deformation is given by

[TABLE]

Therefore, the second variation of the energy is

[TABLE]

In general, for axisymmetric rods in three-dimensions, $\mathbf{u}_{\perp}$ can point anywhere transverse to the axis of symmetry (here, the $z$-axis) and the buckled, bent configurations spontaneously break a continuous symmetry, namely, rotations about the $z$-axis. For simplicity, we will restrict our attention to symmetry-breaking deflections that keep the rod configuration confined to a plane, namely the $xz$-plane, so that $\mathbf{u}_{\perp}=u_{\perp}\mathbf{\hat{x}}$. Now, if we assume that the rod is fixed at its two ends such that $u_{\perp}(0)=u_{\perp}(L)=0$ then we can expand the deflection field as a sum of sine functions, namely

[TABLE]

and the second variation of the energy is given by

[TABLE]

Note that the form of the bending energy used assumes a slowly varying curvature along the length of the rod, so that the second variation of the energy, as shown above, is only valid when the wavenumber of the bending mode, $(n\pi)/L$ is much less than $1/a$. Therefore, for low compressive force $T$, the energy remains positive for small transverse deflections. However, for $T>T^{(n)}_{c}$, where the critical compressive force is $T^{(n)}_{c}=Bn^{2}\pi^{2}/L^{2}$, the total energy may decrease by bending, forming $n/2$ wavelengths. Hence, for $T>T_{c}\equiv T^{(1)}_{c}$, the rod is susceptible to buckling, with the critical compressive force $T_{c}$ marking the onset of the instability.

5.2 Extreme thermodynamics

There is an analogy that can be drawn between buckled rods and thermodynamic phases [82]. In particular, the onset of mechanical instability resembles the onset of thermodynamic instability at a critical point. Indeed, we can even take the Landau phenomenological approach and model the equilibrium of the bent rod, restricted to the lowest wavenumber mode, here $n=1$, via a model elastic energy of the form

[TABLE]

where $s>0$ is a parameter that stabilizes the amplitude of the deflection and $h$ is an external force that couples to the deflection. This simple model recovers the spontaneous symmetry breaking due to Euler buckling at $h=0$ [see figure 15(c)] but also is able to model the action due to an external force $h$ that is applied to the center of the rod (at $z=L/2$) oriented in the $\mathbf{\hat{x}}$ direction. Interestingly, whereas $u_{\perp}\propto h$ at $T<T_{c}$, the response is much more complicated at $T>T_{c}$, where there is a discontinuous snap-through from one curvature to the opposite, as depicted in figure 15(c). This snap-through appears as the mechanical analogy of a first-order phase transition. Interestingly, whereas the rod supports transverse elastic waves in both its straight and buckled configurations, the vibrational frequency of such waves goes to zero as the buckling threshold is reached, due to loss of elastic stability. This is an elastic analogy of the “critical slowing down” of the relaxation time of perturbations to thermodynamic systems near a critical point [83, 84].

Euler buckling and snap-throughs of a slender rod are examples of how bifurcations and limit points in the equilibrium phase diagram of an elastic body exhibit parallels with phase transitions [82]. Indeed, as pointed out in a recent review by Douglas Holmes [25], the analogy is useful for understanding mechanical instability for a wide variety of systems. Additionally, Landau theory has proven to be a useful tool for translating mechanical instability into the language of phase transitions, providing additional insight, for example, into spin arrangements in magnetic dots [85]. In the case of 2D elastic sheets that undergo a wrinkling instability, a similar transfer of ideas has revealed the smectic-like behavior of the wrinkle patterns [79].

It is thus natural to inquire as to how far these analogies between mechanics and thermodynamics may be taken. In particular, if mechanical instability can be used to design new material response, can we then harness thermodynamic instability in a similar way? Can a case be made for extreme thermodynamics, in which materials, such as polymer gels, are tuned near a point of thermodynamic instability in a manner that yields new, interesting behavior?

Thermodynamic instability, like mechanical instability, signals the development of multiple free energy minima, as opposed to a single, as well as the possibility of spontaneous symmetry breaking. A hallmark of the appearance of multiple metastable equilibrium states is the ability to support coexistent phases within the same sample. As we have discussed, phase coexistence can be achieved in polymer gels via rapid heating from the swollen phase to the deswollen phase. In experiments on spheres, such phase coexistence leads to the formation of surface crease patterns; experiments on tori reveal additional buckling behavior of the toroidal shape, as shown in figure 16(e-i). The source of these patterns is an internal stress distribution due to separation into coexistent swollen and deswollen phases. Therefore, one might regard the appearance of the patterns as mechanical buckling due to internal stress generated by the allocation of solvent in the gel [66].

However, there is a feedback: the arrangement of the coexistent phases, or the distribution of the solvent within the gel, depends on the mechanical stress distribution. This is illustrated by the example of the phase-coexistent equilibria of a sphere in section 4.1, in which anisotropic deformations of the polymer network due to coherency strain between the two interfaces altered the equilibrium polymer volume fractions of each of the two phases. As we shall demonstrate in the next section, this coherency strain has a real effect on the distribution of solvent in polymer gel tori, leading to internal stresses that cause buckling. We can therefore conclude that such buckling is the result of phase coexistence, and thus the presence of thermodynamic instability, showing that thermodynamic instability can be harnessed to produce new behavior (buckling and pattern formation) that is distinct from the normal behavior in the thermodynamically stable regime (swelling and deswelling).

5.3 The case of toroidal polymer gels

Tori are characterized by two length scales; it is convenient to use the tube radius $a$ and the ring radius $R$, as depicted in figure 17. In order to address phase coexistence in a toroidal sample of polymer gel as well as the buckling instability in a way that is analytically tractable, we will take the slender rod approximation $a\ll L$ that was used to study Euler buckling. This is equivalent to working in the limit where $\kappa a\ll 1$, where $\kappa=1/R$ is the initial, fabricated curvature of the toroidal centerline. We can therefore approximate the toroid as a slender cylinder of radius $a$ and length $L=2\pi R$, whose ends are identified, ensuring that it still has the topology of a solid torus. Curvature is then incorporated as a perturbative correction of size $\kappa a$ to this flat limit of the torus.

Under similar conditions to the sphere, phase-coexistent equilibria of the flat limit of the torus can be treated in a similar manner. We utilize the same approximations as those used in section 4.1. Proceeding in cylindrical coordinates, the reference configuration is parametrized by $(\rho,\theta,z)$ and the phase-separated target configuration by $(P,\Theta,Z)$. In the limit in which the deswollen shell is thin, we can approximate the radial coordinate $P$ as a piecewise, linear function of $\rho$, such that $P(\rho)\approx\Lambda_{1}\rho$ for $\rho<b$ and $P(\rho)\approx\Lambda_{1}b+\Lambda_{2}(\rho-b)$ for $b\leq\rho<a$, where $\Lambda_{1}$ and $\Lambda_{2}$ represent re-scaling of points transverse to the $z$-axis. Next, translational symmetry along the $z$-axis ensures that the deformation matrix is not a function of $z$. Since points $z=0$ and $z=L$ are identified under the periodic boundary conditions of the torus, translational symmetry requires that $Z=\Lambda_{\ell}z$, where $\Lambda_{\ell}$ describes the longitudinal stretch or compression of the torus. Therefore, in the solvent-rich core, the polymer volume fraction $\phi_{r}=\phi_{0}/(\Lambda_{1}^{2}\Lambda_{\ell})$ and in the solvent-poor shell, $\phi_{p}=\phi_{0}/(\Lambda_{1}\Lambda_{2}\Lambda_{\ell})$. Accordingly, the total free energy $F^{0}_{\rm torus}$ of the phase-separated, unbent torus is given by

[TABLE]

Values of $\phi_{p}$, $\phi_{r}$, and $f$ that minimize this free energy are plotted in figure 18, assuming the Flory-Rehner model. Comparing the results with figure 14(b,c), notice there are some marked differences between the variation in $\phi_{r}$ and $f$ with $\chi_{1}$ between a sphere and a torus. For example, whereas $\phi_{r}$ for the sphere steadily decreases with $\chi_{1}$, there is not much change in $\phi_{r}$ for the flat torus until larger values of $\chi_{1}$. This is due to the additional length-change degree of freedom of the flat torus: for sufficiently small values of $\chi_{1}$, the deformation is mainly of the cross-section, so that $\Lambda_{\ell}\approx 1$, and involves a large shrinking of the thin shell with a modest stretch of the core, resulting in $\phi_{p}\sim\Lambda_{2}^{-1}$ and $\phi_{r}$ roughly constant; for larger values of $\chi_{1}$, deformations of the length become more important as the shell radius change becomes costly, so that $\phi_{p}\sim\phi_{r}\sim\Lambda_{\ell}^{-1}$, shown in the inset of figure 18(b). There is also a drastic difference in $f(\chi_{1})$ for the two geometries, where the fraction of the solvent-poor phase appears to level-off with increasing $\chi_{1}$ for the sphere but rapidly grows for the toroid. This is also due to the additional degree of freedom present in the toroid and there is eventually a similar level-off of $f$ for higher values of $\chi_{1}$ (not shown), due to the activation of the length change.

Next, we consider the effects of curvature on the phase-coexistent equilibria. The curvature of a toroidal centerline lifts the rotational symmetry of the flat, cylinder-like torus. Consequently, we must revise the assumption that the interface between the solvent-rich core and the solvent-poor shell is axisymmetric. Furthermore, this also breaks the assumption of axisymmetric deformation, requiring a more general representation of the deformation matrix $\Lambda$. In order to address this asymmetry, we move from the cylindrical coordinate system to a more general one in which the straight $z$-axis is replaced with a closed parametric curve $\bm{\gamma}(s)$, where $s\in[0,L)$ is the arclength parameter of the centerline, defined such that $|\partial_{s}\bm{\gamma}|=1$. As long as the centerline is everywhere curved, we can uniquely construct the Frenet-Serret frame $\{\mathbf{\hat{t}},\mathbf{\hat{n}},\mathbf{\hat{b}}\}$, where $\mathbf{\hat{t}}(s)\equiv\partial_{s}\bm{\gamma}$ is the unit tangent vector, $\mathbf{\hat{n}}(s)\equiv(\partial_{s}\mathbf{\hat{t}})/|\partial_{s}\mathbf{\hat{t}}|$ is the unit normal vector, and $\mathbf{\hat{b}}(s)\equiv\mathbf{\hat{t}}(s)\times\mathbf{\hat{n}}(s)$ is the unit binormal vector, which completes an orthonormal triad at all points on the centerline; this frame is illustrated in figure 19. The rotation rate of this frame along the centerline depends on the curvature $\kappa$ and the torsion $\tau$ of the centerline, via

[TABLE]

a set of geometric relations known as the Frenet-Serret equations [86]. For the fabricated, planar torus of constant curvature $\kappa$, there is no torsion, $\tau\equiv 0$. However, to study buckling-type deformations of the torus, where the ring can adopt non-planar deformations, we require this full geometry of curves where $\tau\neq 0$. Furthermore, the cross-section of the gel can, in general, twist independently of the centerline. We therefore define a material frame $\{\mathbf{\hat{d}_{1}},\mathbf{\hat{d}_{2}},\mathbf{\hat{d}_{3}}\}$ that is adapted to the centerline and can therefore be expressed in terms of the Frenet-Serret frame. It is convenient to choose

[TABLE]

where $\{\mathbf{\hat{d}_{1}},\mathbf{\hat{d}_{2}}\}$ define the transverse frame to the centerline and $\varphi(s)$ is the angle of rotation from the Frenet-Serret frame [see figure 19]. There is a more general set of equations describing the rotation rate of the material frame, namely

[TABLE]

where $\mu\in\{1,2,3\}$ and $\bm{\omega}(s)$ is known as the Darboux vector and is given by

[TABLE]

Therefore, we can represent points $\mathbf{r}$ near an arbitrarily curved centerline as

[TABLE]

where $x_{1},x_{2}$ give coordinates in the cross-section.

We will therefore represent the reference configuration $\mathcal{R}$ of the torus with centerline $\bm{\gamma}$, frame $\{\mathbf{\hat{d}_{1}},\mathbf{\hat{d}_{2}},\mathbf{\hat{d}_{3}}\}$, and coordinates $(s,x_{1},x_{2})$. The cross-section coordinates $(x_{1},x_{2})$ can be written in polar form $(\rho,\theta)$ via $x_{1}=\rho\cos\theta$ and $x_{2}=\rho\sin\theta$. While $\rho=a$ is the boundary of the torus, the interface between the solvent-rich and solvent-poor regions is generally not axisymmetric but adopts a more general form $\rho=b(\theta)$ in these polar coordinates. There are many shapes that the interface can adopt and there is evidence, namely the ballooning patterns that forms on the surface of toroidal gels [see figure 16(f-j)], that the deformation may have a relatively low wavelength along both the toroidal meridian (the $\bm{\hat{\theta}}$ direction) and the toroidal ring (the $\mathbf{\hat{d}_{3}}$ direction). However, in examining the effect of curvature on the phase-coexistent equilibrium, we will focus on longest-wavelength alterations to the coexistence pattern that may arise due to broken axial symmetry. Therefore, let the interface shape $b(\theta)$ adopt the simple form

[TABLE]

which is a single-wavelength deformation of the interface in the cross-section. As long as the amplitude $|\mathbf{p}|$ of this deformation is small, its overall effect is a simple translation of the circular interface such that it is no longer centered on the toroidal centerline as shown in figure 20(a). In effect, it describes a polarization of the distribution of solvent mass within the cross-section of the toroid; we shall refer to $\mathbf{p}$ as the polarization vector. It is possible to describe other moments of the solvent mass distribution, such as an inertia tensor, which corresponds to a two-wavelength deformation of the interface shape in the cross-section. However, the coupling between curvature and polarization is the simplest, as dictated by symmetry; the curvature coupling to other moments occurs at higher powers of the free energy $F$ and thus may be neglected in the low-curvature regime.

In order to represent the equilibrium target configuration $\mathcal{T}$, let us return to the idea of incorporating the effect of curvature as a small correction to the equilibrium phase coexistence in the flat torus. To do this, we will add a step into the deformation process $\mathcal{R}\rightarrow\mathcal{T}$, and call this intermediate step $\mathcal{I}$. In the intermediate configuration $\mathcal{I}$, the toroid consists of points $\mathbf{r^{\prime}}$, given by

[TABLE]

where $\bm{\gamma}^{\prime}$ is the deflected centerline, e.g., due to buckling, with new frame $\{\mathbf{\hat{d}^{\prime}_{1}},\mathbf{\hat{d}^{\prime}_{2}}\}$ [see figure 19]. Furthermore, assume that in going from $\mathcal{R}\rightarrow\mathcal{I}$, the toroid undergoes phase-separation assuming the same deformation matrices as in the axisymmetric, flat torus limit. Therefore, in the solvent-rich region,

[TABLE]

where $x^{\prime}_{1}=\rho^{\prime}\cos\theta^{\prime}$ and $x^{\prime}_{2}=\rho^{\prime}\sin\theta^{\prime}$ yields the polar representation $(\rho^{\prime},\theta^{\prime})$ of the cross-section coordinates in $\mathcal{I}$. The anisotropic deformation of the solvent-poor shell is complicated by a more general, thickness-varying shell; if the solvent-rich core maintains the above isotropic deformation then the continuity requirement across the phase interface results in a deformation $\Lambda_{1}$ tangential to the interface and $\Lambda_{2}$ normal to the interface in the cross-section. In the reference configuration $\mathcal{R}$, the interface tangent is given by $\mathbf{\hat{T}}\equiv\partial_{\theta}(b\bm{\hat{\rho}})/|\partial_{\theta}(b\bm{\hat{\rho}})|$ and normal by $\mathbf{\hat{N}}\equiv\mathbf{\hat{T}}\times\mathbf{\hat{d}_{3}}$; to leading order in $p=|\mathbf{p}|$, defined in equation (1ajbvbzckco), these are given by

[TABLE]

The interface tangent and normal, $\mathbf{\hat{T}^{\prime}}$ and $\mathbf{\hat{N}^{\prime}}$, in the intermediate configuration $\mathcal{I}$ have the same form due to the isotropic deformation of the core, replacing $\bm{\hat{\rho}}\mapsto\bm{\hat{\rho}^{\prime}}$ and $\bm{\hat{\theta}}\mapsto\bm{\hat{\theta}^{\prime}}$. Therefore, in the solvent-poor shell,

[TABLE]

as depicted in figure 20(b). Let $\Lambda^{\prime}$ be the deformation matrix field describing position-dependent changes in length that occur in going from $\mathcal{R}\rightarrow\mathcal{I}$. Then

[TABLE]

where $\Lambda^{0}$ is the deformation matrix corresponding to axisymmetric deformation, namely

[TABLE]

where the top holds for points in the solvent-rich core $\rho<b(\theta)$ and the bottom holds for points in the solvent-poor shell $b(\theta)\leq\rho<a$. The other part $(\mathbbm{1}+u^{\rm ext})$ of the deformation matrix $\Lambda^{\prime}$ describes small changes in length due to deformations of the centerline and the material frame. These small changes are encoded in an “external strain” $u^{\rm ext}$, which is given by

[TABLE]

where $\Delta\omega_{\alpha}\equiv\omega^{\prime}_{\alpha}-\omega_{\alpha}$ is the change in the Darboux vector, i.e., change in the curvature, torsion, and twist of the framed curve degrees of freedom of the gel. Here, we have taken the approximation that the arclength $s^{\prime}$ after phase-separation is proportional to the arclength $s$ before phase-separation and that this proportionality is the longitudinal deformation matrix $\Lambda_{\ell}$, describing the change in length of the gel. Even though it can be expected that curvature and torsion of the toroidal centerline results in a renormalization of $\Lambda_{\ell}$, this effect should be small since changes in the length of the rod are described by deformations of the gel that respect axial symmetry; as curvature breaks axial symmetry, we can expect such renormalization to be a higher order effect than what we seek to describe. The external strain also describes the non-axisymmetric changes in length that occur due to axisymmetric phase-separation when curvature is present; this effect disappears when $\Lambda_{1}\rightarrow 1$, i.e., when there is no transverse deformation of the gel in the solvent-rich region.

The intermediate configuration $\mathcal{I}$ is a poor representation of the phase-separated gel as it assumes the axisymmetric solution in a geometry without axial symmetry. To rectify this, points $\mathbf{R}$ in the actual equilibrium state, that is, in the target configuration $\mathcal{T}$ are obtained from points $\mathbf{r^{\prime}}$ in $\mathcal{I}$ via an “internal” displacement field $\mathbf{u}$ such that

[TABLE]

This displacement field $\mathbf{u}$ represents a “correction” to the shape of the equilibrium gel, resulting in a small, non-axisymmetric strain that is added to the large axisymmetric deformation. By including the intermediate configuration $\mathcal{I}$, we are able to factor the deformation matrix $\Lambda$ into an axisymmetric part $\Lambda^{0}$ and a part that describes non-axisymmetric deformations, namely

[TABLE]

where we have retained only the leading order non-axisymmetric terms, $\partial_{i}u_{j}$ and $u^{\rm ext}_{il}$. Whereas the external strain encodes the shape of the centerline of the toroid, the remaining strain $\partial_{l}u_{i}$ describes degrees of freedom that we can regard as “internal” to the gel and that are allowed to equilibrate given a certain centerline shape. Thus, we call $\partial_{i}u_{j}\equiv u^{\rm int}_{ij}$ the “internal strain.” Finally, it is useful to fold the external and internal strain into a total strain $u\equiv u^{\rm ext}+u^{\rm int}$.

The total free energy $F$ of the gel is then given by

[TABLE]

where $F_{p}$ and $F_{r}$ are the free-energies of the solvent-poor and solvent-rich region, respectively, given by

[TABLE]

with $\int_{p}$ and $\int_{r}$ representing integration over solvent-poor and solvent-rich regions of the gel in the reference configuration $\mathcal{R}$. The third term enforces the volume constraint, with $p$ the Lagrange multiplier. Thus, we can compute the total free energy $F$, given the deformation matrix expressed in equation (1ajbvbzckcrcx) and using the relation $\phi=\phi_{0}/({\rm det}\,\Lambda)$. However, the amount of solvent in each region of the gel, described by the polymer volume fraction $\phi(\mathbf{r})$, should be largely unaffected by the presence of a small, non-axisymmetric strain $u$; only the shape of the regions is affected by such strain. This is a common assumption adopted in theories of elasticity: the “elastic” degrees of freedom are independent of any order parameter describing the phase of the material. For example, in the elastic theory of nematic-phase liquid crystals, gradients in the director field do not greatly affect how “nematic” the material is—such elastic terms are considered transverse to order parameter (see e.g., [32]). This assumption holds as long as the thermodynamic phase of the material is well-defined, failing near the critical point. Similarly, even though the gel in study is in the coexistence region, it consists of two well-defined phases, assuming that it is far enough from the critical point. Therefore, the polymer volume fraction $\phi$, which describes the phase of the gel, has two discrete values, $\phi_{p}$ and $\phi_{r}$. Referring to the decomposition of the deformation matrix (1ajbvbzckcrcx), we therefore require that the strain $u$ satisfies the incompressibility constraint, ${\rm det}(\mathbbm{1}+u)=1$, everywhere within the gel. If this holds, then the volume constraint is maintained independently of $u$. Importantly, while the non-axisymmetric strain $u$ does not affect the volume fraction $\phi$ of the two phases, it does affect the spatial distribution of the two coexistent phases, which is described by the interface shape $b(\theta)$. Any interface shape $b(\theta)$ that does not preserve the continuous rotational symmetry of the gel about its centerline axis results in a non-axisymmetric coherency strain, which is encoded in $u$. We seek to determine the change $\Delta F=F-F^{0}_{\rm torus}$ in free energy due to changing the arrangement of phases in the presence of curvature.

To proceed, we furthermore require that only the symmetric part $\epsilon\equiv(u+u^{T})/2$ of the strain $u$ appears in the free energy in order to properly describe the cost of elastic deformations. Furthermore, in order to enforce the incompressibility constraint ${\rm det}(\mathbbm{1}+u)=1$, the determinant can be expanded using the relation ${\rm det}\,A={\rm exp}({\rm ln}\,A)$ for any matrix $A$, resulting in the condition

[TABLE]

As a result, the free energy change $\Delta F$ due to lifting the axial symmetry is given by

[TABLE]

where $c^{p}$ and $c^{r}$ are anisotropic elasticity tensors that characterize the linear elasticity of the phase-separated “flat” toroid (see [87] for details). Integrals over the solvent-poor and solvent-rich regions are given by

[TABLE]

where the portion of the Jacobian $|\partial\mathbf{r}/\partial(s,\rho,\theta)|$ that depends on curvature contributes a higher-order correction to the free energy; thus these integrals are over cylindrical regions of the gel. Next, the internal displacement field $\mathbf{u}$ is found by finding conditions under which the free energy is an extremum, $\delta\Delta F=0$, at fixed centerline geometry and polarization $\mathbf{p}$, i.e.,

[TABLE]

which yields, to leading order, a displacement field that is linear in $\Delta\omega$, $\omega$, and $b(\theta)$. The result is an effective curve elastic energy

[TABLE]

where $B$ and $C$ are effective bending and twisting moduli, $k>0$ is a constant that couples toroid curvature to interface shape, and $r>0$ characterizes the free energy change due to polarization of the solvent distribution (see [87] for values of these parameters in terms of $\mu_{0}$, $\Lambda_{1}$, $\Lambda_{2}$, and $\Lambda_{\ell}$). The positive value of $r$ implies that there is always a free energy decrease that can be achieved for a polarized solvent distribution. While this is perhaps unexpected, note that a polarized solvent distribution results in the formation of a solvent-poor shell with a thinner region and a thicker region, resulting in a non-uniform stress transmitted across the phase interface. This stress results in a non-axisymmetric deformation of the solvent-rich core that, due to coherency strain, also stretches parts of the shell, leading to a local free energy density increase, whilst compressing other parts of the shell, leading to a local free energy decrease. However, the part of the shell that is stretched is the thinner part, whereas the part that is compressed is thicker, leading to a net free energy decrease, as compared to the free energy cost of interfacial strain with a shell of uniform thickness. As long as the toroid maintains this solvent-poor shell, however, the magnitude $p$ of the polarization is ultimately limited: for large $p$, the thin part of the shell is greatly stretched, resulting in a large elastic penalty. Still, the polarization, and other higher moments of the solvent distribution, represent a coarsening process where the initial form of the phase-separation into a coexistent axisymmetric core-shell geometry evolves over time. Within the plateau period of the equilibration dynamics, the solvent distribution evolves to minimize the elastic free energy cost due to coherency strain.

The curvature-polarization coupling term in (1ajbvbzckcrdcde) shows that the polarization direction is controlled by the curvature of the toroid, both the curvature ($\omega_{1},\omega_{2}$) at which it was fabricated and the change in curvature ($\Delta\omega_{1},\Delta\omega_{2}$) after deformation. A toroid fabricated with initial curvature $\kappa=1/R$ and without twist, so that the transverse material frame is defined by $\varphi={\rm constant}\equiv 0$, has an initial Darboux vector $\bm{\omega}=\kappa\mathbf{\hat{b}}$, where $\mathbf{\hat{b}}$ is the binormal of the centerline of the torus prior to deformation. The resulting change in free energy due to coupling between initial curvature and polarization is therefore

[TABLE]

where $\mathbf{\hat{n}}$ is the normal of the centerline. The free energy therefore decreases when $\mathbf{p}$ points along $-\mathbf{\hat{n}}$, describing a situation in which solvent mass is greater towards the outer radii of the torus and there is more polymer mass near the “hole” region; this is confirmed by experiment [see figure 16(f-j)]. This curvature-solvent distribution coupling has been observed in bent rubber and is due to the effect of internal stress on a material’s ability to swell [88]; swelling is promoted in regions under tension and impeded in regions under compression. Furthermore, with such a polarization of the solvent distribution, the similar coupling with $\Delta\omega_{m}$ implies that free energy is decreased if the torus deforms such that curvature increases. This coupling between bending deformations and the curvature of the centerline is due to an internal stress distribution, causing a torque about the centerline. A simplified picture of the situation is shown in 21(a), where the solvent-rich region is under compressive stress due to lamination to the solvent-poor region; if axial symmetry is broken, then this stress is centered along a region that is offset from the centerline. The effect is then similar to the classical problem of Timoshenko’s heated bimetallic strip [89], in which two metals with different thermal expansion coefficients are laminated together in a strip-like geometry; under heating, the constraint imposed by lamination results in a coherency strain that causes the strip to bend in the direction of the strip that expands less. Here, the gel bends in order to compress the solvent-poor region whilst expanding the solvent-rich region, thus increasing the curvature of the centerline; this swelling version of the bimetallic strip is ubiquitous in shape-changing soft materials [90, 91, 88, 16, 21]. We may therefore define a swelling moment $M_{m}=-k\epsilon_{mn3}p_{n}$ that characterizes the local internal torque that is applied transverse to the centerline, bending it. Since this moment is uniform around the central ring of the torus, it acts to uniformly increase the curvature $\kappa$.

However, a fundamental result of the differential geometry of curves (see e.g., [86]) implies that the curvature $\kappa$ of any planar closed curve, when integrated over the curve’s arclength, remains constant. Therefore, unless the magnitude and direction of the polarization changes, the only way for the toroid to deform such that the total curvature $\kappa$ increases is for the toroid to deform out-of-plane. In order to bend out-of-plane, however, the deformation must overcome the cost of bending and twisting. To study this buckling transition, we fix the magnitude and direction of the polarization $\mathbf{p}$ and consider only the elastic part of the free energy change, namely

[TABLE]

where we fix $\mathbf{p}$ along the $-\mathbf{\hat{n}}$ direction, which then fixes the swelling moment $\mathbf{M}$ along $-\mathbf{\hat{b}}$.

To determine the critical swelling moment $M_{c}$ needed to buckle the toroid from its planar configuration, parametrize a deformed centerline $\bm{\gamma^{\prime}}$ by

[TABLE]

where $\bm{\gamma}(s)$ is the original centerline

[TABLE]

which has a corresponding Frenet-Serret basis

[TABLE]

as pictured in the inset of 21(b). To second order in the out-of-plane deflection $\zeta$, the Frenet-Serret basis of the deformed centerline is

[TABLE]

from which we can define the deformed material frame $\{\mathbf{\hat{d}^{\prime}_{1}},\mathbf{\hat{d}^{\prime}_{2}}\}$ via a rotation by $\varphi^{\prime}(s)$ as shown in (1ajbvbzcka-1ajbvbzckc). To second order in $\zeta$ and $\varphi^{\prime}$, the elastic part of the free energy change is given by

[TABLE]

which may be simplified via the substitution555Since there is an ambiguity in how to define the material frame, the field $\varphi$ represents a gauge degree of freedom of the framed curve and this substitution is a gauge transformation. $\varphi^{\prime}=-\partial_{ss}\zeta-\partial_{s}\zeta+\tilde{\varphi}^{\prime}$ resulting in a transformed form of the free energy change

[TABLE]

Next, the two perturbing fields $\zeta(s)$ and $\tilde{\varphi}^{\prime}(s)$ can be expanded in Fourier modes,

[TABLE]

which diagonalizes the free energy change, yielding a quadratic form

[TABLE]

with

[TABLE]

The planar torus is therefore unstable to buckling out-of-plane when the free energy change of a certain mode $n$ becomes negative. This stability threshold occurs when ${\rm det}\,\mathcal{A}_{n}=0$, i.e.,

[TABLE]

which results in a critical swelling moment $M^{(n)}_{c}$ for each mode $n$ as a function of bending modulus $B$, twisting modulus $C$, and initial curvature $\kappa$, as plotted in figure 21(b) [66]. For set values of bending and twisting moduli and initial curvature, there is a finite value of $M$, above which the planar toroid is unstable to buckling out of the plane. The first mode that becomes unstable is the $n=2$, corresponding to a saddle or PringleTM-like morphology; higher modes become unstable for larger values of $M$. While the values of the effective bending and twisting moduli depend on factors such as the composition of the torus and the thickness of the shell, the ratio $C/B$ is well-approximated by the result for a uniform elastic rod with circular cross-section, namely $C/B\approx 1/(1+\nu)$, where $\nu$ is the Poisson ratio of the elastic material. Since the gel at fixed volume fraction is similar to an incompressible rubber, we take $\nu=1/2$, yielding $C/B\approx 2/3$, which marks the lower limit of $C/B$ for rods of circular cross-section, as predicted by classical elasticity theory [81]. Note that threshold value of $M$ for buckling decreases as the curvature $\kappa$ of the torus decreases at fixed bending modulus $B$. This is completely analogous to the Euler buckling prediction of smaller critical compression $T_{c}$ for longer rods at fixed bending modulus. Given a simple estimation of $M/(B\kappa)$, we have found that the predicted buckling threshold at $C/B\approx 2/3$ agrees with experiments, supporting the phase-separated ring model [66].

The form of the effective elastic energy $\Delta F$ given in equation (1ajbvbzckcrdcde) hints at a description of the phase-separated gel in terms of Landau theory, namely

[TABLE]

where $C_{p}$ and $u$ are positive coefficients that stabilize the spatial variations and magnitude of the polarization order parameter $\mathbf{p}$. Note that in the case where the gel is constrained to lie straight, i.e. where $\omega_{m}=0$, the equilibrium configuration of the solvent polarization field $\mathbf{p}$ is ordered in a ferromagnetic arrangement, with fixed magnitude and a spontaneously selected alignment direction transverse to the centerline of the gel. Consequently, uniform rotations of this alignment direction do not increase the free energy and thus disturbances in the configuration of the gel, such a material inhomogeneity or an interruption in the gel’s uniform shape, can easily cause long-wavelength modulation of the alignment direction. These Nambu-Goldstone modes, which in the context of the ferromagnetic order that we expect of the polarization field, are similar to spin waves in ferromagnetic materials [32, 92]. If the gel is then allowed to bend in response to the polarized solvent distribution, re-introducing the coupling between centerline shape and solvent distribution, the ferromagnetic order causes a uniform bending moment, resulting in a uniformly curved gel ring; again, the direction of the ring curvature is spontaneously selected, much like the magnetic field of a ferromagnetic material in the absence of an externally applied field. As an aside, note that our discussion of toroidal gels carries through here: if the gel was initially formed in a ring shape, then the manufactured curvature acts as an applied field, aligning the solvent polarization vector in a preferred direction. The spin-wave excitations of the polarization field have a rather interesting consequence for the shape of the gel. A long-wavelength rotation of the polarization field causes a similarly long-wavelength rotation of the curvature direction. This means that the shape of the centerline is no longer confined to a plane; it adopts a helical shape rather than a circular one. In other words, the gel has a soft torsion mode so that the Nambu-Goldstone modes of this theory are perhaps better referred to as twist or torsion waves.

Aside from a theoretical curiosity, this observation suggest some potentially interesting experiments. Using polymer gel printing techniques [66] it is possible to create gel samples with a variety of different shapes. Continuing our analogy with ferromagnetic materials, a gel drawn in the shape of the letter “$S$,” as illustrated in figure 22, should, under rapid heating, adopt a polarization field that undergoes a reversal in direction due to the flip in the curvature direction. Since the two arcs of the “$S$”-shape are joined by a straight segment, our model suggests that the polarization should undergo a $\pi$ rotation, actuating an out-of-plane twist of the $S$-shape. This is an investigation for the future. For now, we conclude there are remarkable similarities between the solvent-stress coupling in the phase-separated polymer gel and magnetoelastic effects, which are currently being investigated for use in so-called “shape-programmable magnetic soft matter” [93], including recent work with flexible ferromagnetic rings [94].

6 Conclusions & Outlook

We have discussed a small subset of the rich array of phenomena that polymer gels can exhibit. Starting with a discussion of the derivation and assumptions that are present in the Flory-Rehner model of the equation of state for isotropic polymer gels, we have attempted to discuss some of the swelling behavior of gels in a manner that is model-agnostic, only using the Flory-Rehner model to illustrate certain predictions. In particular, we have highlighted the swollen-deswollen phase transition, situating it within the classical theory of phase transitions of fluids. Departing from the traditional discussion surrounding this topic, we have paid particular attention on how intuition derived from the phase behavior of fluids fails, particularly at the critical point and along the phase coexistence curve, due to shear rigidity. In a spectacular departure from quasistatic processes, we have seen that rapid quenches across the first-order swelling transition can lead to arrested deswelling, trapping the gel in the coexistence region for a prolonged period, where it reaches an equilibrium state of coexisting phases. We have shown that the shapes adopted by toroidal gels when forced to coexist in this manner are dramatically distinct from the shapes adopted in quasistatic processes. By developing a description of the rapidly-heated gel that couples the spatial distribution of solvent within the gel to its elasticity, we have shown that the observed buckling arises from phase coexistence, and is thus linked to thermodynamic instability of single-phase gel. This demonstrates that, in the context of polymer gels, thermodynamic instability can be used to achieve a shape change that cannot be normally accessed in the thermodynamically stable regime. There are other interesting changes in material properties that can occur due to thermodynamic instability, such as possible auxetic behavior [95] and microstructure formation [49, 96]. We have conjectured that such instability may be used as part of material design, an idea that we describe as extreme thermodynamics.

Recently, there has been interest in designing equilibrium shape change in polymer materials [97, 98]. In one approach, the swelling response of the bulk gel is tuned by spatially modulating the density of cross-links or the gel’s chemistry, which impacts both the elasticity of the gel and its equilibrium volume fraction at constant temperature [75, 99, 100, 101]. Since different portions of the gel equilibrate to different volume fractions, the result is an internal stress distribution due to coherency strain that frustrates the original shape of the gel, leading to interesting shape change. Contrasting this design of the equilibrium gel shape, it has been shown [102] that if the gel is taken out of equilibrium by exposing different parts to different temperatures or by only exposing part of the gel to solvent whilst keeping other parts dry, the gel undergoes a dramatic set of shape transformations.

Another approach that yields the ability to design complex geometries from initially planar gels has its roots in certain processes observed in nature [103, 104]. It has been observed that pine cones are able to actuate shape change in response to changes in humidity, opening and closing their scales. The reason for this is a bilayer structure built from plant tissue that yields an anisotropic response upon swelling [90]. Following this realization, it was found that certain seed pods [105] split open from a flat state, forming two helical halves of opposite chirality, via an anisotropic shrinking process where two layers of tissue shrink in different directions, changing the intrinsic curvature of the seedpod. This layered anisotropy has been adopted for use in a novel additive manufacturing technique [98, 21], where polymer gel, made anisotropic through the use of aligned cellulose fibers within the gel, is printed in layers of different swelling-direction. When swelling is actuated by immersion in a solvent, these gel structures undergo morphological evolution that mimics natural processes, such as the opening of orchids.

We suggest that the shape changes accessed by rapidly heating gels through their phase transition may be considered in the context of these examples of designed shape change. As we have shown, the solvent distribution within the torus that is brought to a state of phase-coexistent equilibrium is set by the curvature of the toroidal centerline. While we have focused our attention on toroidal gels with a single curvature, the form of the solvent distribution polarization-curvature coupling in the free energy change (1ajbvbzckcrdcde) shows that the solvent polarization direction can be guided by local curvature of the centerline. We therefore conjecture that for more general polymer gel rings, where the curvature can vary continuously along the centerline, the solvent distribution will be polarized according to the local curvature direction. The result is that after rapid heating, these rings should deform in a manner that increases the magnitude of the local curvature. Moreover, in the case where the curvature direction undergoes a rapid reversal, such as in the letter “$S$,” we speculate that in order to interpolate between the opposite curvature directions, the polarization vector will rotate, actuating a twist of the gel. It is worth noting that the interesting deformations observed in experiments on tori require a simple actuation, namely rapid heating, without any prior patterning of the gel: the only feature that guides shape change is the initial shape of the gel. Thus, a single toroid can undergo at least two very different types of shape change depending only on heating protocol, namely isotropic deswelling under slow heating and buckling under rapid heating. This provides access to a much richer array of possible material responses that extends beyond the current regime of prescribed buckling, leading to a possibility of feedback between the material shape, its phase, and its response to applied stress.

Naturally, there is much more that can be explored regarding the physics of polymer gels. For example, we have neglected the topic of polymer gel dynamics entirely. Continuum hydrodynamic models have been developed, based on small deviations from equilibrium conditions for the gel [106, 107, 67, 108, 109, 48]. There have been descriptions of the coarsening of gels that have undergone spinodal decomposition [96] as well as of the development of surface patterns [110, 111, 65, 99]. However, many of these models either represent important idealizations of the gels or are extremely complex, requiring considerable computational resources. Furthermore, even robust hydrodynamic descriptions based on the Flory-Rehner model may not adequately describe the kinetics of the phase-transition, due to the limitations of the model. In particular, the equilibration kinetics that arise in quench experiments on polymer gels represent a considerable challenge for dynamical studies, due to the multiple timescales involved. For example, it is difficult to determine or predict the thickness of the solvent-poor skin as a function of time after the quench; this information is essential for a full understanding of the process as it determines the flow rate of solvent out of the swollen interior. It is not clear that the Flory-Rehner model, which provides a clear description of a homogeneous gel in or near equilibrium, is the appropriate equation of state for studying the dynamics in this regime; rather, a separate kinetic description may be necessary. In addition, the development of patterns, such as balloon and bamboo-like structures that appeared in both cylindrical [12, 65] and toroidal gels during deswelling, as shown in figure 16(j), is complicated by the necessity of a full nonlinear elasticity description of the gel that is able to incorporate large strains.

Acknowledgments

This work was supported by the National Science Foundation (DMR-1609841).