# Computer simulations of the mechanism of thickness selection in polymer crystals

## I Introduction

In 1957 Andrew Keller reported that polyethylene formed chain-folded lamellar crystals from solution. This discovery was followed by the confirmation of the generality of this morphology—lamellar crystals are formed on crystallization from both solution and the melt for a wide variety of polymers—and the basic phenomenological laws describing such properties as the thickness and growth rate. In particular, the crystal thickness, $`l`$, has been found to be inversely proportional to the supercooling, which is interpreted as resulting from $`l`$ being slightly larger than $`l_{\mathrm{min}}`$, the minimum thickness for which a lamellar crystal is stable with respect to the solution or melt, i.e. $`l=l_{\mathrm{min}}+\delta l`$, where $`\delta l`$ is small.

Surprisingly, however, no theoretical consensus has yet been reached as to the mechanism of this seemingly simple behaviour. In particular, two of the most well-known theories—the Lauritzen-Hoffman (LH) surface nucleation theory and the Sadler-Gilmer (SG) entropic barrier model—present very different explanations of thickness selection. Of course, in such a situation, one would like to determine which of the theories, if any, is closest to the truth. There are two aspects to such a task. Firstly, the predictions of the theories should be critically compared with experimental results. In the case of polymer crystallization both the LH and SG theories are able to reproduce the basic behaviour: the observed temperature dependence of the thickness and the growth rate. Additionally, Hoffman and coworkers have further developed the surface nucleation approach in order to explain some of the more detailed behaviour of crystallizing polymers, for example the regime transitions in the growth rate. However, this comparison does not conclusively favour one of the theories. This situation illustrates the fact that although consistency with experiment is an important first hurdle for any theory, it does not automatically imply the correctness of a theory. There may be a number of different ways of generating a particular experimental law. Furthermore, the number of parameters in a complex theory may give the theory sufficient plasticity to fit a wide variety of scenarios.

Secondly, it is important that the assumptions of a theory, particularly those about the microscopic mechanisms, are critically examined. However, in the case of polymer crystallization this task is very difficult to achieve experimentally. By addressing this gap, computer simulations can potentially play an important role in this field. Such simulations could range from examining simple models to performing realistic atomistic simulations of the crystal growth process. The former could allow the effects of relaxing some of the theoretical assumptions to be determined and the latter could provide a detailed molecular picture of the growth process. Indeed, there has been an increasing number of computational studies pursuing these aims. In this paper I will review my efforts in this direction and hope to illustrate the positive role that computer simulations can play in helping to understand polymer crystallization. In particular, the aim of my simulations has been to critically examine the LH and SG theories.

## II Free Energy profiles

In the LH theory the growth of a new layer is modelled as the deposition of a succession of stems (straight sections of the chain that traverse the growth face) along the growth face from an initial nucleus, where the length of each stem is the same as the thickness of the lamella. The inset of Figure 1 illustrates the geometry of this mechanism. To analyse the kinetics of growth, a thermodynamic description of the nucleation and growth of a new layer is first required. The free energy of a configuration with $`N_{\mathrm{stem}}`$ complete stems is taken to be

$$A(N_{\mathrm{stem}})=2bl\sigma +2(N_{\mathrm{stem}}1)ab\sigma _fN_{\mathrm{stem}}abl\mathrm{\Delta }F,$$
(1)   
where $`a`$ and $`b`$ are the width and depth of a stem, $`l`$ is the thickness of the lamella, $`\sigma `$ is the lateral surface free energy, $`\sigma _f`$ is the fold surface free energy, and $`\mathrm{\Delta }F`$ is the free energy of crystallization. The first term corresponds to the free energy of the two lateral surfaces created on the deposition of the first stem and is proportional to $`l`$. The second term is the free energy of the new fold surface created on the deposition of subsequent stems. It is then assumed that at the barrier between configurations with different numbers of stems all the new surfaces have been created and that a fraction $`\mathrm{\Psi }`$ of the free energy of crystallization is released. This then gives the LH free energy profile that is illustrated in Figure 1.

From this free energy profile, $`S(l)`$, the flux over the barrier, can be obtained. The observed crystal thickness is then taken to correspond to the average

$$\overline{l}=_{l_{\mathrm{min}}}^{\mathrm{}}lS(l)𝑑l.$$
(2)   
This average thickness is close to the value of $`l`$ at the maximum in $`S(l)`$, which in turn is close to, but slightly above $`l_{\mathrm{min}}`$, thus reproducing the observed behaviour of $`l`$. The maximum in $`S(l)`$ is the result of two competing factors. The free energy barrier for deposition of the first stem increases with $`l`$, thus making the growth of thick crystals prohibitively slow. However, as $`l_{\mathrm{min}}`$ is approached from above, the thermodynamic driving force for crystallization goes to zero.

It is important to note that by integrating over $`l`$, Equation (2) assumes that there are crystals with all values of $`l`$ greater than $`l_{\mathrm{min}}`$ which all grow with constant thickness and contribute to the average $`\overline{l}`$. Those crystals with a thickness close to the maximum in $`S(l)`$ dominate this ensemble and contribute more to Equation (2) because of their rapid growth. As was realized by Frank and Tosi, the results of experiments where the temperature is changed during crystallization argue against such an ensemble. The temperature jumps give rise to steps on the lamellae, showing that a crystal need not necessarily grow at constant thickness.

We will come back to this issue later, but in this section we focus on the LH free energy profile. In particular, we compare this theoretical profile with ones computed from simulations of a simple polymer. In our model the polymer is represented by a self-avoiding walk on a simple cubic lattice. There is an attractive energy, -$`ϵ`$, between non-bonded polymer units on adjacent lattice sites and between polymer units and the surface, and an energetic penalty, $`ϵ_g`$, for kinks (or ‘gauche bonds’) in the chain. The parameter $`ϵ_g`$ determines the stiffness of the chains. In our simulations we have included a surface which represents the growth face of a polymer crystal.

To follow the crystallization of the polymer on the surface, we need to define an order parameter which determines the degree of crystallinity. We use $`N_{\mathrm{xtal}}`$, the largest fragment of the polymer with the structure of a target crystalline configuration. In our case, we examine the crystallization of a 200-unit chain into a structure with 5 stems of length 40 units. In order to compare with the theoretical profiles we have to constrain the other $`NN_{\mathrm{xtal}}`$ units in the chain to be disordered. The simulations were carried out using configurational-bias Monte Carlo, and the umbrella sampling technique was used to calculate the free energy profiles.

The free energy profiles that we obtained are shown in Figure 2a. They show the expected temperature dependence: at low temperature the crystal is most stable and at high temperature the disordered state is most stable. Note that the value of $`N_{\mathrm{xtal}}`$ for the disordered state is non-zero, because the disordered polymer is adsorbed on the surface. The adsorbed polymer is bound to have some short straight sections that qualify as crystalline by the definition of $`N_{\mathrm{xtal}}`$. The free energy profiles also have a sawtooth structure resembling that of the theoretical profile. The barriers occur immediately after the previous stem has been completed, and correspond to the formation of a new fold. They are followed by a monotonic decrease in energy as this new stem grows to completion. In the language of the LH theory $`\mathrm{\Psi }(N_{\mathrm{stem}}N_{\mathrm{stem}}+1)0`$ for $`N_{\mathrm{stem}}2`$. However, there is no feature in the simulation profiles that corresponds to the formation of the first fold. This is because the initial nucleus is not a single stem, but two stems connected by a fold that grow simultaneously. Such a possibility had previously been suggested by Point.

Confirmation of a two-stem nucleus comes from a simple model calculation of the free energy profile. We can write the free energy as

$$A(N_{\mathrm{xtal}})=A_{\mathrm{coil}}(NN_{\mathrm{xtal}})+kT\mathrm{exp}\left(E_{\mathrm{xtal}}/kT\right),$$
(3)   
where the sum is over all possible crystalline configurations which are $`N_{\mathrm{xtal}}`$ units long, $`E_{\mathrm{xtal}}`$ is the energy of the crystalline configuration, and $`A_{\mathrm{coil}}`$ is the free energy of an ideal two-dimensional coil. The resulting profile is very similar to the simulation profile (Figure 2b). In particular, there is no feature due to the formation of the first fold. However, when we force the initial nucleus to be a single stem by restricting the sum in the above equation to only those crystalline configurations with one incomplete stem, a free energy barrier associated with the formation of the first fold appears. The reason for the preference for a two-stem nucleus is simply energetic. For $`N_{\mathrm{xtal}}>4ϵ_g/ϵ+2`$ the two-stem nucleus is lower in energy because of the interaction between the two stems. Our simulations were performed on a surface that was infinite. Whether a two-stem nucleus would be expected, when, as with a lamellar crystal, the thickness of the growth face is finite, depends upon how this critical size compares to the thickness of the lamella.

It can be clearly be seen from Figure 2b that the two-stem nucleus significantly reduces the nucleation barrier. In particular, it will no longer be proportional to $`l`$. This has significant implications for the LH theory given the key role played by this initial free energy barrier in constraining $`\overline{l}`$ to a value close to $`l_{\mathrm{min}}`$.

Before we move on we should make a number of comments. First, the polymer model is very simple, and although there is no obvious reason why the thermodynamic reasons behind the two-stem nucleus should not also apply to a real polymer, there may be factors that are not included in our model that come into play.

Second, the profiles reflect our choice of order parameter. As we monitor crystallization unit by unit, during the growth of the first two stems the lateral surface energy is paid for at the same time as the free energy of crystallization is released. Therefore, in this size range $`\mathrm{\Psi }`$ is effectively equal to 1, albeit with the possibility of a two-stem nucleus. Hoffman, however, advocates a $`\mathrm{\Psi }`$=0 version of the LH theory—it has the advantage that it avoids a ‘$`\delta l`$-catastrophe’ (a divergence of the lamellar thickness) at large supercooling—in which he postulates that prior to crystallization an aligned physisorbed state is formed that has lost its entropy but not yet gained the free energy of crystallization. Such a state cannot occur in our lattice model because there is no difference in the interaction with the surface for a disordered chain adsorbed on the surface and a crystalline layer. In an off-lattice model the interaction energy for the crystalline layer would be greater because the stems would fit into the grooves provided by the stems of the previous layer.

Third, a good order parameter must pass continuously through intermediate values when the system goes between two states. However, one can imagine a number of mechanisms by which this criterion for $`N_{\mathrm{xtal}}`$ is broken. For example, in a realistic simulation of the surface crystallization of a long alkane into a once-folded configuration, the chain first formed non-adjacent crystalline stems connected by a loose fold which then came together by the propagation of a defect through one of the stems. Another possibility that has been observed in simulations is the formation of crystallites in different portions of a chain that subsequently coalesce to form a single crystallite.

## III A multi-pathway model

In the previous section, in order to compare the LH free energy profile with those from simulation, we had to constrain the $`NN_{\mathrm{xtal}}`$ units not having the target structure to be disordered. If we had not done this, at temperatures where the crystal is most stable the rest of the chain would have formed a crystalline configuration with stem lengths different from the target configuration. This naturally raises questions about the LH assumption that the stems in a new layer must all have the same thickness as the previous layer. In this section, we examine the effects of relaxing some of the LH assumptions by studying a model in which the stems grow unit by unit and the length of a stem is unconstrained. We term it a multi-pathway model because it can take into account the many possible ways that a new crystalline layer can form.

This idea is not new. Frank and Tosi, Price and Lauritzen and Passaglia considered models where the stem length is not always constant, and Point, and DiMarzio and Guttman studied models where the stems could grow unit by unit. All these studies were performed at a time when computational resources were much less, so approximations and simplifications had to be made in order to render the models tractable. The natural way to solve such problems, though, is through the use of computational techniques, such as kinetic Monte Carlo. However, the only applications of computational methods to this problem were in a short note by Point and the continuation of this work in the PhD thesis of Dupire. Some of the results presented in these earlier studies are similar to those we report here.

In our model we grow a single new crystalline layer by the successive growth of stems across a surface that represents the growth face of a polymer crystal. The polymer interactions are the same as used in the previous section, and we only model the crystalline portion of the polymer explicitly—the rest is assumed to behave like an ideal coil. An example configuration is illustrated in Figure 3 along with possible changes of configuration. These changes can only occur at the ends of the crystalline portion, and are selected using the kinetic Monte Carlo algorithm, in which a move is chosen with a probability proportional to the rate for that process.

First, we shall examine the effect of the initial nucleus on the thickness of the layers grown. If the stem lengths are unconstrained and the initial nucleus is a single stem, one might imagine that one way of reducing the large initial free energy barrier in Figure 1 (and achieving faster initial growth) would be for the stem length to increase gradually to its average value as crystallization progresses. For this pathway, the lateral surface free energy is paid for ‘in installments’ rather than all initially. This is exactly what we observe when we force the initial nucleus to be a single stem by only allowing growth from one end of the crystalline portion of the chain (Figure 4a). When a double-stem nucleus is allowed the initial growth is very different because there is now no longer a large initial free energy barrier to circumvent. The most important thing to note from these results is that, contrary to the LH theory, the thickness of the inital nucleus does not determine the thickness of the layer. Further confirmation of this can be obtained when we examine the growth from initial seed crystals. Whatever the thickness of the initial seed the thickness of the growing crystal converges to the same value (Figure 4b). This implies that the thickness of a crystalline layer must be determined by factors which are operating on the deposition of each stem and not those specific to the initial stems.

To determine what these factors might be, in Figure 5 we show how the thickness of a new layer depends on temperature. First, it is immediately obvious that the thickness of a new layer is not necessarily the same as that of the growth face. Second, all the curves increase as the temperature approaches $`T_m`$, the melting or dissolution temperature, because of the rise of $`l_{\mathrm{min}}`$. Third, the thickness also increases at low temperature, in this instance because it becomes increasingly difficult to scale the free energy barrier for forming a fold and so on average the stems continue to grow for longer. However, this rise is checked by the thickness of the growth face. It is unfavourable for the polymer to overhang the edge of the growth face because these units do not interact with the surface.

Figure 5 only describes the growth of a single layer. However, as the thickness of the new layer is not generally the same as the thickness of the growth face, one needs to consider the addition of a succession of layers. If we assume that all the variations in the stem length within a layer are annealed out before a new layer begins to grow, this can be achieved using Figure 6a, in which we have plotted for a single temperature the thickness of the new layer against the thickness of the growth face. By following the dotted lines one can see what would happen for growth on a growth face that is 50 units thick: the first layer is 36 units thick, the second 28, the third 23, …Thus, the thickness converges to the value $`l^{}`$ at which the curve crosses $`y=x`$, i.e. to the point where the thickness of the new layer is the same as the previous, and then the crystal continues to grow at this thickness. The mapping represented in Figure 6a is a fixed-point attractor.

A similar picture emerges if we explicitly perform simulations of multi-layer growth. Figure 7 shows a cut through a typical configuration that results. Within 5–10 layers the thickness of the crystal converges to its steady-state value $`l^{}`$ and then growth continues at that thickness.

The mechanism of thickness selection that occurs in our multi-pathway is at odds with the LH theory. It shows that it is inappropriate to compare the growth rates of crystals of different thickness because the thickness has only one dynamically stable value for which growth at constant thickness occurs. The ensemble of crystals assumed by Equation (2) is fictitious. Furthermore, the growth rate of a new layer slows down as $`l^{}`$ is approached from above (inset of Figure 6a). However, we should note that in some of the multiple-pathway studies mentioned earlier, it was realized that stable growth can only occur at the one thickness where a new layer has the same thickness as the previous. Since then this insight has for the most part been neglected.

To analyse the reasons for the dynamical convergence of the thickness to the value $`l^{}`$ we examine how the probability distributions for the stem length depend on the thickness of the growth face (Figure 8). $`l_{\mathrm{min}}`$ places one constraint on the stem length; only a small fraction of the stems can be shorter than $`l_{\mathrm{min}}`$ if the layer is to be thermodynamically stable. The thickness of the growth face places the second constraint on the stem length; it is energetically unfavourable for the polymer to extend beyond the edges of the growth face. There is also a third weaker kinetic constraint on the stem length. At every step there is always a certain probability that a fold will be formed. Therefore, even in the absence of the second constraint, i.e. an infinitely thick growth face, the probability distribution will decay exponentially to zero at large stem length (Figure 8a). Although, this effect prevents the thickness from ever diverging in a $`\delta l`$-catastrophe, it does not stop the thickness becoming very large.

When the growth face is significantly thicker than $`l_{\mathrm{min}}`$ there is a range of stem lengths between $`l_{\mathrm{min}}`$ and the thickness of the growth face that are viable, and therefore the new layer will be thinner than the previous layer. However, as the thickness of the growth face decreases, the probability distributions of the stem length becomes increasingly narrow and the difference in probability between the stem length being greater or less than the surface thickness diminishes. Finally, at $`l^{}`$, as the thickness of the growth face approaches $`l_{\mathrm{min}}`$, the probability distribution become symmetrical about the surface thickness and the thickness of the new layer becomes equal to the thickness of the growth face (Figure 8e). When the thickness is less than $`l^{}`$, the asymmetry of the probability distribution is reversed (Figure 8f). It is, therefore, through the combined action of the two thermodynamic constraints on the stem length that the thickness converges to a value close to $`l_{\mathrm{min}}`$.

The picture is not quite this simple at all temperatures. As the supercooling decreases, it becomes increasingly unfavourable for a stem to overhang the edge of the growth face. Indeed, for sufficiently small supercoolings the probability distribution for the stem length never becomes symmetrical about the thickness of the growth face, not even when the thickness of the growth face is close to $`l_{\mathrm{min}}`$. This situation is illustrated in Figure 6b. After the growth of two layers on a 50-unit thick surface, the crystal stops growing because the outer layer is too thin for a new layer to form. For these supercoolings, as in the SG model, the rounding of the crystal profile inhibits growth. To overcome this barrier requires a cooperative mechanism whereby a new layer takes advantage of (and then locks in) dynamic fluctuations in the outer layer to larger thickness. However, unlike the SG model, the current model has no interlayer dynamics—we attempt to grow a new layer on an outer layer that is static—and so growth stops. Despite this it is clear that if this interlayer dynamics could be included, it would again lead to steady-state growth close to $`l_{\mathrm{min}}`$.

We should note that this cessation of growth was also found in the model of Frank and Tosi at low supercoolings. Lauritzen and Passaglia were also aware of this effect, but they introduced an ad hoc energetic term in their rate constants to prevent it. However, in the restricted equilibrium model of Price this effect was absent. In this study each new layer, but not the crystal as a whole, was allowed to reach equilibrium and so the kinetic constraint on the stem length is absent.

Finally, we should note that our multi-pathway model is not parameter-free, and that, like most other models of polymer crystallization (including the LH theory, the SG model and the earlier multi-pathway models), for some choices of parameters (not those used here) the lamellar thickness begins to increase at sufficiently large supercooling. This effect occurs because the large driving force for crystallization at large supercoolings reduces the effect that the thickness of the growth face has in constraining the stem lengths.

## IV The Sadler-Gilmer model

In this section we re-examine the model used by Sadler and Gilmer in order to see whether the mechanism of thickness selection that we found in the previous section for our multi-pathway model also occurs in the SG model. Sadler and Gilmer interpreted this model in terms of an entropic barrier. In particular, they argued that the rounding of the crystal profile gives rise to an entropic barrier, which can only be surmounted by a fluctuation to a squarer profile before growth can continue. As this barrier increases with lamellar thickness it constrains the thickness to a value close to $`l_{\mathrm{min}}`$. However, we shall not dwell on this interpretation here, but instead direct the interested reader to a critique of this argument in Ref. .

In the SG model the growth of a polymer crystal results from the attachment and detachment of polymer units at the growth face. The rules that govern the sites at which these processes can occur are designed to mimic the effects of the chain connectivity. In the original three-dimensional version of the model, under many conditions the growth face is rough and the correlations between stems in the direction parallel to the growth face are weak. Therefore, an even simpler two-dimensional version of the model was developed in which lateral correlations are neglected entirely, and only a slice through the polymer crystal perpendicular to the growth face is considered.

The geometry of the model is shown in Figure 9. Changes in configuration can only occur at the outermost stem and stems behind the growth face are ‘pinned’ because of the chain connectivity. At each step, there are three possible changes in configuration: the outermost stem can increase in length, a new stem can be initiated and a polymer unit can be removed from the outermost stem. The model can be formulated in terms of a set of rate equations that can be easily solved by numerical integration.

When we examine the dependence of the thickness of a layer on the previous, we again find a fixed-point attractor describing the convergence of the thickness to its steady-state value (Figure 10). Moreover, when we examine the probability distributions for the stem length we find evidence for the same three constraints as for the multi-pathway model (Figure 11b). The weaker nature of the kinetic constraint is particularly clear from the much more rapid exponential decay of the probability for stems that extend beyond the growth face. The role played by the two thermodynamic constraints in the mechanism of thickness selection is particularly clear from Figure 11b. As the thickness of the growth face decreases the viable range of stem lengths decreases until the the thickness of the growth face meets $`l_{\mathrm{min}}`$ at the fixed point.

## V Discussion

In this paper we have outlined evidence from computer simulations for a mechanism of thickness selection in lamellar polymer crystal that differs from the theories of Lauritzen and Hoffman, and Sadler and Gilmer. Instead, the mechanism has much more in common with the results of earlier multi-pathway models. We find that a fixed-point attractor which describes the dynamical convergence of the crystal thickness to a value just larger than the minimum stable thickness, $`l_{\mathrm{min}}`$. This convergence arises from the combined effect of two constraints on the length of stems in a layer: it is unfavourable for a stem to be shorter than $`l_{\mathrm{min}}`$ and for a stem to overhang the edge of the previous layer. It is encouraging to note that we find the same mechanism of thickness selection operating in two models which make very different assumptions about the microscopic growth processes. This provides evidence of the generality of this mechanism, and so suggests that, although the models described here have a very simplified description of the microscopic dynamics, the physical principles behind the mechanism could be general enough to apply to real polymers.

This mechanism of thickness selection is also consistent with experiments where the temperature is changed during crystallization. The steps that result indicate that the thickness of the lamellar crystals dynamically converges to the steady-state thickness for the new temperature by a mechanism similar to that which we observe in our simulations. Furthermore, if the step profiles could be characterized with sufficient resolution by atomic-force microscopy, it may be possible to extract the fixed-point attractor of a real polymers. However, for a temperature decrease the step profiles may also reflect the rounding of the crystal edge and for a temperature increase the roughness of the fold surface. Furthermore, any annealing mechanisms that operate could change the shape of the step profile from its as-formed state.

Although the multi-pathway approach is, in some ways, an extension of the LH theory, the removal of many of the LH constraints leads to significantly different behaviour. In particular, our work undermines the LH assumptions that the initial nucleus determines the thickness of a layer, and shows that the approach embodied in Equation (2) (i.e. a comparison of the growth rates of the crystals in an ensemble of crystals of different thickness all of which grow at constant thickness) is inappropriate because crystals of arbitrary thickness do not necessarily continue to grow at that thickness. Although our results lead us to question the thickness selection mechanism in the LH theory, other aspects of the nucleation approach may not be affected by our critique. For example, the regime transitions are a result of the different functional dependence of the growth rate on the nucleation rate and the substrate completion rate in the different regimes.

Recently, there have been a number of alternative theoretical proposals that have made recourse to metastable phases. Keller and coworkers suggested that crystallization of polyethylene could initially occur into the mobile hexagonal phase. These crystals would then thicken until a critical thickness was reached at which a phase transition to the orthorhombic phase would occur. Olmsted et al. have argued that the density fluctuations resulting from the spinodal decomposition of a polymer melt assist the nucleation of crystals. Strobl and coworkers have argued, on the basis of the thickness dependence of the crystallization and melting temperatures of syndiotactic polypropylene, and the granular texture in AFM images of the same polymer, that the polymer first crystallizes into blocks, which are subsequently stabilized when they fuse into lamellae. Our simulations can say little about these proposals since our polymer models are too simple to be able to capture such features. However, all these approaches are based on behaviour that has been observed in crystallization from the melt, so it is not clear how the ideas can apply to crystallization from solution, where the same basic laws for lamellar polymer crystals apply.