Spatial distribution of entanglements in thin free-standing films
SSpatial distribution of entanglements in thin free-standing films
Daniel M. Sussman ∗ Department of Physics and Astronomy,University of Pennsylvania, 209 South 33rd Street,Philadelphia, Pennsylvania 19104, USA
Abstract
We simulate entangled linear polymers in free-standing thin film geometries where the confiningdimension is on the same scale or smaller than the bulk chain dimensions. We compare bothfilm-averaged and layer-resolved, spatially inhomogeneous measures of the polymer structure andentanglement network with theoretical models. We find that these properties are controlled by theratio of both chain- and entanglement-strand length scales to the film thickness. While the film-averaged entanglement properties can be accurately predicted, we identify outstanding challenges inunderstanding the spatially resolved character of the heterogeneities in the entanglement network,particularly when the scale of both the entanglement strand and the chain end-to-end vector iscomparable to or smaller than the film thickness. ∗ [email protected] a r X i v : . [ c ond - m a t . s o f t ] J u l . INTRODUCTION The properties of polymer melts under strong confinement conditions has been a subject ofintense interest, where much of this attention has focused on elucidating changes in the glassyproperties of nanometrically thin samples [1]. It is known, though, that strong confinementalso affects the properties of thin polymer films in the melt state. For instance, bubble-inflation experiments have revealed a stiffening in the creep compliance of thin films in therubbery plateau regime [2], a feature which is robust to performing the experiments wellabove the glass transition temperatures of the films [3]. Measurements in nanometricallyconfined thin films have reported lower viscosity [4, 5] and greater mobility [6], and thediffusion rate in, e.g., cylindrical nanopores is known to be affected by the relative degree ofchain confinement [7]. Recent work on polymer nanocomposites has revealed that the (local)confinement effects introduced by nanoparticles can lead to interesting chain conformationand entanglement heterogeneities in a sample[8–10], which can in turn affect the globalsystem properties.An emerging theme is the difficulty of understanding entangled polymer melts when thedimensions of the confining volume are comparable to or smaller than the typical end-to-enddistance of the chains, R ee . This difficulty is compounded in the “ultra-thin” films typicallysimulated, where the confining dimensions may be small or comparable not only to the chainconformational scale but also to the typical length scale associated with the entanglementstrands, L e . While techniques to predict chain-scale conformational changes are relativelywell-developed [11, 12], the connections between these chain or sub-chain structural changesand system dynamics are much more poorly understood.Experiments on entangled thin films have interpreted changes in macroscopic systemproperties in terms of an increase in the entanglement molecular weight, M e , and a decreasein interchain entanglement density in response to chain-scale confinement [13, 14]. Consis-tent with this observation, coarse-grained simulations analyzed using “primitive path” (PP)-determining methods directly observed a change in the number of interchain entanglementsunder both geometric [12, 15–17] and phase-induced confinement [18] conditions. Again,though, the connections between structural entanglement statistics – such as the number ofentanglements per chain ( (cid:104) Z (cid:105) ), the tube diameter ( d T ), and the number of monomers perentanglement ( N e ) – and the dynamical properties of the system are unclear. The classic2oi-Edwards model provides a natural framework for addressing this question in bulk, ho-mogeneous systems [19], but some modifications to the theory are needed under conditionsof strong confinement. For instance, recent simulations of polymers confined in nanoporessuggest that the changes in chain diffusion are qualitatively consistent with stronger confine-ment leading to lower (cid:104) Z (cid:105) , but that a straightforward reptation-model estimate of diffusivityseems quantitatively inaccurate [20].We take one possible interpretation of these findings to be that the primitive path mea-sures of entanglement in simulations indeed directly inform system dynamics, but that amodel that explicitly treats any spatial inhomogeneity of entanglement statistics must beused to accurately predict dynamical properties. This challenge is two-fold, since at themoment there are few theories that might be used to predict the inhomogeneities of anyentanglement statistics in the first place, let alone theories that also connect those inhomo-geneities to other measurable single-chain or macroscopic observables of the system. Oneappealing model of entanglements near a free surface simply posits that at the surface thereare only half as many chains with which to entangle. If there are (cid:104) Z (cid:105) inter-chain entangle-ments per chain for an equivalent bulk system, then, at the surface there might be (cid:104) Z (cid:105) / ∼ R ee from aboundary, and a recent theory has been proposed connecting the orientational correlationsthat arise from such conformational perturbations to entanglement statistics [22, 23]. Suchcorrelations could be relevant for not only confined chain systems, but also bulk chains sys-tems undergoing deformations, where it is known that chain-scale orientation occurs togetherwith a decrease in the effective entanglement number [24, 25].In this paper we perform simulations with the goal of collecting statistics on the averageand spatially resolved distribution of both chain conformational properties and entanglementstatistics of confined polymer chains. We hope that this data both provides a robust testof the above theories and encourages the development of other theoretical models of chainentanglement in confined and non-equilibrium polymer systems. We focus on preparingequilibrated free-standing thin films, using a standard bead-spring coarse graining, equi-libration protocol, and PP-extraction method, described in more detail below. We havechosen to work with monodisperse polymers composed of N = 2000 monomers, in films ofthickness ∼ , , ,
48, and 58 in units of the monomer diameter. We also present results3n shorter-chain films derived from these base data sets. These parameters were chosen toprobe regimes in which the film thickness is comparable to or larger than the bulk tube di-ameter, while being simultaneously comparable to or smaller than characteristic end-to-endchain scale.The remainder of the paper is organized as follows. Section II describes in detail thesimulation model used and the protocol followed to equilibrate the systems studied in thiswork. Section III presents the numerical results for the film-averaged and the spatiallyresolved chain entanglement statistics found in our free-standing films, as well as spatiallyresolved orientational correlations on both the chain end-to-end and PP length scales. Weclose in Section IV with a discussion of our results.
II. METHODS
We perform molecular dynamics simulations of a commonly used coarse-grained modelpolymer [26] prepared in thin film geometries. Our systems have N T = 5 × total particlesof diameter σ , composed of M = 250 monodisperse chains of length N = 2000. The non-bonded interactions between particles i and j are specified by V nbij = 4 (cid:15) (cid:32)(cid:18) σr ij (cid:19) − (cid:18) σr ij (cid:19) (cid:33) − (cid:15) (cid:32)(cid:18) σr c (cid:19) − (cid:18) σr c (cid:19) (cid:33) (1)for r ij < r c and V ij = 0 for r ij > r c . Here (cid:15) sets the energy scale, and we take the range ofthe non-bonded interactions to be r c = 2 .
5. For the bonded interactions we use a very stiffharmonic potential, V bij = k r ij − σ ) , (2)where k = 2000 (cid:15)/σ . In this paper we report results in reduced LJ simulation units, e.g.,the temperature T = kT ∗ /(cid:15) and the time τ = τ ∗ (cid:112) (cid:15)/mσ , where T ∗ and t ∗ are defined inlaboratory units and m is the mass of a particle. All of our simulations were run in the NVTensemble at T = 1 . δt = 0 . τ using the HOOMD-blue simulationpackage [27, 28] .Five base systems are considered in this work, corresponding to different effective filmthicknesses, h eff . The films were prepared using a common protocol, which we documenthere for completeness, given that for these large systems the overwhelming bulk of process-ing time is spent preparing independent equilibrated samples. We begin by pre-packing4on-interacting random walks with the correct single-chain statistics [29] in a simulationcell with the dimensions of the desired film. By treating two of the three dimensions asperiodic and the other ( z -) dimension as a reflecting wall at this step we are able to startwith a configuration of chains that will not be greatly perturbed when the z direction of thebox is expanded to expose the surfaces of the film to vacuum. We then use a collection ofchain-altering Monte Carlo (MC) moves that simultaneously reduce the density fluctuationsof the pre-packed chains[29] and respect the single-chain statistics imposed by the modelpotentials and by the thin film geometry [30]. From these non-interacting chain configura-tions the box was expanded in the z direction to allow for a free-standing film, the LJ andbonded interactions were slowly turned on, and molecular dynamics using the potentials andparameters above were begun. The film thickness, h eff , is determined by fitting the densityprofile of the films in the z direction, ρ ( z ), to ρ ( z ) = c (cid:20) (cid:18) | z | − ξ /λ (cid:19)(cid:21) . (3)This mean-field-like functional form fits the density profiles very well; here c is determinedby the overall density and the simulation box size, ξ is a measure of film thickness, and λ describes the extent of density fluctuations near the surface. These density fluctuationscombine the effect of equilibrium fluctuations whose scale depends on the bulk correlationlength and surface capillary waves whose contribution to the interfacial width would growlogarithmically with the size of the simulation box [31, 32]. For our films at T = 1 . λ ≈ h eff = 2( ξ + λ ), a spatial extent which for ourfilms encompasses at least 99% of the particles in the simulation box. By this measure, thefree-standing films simulated in this work have thicknesses h eff /σ = 13 . , , ,
48, and58.We are interested in the spatial distribution of entanglements in these simulated films,and we use the Z1 algorithm to extract the entanglement statistics of our polymer melts.For bulk systems the two geometric methods of extracting primitive paths (Z1 [35–38] andCReTA [39]) are expected to give similar global entanglement measures, and for confinedthin films this expectation has been explicitly verified [12]. We note that the primitive pathsrevealed by the Z1 algorithm vary slightly from the ideal random-walk statistics envisionedby the conventional tube model, with non-Gaussian distributions of primitive path segment5izes and weak angular correlations between consecutive primitive path steps. However, atthe level of theoretical comparisons that will be made in this paper we do not expect thenon-Gaussianities to be important.To speed up equilibration we implemented double-rebridging, connectivity altering moves[29, 33] on graphics processing units in the HOOMD-blue framework [30]. Modeled afterthe protocol outlined in Ref. [29], we check roughly half of the particles for chain-length-preserving bond swaps every 0 . τ ; this longer time leads to fewer proposed reverse swapscompared to checking every few time steps, and reduces the impact of the bond-swappingalgorithm on overall simulation performance. We note that we have expanded the algorithmto be able to handle not only bonds and angular potentials, but also dihedral potentials;results generated by treating chains with dihedral potentials will be discussed in a futurepublication.As an initial consideration for system equilibration, Auhl et al. noted that one needs O (10) pivoting moves to decorrelate the end-to-end vector of an isolated chain [29]. Inthe context of a dense melt, bond-swapping alone decorrelates subchains on which at least O (10) bond-swapping moves have been performed. To equilibrate the chains on all lengthscales, then, the molecular dynamics must be run up to the decorrelation time of the longestsubchains that have not participated in O (10) connectivity-altering Monte Carlo moves.Another consideration is that individual monomers should be able to sample bond swapswith a suitable number of exchange partners. For these linear chains the typical distancebetween monomers of the same index is ∼ ( ρ/ ( N/ − / ; if the monomers were to diffusewith Rouse-like dynamics the typical time to diffuse this distance, t s , corresponds to theRouse time of a chain of length N eff = ( N/ / [29]. For our chains of length N = 2000this corresponds to diffusive time scales of order t s ∼ τ . This use of Rouse-like dynamicsto estimate the relevant time scale arises because the Monte Carlo bond-swapping procedureallows chain reconfigurations that remove previously existing entanglement constraints. Thiseffect is seen in the monomer mean-square displacements, which under the action of thebond-swapping moves cross over directly from a 1 / / t eq = 7 . × τ , a time long enough that every monomeris involved in O (10) bond-swapping moves and long compared to the diffusive time for amonomer to travel the typical distance separating monomers of the same index. Since itis impractical to simulate the system for long enough that a typical monomer could diffusemore than the end-to-end distance of a chain, we check equilibration by (1) looking at thedistribution of single-chain conformational statistics [29] and (2) using the Z1 algorithm tolook at the entanglement statistics of our systems [34]. Figure 1 shows illustrative examplesof the time evolution of these measures during the equilibration phase. We note that theconformational comparison in Fig. 1 is to polymers using finite extensible non-linear elasticpotentials, V bij,F ENE = − . kR log (1 − ( r ij /R ) ) r ij ≤ R ∞ r ij > R , (4)with parameters k = 30 (cid:15)/σ and R = 1 . σ , rather than the stiff harmonic springs usedin this work. Additionally, to treat cohesive thin films we use a LJ potential with bothrepulsive and attractive regions, whereas Auhl et al. use a purely repulsive LJ potential forthe non-bonded interactions.. These differences lead to a slight change the expected bondlength (cid:104) b (cid:105) / for the single chains, here (cid:104) b (cid:105) / = 0 . σ whereas the system of Auhl etal. has (cid:104) b (cid:105) / = 0 . σ at the same density and temperature [29]. The angular statistics areidentical, though, and on the scale of the plot the effect of this difference for the mean-squareinternal distance statistics is negligible.For the entanglement statistics, in our equilibrium bulk N = 2000 systems of monomerdensity ρ s = 0 .
86 we find that the primitive path segment size, a pp ≡ (cid:104) R ee (cid:105) /L pp where L pp is the primitive path contour length, is a pp /σ = 12 . ± .
4, the mean number of kinksper chain is (cid:104) Z (cid:105) = 40 . ± .
20, and N e = 49 . ± .
28. These numbers are consistent withvalues reported in the literature for this model polymer [34]. In the remainder of the paperwe will take the PP segment size as our characteristic scale of the entanglement strands, L e ≡ a pp,bulk ; again, the slight non-Gaussianities of the primitive paths imply that theselength scales are, in principle, close but distinct quantities.Once the films are equilibrated as described above, we continue to run the moleculardynamics with connectivity-altering MC moves. We record and analyze our system every5000 τ , which corresponds to each chain experiencing ≈
100 successful bond swaps. By the7 · · · · /τ < Z > < R ( n ) > / n FIG. 1. Average number of kinks per chain measured by the Z1 algorithm during the equilibrationof systems with M = 250 chains each of length N = 2000. From top to bottom the curvescorrespond to a bulk system followed by thin films of size h eff = 35 σ , h eff = 20 σ , and h eff = 13 . σ .Inset. The mean-square internal distance along a chain [29] as a function of monomer indexseparation along a chain. Shown are curves for the equilibration of a bulk system at (top tobottom, dotted curves) t/τ = 0 , . × , . × , . × , . × , . × . The solidcurve is the “target function” from Auhl et al.[29] criteria described above our N = 2000 chains are still correlated: on average subchainsof length O (500) and smaller have yet to experience O (10) bond swaps, and since 5000 τ is much less than the Rouse time of N = 500 chains the molecular dynamics does notequilibrate these sub-chains. Additionally, individual monomers have not diffused enoughto sample a reasonable number of swap partners. By these criteria a better choice wouldbe to perform our analysis on frames separated by at least 4 × τ . This is greater than t s , and on this time scale bond-swapping would decorrelate the subchains up to length O (120), for which the Rouse time is only ∼ × . Even less frequent sampling wouldfurther decrease the correlations between the samples. However, the primitive path networkfound in the simulations is very sensitive to both the global chain features (that is, a largechange in the positions of the chain ends can affect entanglements even in the middle ofthat chain) and even a small number of connectivity altering swaps along subchains of size N e . Thus, we find that the primitive paths are decorrelated even on the shorter time scaleat which we save our configurations. Furthermore, by separately averaging our results overall saved frames and over only a subset, we have checked that our entanglement results arequalitatively insensitive to the fact that we are not sampling truly independent underlying8onfigurations. III. RESULTSA. Global entanglement statistics
With the film configurations described above in hand, we first report the film-averagedentanglement statistics of our systems. Using the Z1 algorithm, we find that our free-standing thin films indeed show a smaller number of entanglements per chain, (cid:104) Z (cid:105) , than ourbulk system. Figure 2 shows this reduction as a function of the film confinement parameter δ ≡ h eff /R ee,bulk , the effective film thickness divided by the bulk end-to-end chain distance,that has previously been used to collapse similar data [12, 13]. The figure also shows atheoretical prediction (described below) and simulation data from Ref. 12.To increase the number of data points we can investigate in this representation, wenote that we can subdivide the chains in our simulations of length N = 2000 polymersinto subchains of shorter length. For instance, severing the bonds in the center of eachchain creates a system of M = 500 chains of length N = 1000. Since the N = 2000chains have correct intermediate internal separation single-chain statistics, this derived sub-chain system will be in a nearly equilibrium configuration (excluding chain-end effects). Wenumerical test this assumption of equilibration in the same way as for the full chain systems,i.e., by comparing single-chain conformational statistics and by looking at the number ofentanglements per chain detected by the Z1 algorithm. We have confirmed that both themean-square internal distance statistics of single chains and the entanglement statisticsobtained in this way are indeed representative of the entanglement statistics shorter-chainsystems by independently equilibrating films of thickness h eff = 35 σ and h eff = 20 σ withthe same number of total monomers ( N T = 5 × ) grouped in chains of length N = 500.We confirm that the reduction in entanglement density is the same as that obtained byrunning the Z1 algorithm on systems created by the subdivision method. For the many othersubdivided chain systems considered below, we do not repeat the process of independentlyequilibrating films of the same chain length; we confirm that the single-chain conformationalstatistics fall on the equilibrium master curve, and take the entanlement data from runningthe Z1 algorithm on the sub-chain systems to represent the equilibrated value.9 .0 0.5 1.0 1.5 2.0 2.5 3.00.700.750.800.850.900.951.00 δ〈 Z 〉〈 Z 〉 bulk FIG. 2. System-averaged reduction of entanglements as a function of confinement. The largediamonds correspond to free-standing films of N = 2000 chains. Small diamonds correspond tofilms generated by taking equilibrium configurations of the N = 2000 chains and subdividing eachchain into subchains of length N = 1000 , , , , , , N = 350 chains confined between hard amorphous walls from Ref. [12]. The solid curveis the theoretical prediction of Ref. [12] Using the subdivision method allows us to infer the entanglement statistics in films ofthe same thickness but for chains whose length is (approximately) an integer divisor of N = 2000. Although the lack of collapse of the N = 2000 chains simulated here with thedata from Ref. [12] could be potentially attributed to the difference between freestandingand capped thin films, it is clear that the entanglement statistics of these subdivided chainsdo not collapse cleanly when the film thickness is reported in units of δ . Our independentsimulations of N = 500 chains confirms that this lack of collapse cannot be attributed toan artifact arising from studying subdivided-chain systems derived from equilibrated filmsof longer chains.One theoretical expectation for collapse of thin film entanglement data with δ , and thesource of the theoretical prediction in Fig. 2, arises from a statistical mechanical frameworkoriginally developed for fluids of entangled needles [40]. The extension of this theory toflexible chains argues that orientational ordering on the chain scale could lead to a largertube diameter and hence a reduction in the number of entanglements per chain [22, 23, 41].That theory assumes that orientational correlations present on the chain end-to-end scalewere directly communicated to the PP scale, and that PP orientational correlations, averaged10ver the entire system, affected the tube diameter as d T d T,bulk ∝ (cid:18)(cid:90) d(cid:126)u d(cid:126)u g ( (cid:126)u ) g ( (cid:126)u ) (cid:113) − ( (cid:126)u · (cid:126)u ) (cid:19) − , (5)where g ( (cid:126)u i ) is the orientational distribution of PP i . We note that the non-Gaussianities ofthe primitive paths revealed in simulations were shown to be reasonably well captured bythe theory that leads to Eq. refeq:theory [22]. Since it is possible to theoretically predictthe conformational properties of an ideal melt of chains in the presence of hard, reflectingwalls [11, 12], Eq. 5 can be used to predict from first principles the change in the tubediameter as a function of confinement. By assuming standard Doi-Edwards relations amongthe properties of the tube, one obtains (cid:104) Z ( δ ) (cid:105)(cid:104) Z (cid:105) bulk = (cid:104) N e (cid:105) bulk (cid:104) N e ( δ ) (cid:105) = (cid:104) d T (cid:105) bulk (cid:104) d T ( δ ) (cid:105) , (6)where d T /d T,bulk is given by Eq. 5. As seen in Fig. 2, this theory was found to be a fairlygood model for the film-averaged changes in (cid:104) Z (cid:105) as a function of film thickness for systemsof N = 350 chains confined by hard, amorphous walls.In contrast, the theory systematically over-estimates the changes in entanglement for the N = 2000 free-standing films simulated here. On the one hand, this over-estimate is consis-tent with the expectation that free-standing films likely induce less chain-scale orientationalcorrelations than hard, reflecting walls (although it is not consistent with the expectationthat the film dimensions we consider might be expected to have larger changes to the distri-bution of chain end-to-end vectors [11] than are considered in the simplified conformationalmodel of Ref. [12]). On the other hand, more concerning is the apparent lack of collapsewhen the film sizes are reported in units of δ for our independent simulations of N = 500chains and subchain-derived data sets for chains of length N = 2000 , , , , , , , and 250. We emphasize, though, that for the range of film thicknesses studiedhere (and throughout much of the simulation literature) there is another crucial length scalecontributing to the problem: for such thin films the tube diameter itself may be comparableto the film thickness, and the ratio h eff /L e should also enter any attempted scaling collapseof the entanglement data.Although our data are quite limited, in Fig. 3a we see that the data for entanglement lossplausibly fall on a single surface when plotted simultaneously as a function of δ and h eff /L e ,a reflection that in very thin films both length scales are of importance in determining the11ystem-averaged reduction in the number of entanglements per chain. We propose an ansatzin which the dependence is separable: (cid:104) Z (cid:105)(cid:104) Z (cid:105) bulk = F ( δ, h eff /L e ) ≈ F ( δ ) F ( h eff /L e ) . (7)To test this, we approximate F ( h eff /L e ) by the value needed to collapse the data simulatedin this work as a function of δ . As shown in Fig. 3b, smoothly extrapolating this functionalso approximately collapses the data for capped films of N = 350 chains simulated in Ref.[12]. The values of F ( h eff /L e ) needed to collapse the data are shown in the inset of Fig.3b. Thus, we view the ansatz of Eq. 7 as plausible in light of the data available in thiswork, but caution that this collapse is over a relatively limited range of parameter space.This encourages a future systematic study of entanglement reduction as a function of δ and h eff /L e . We note that the experimental data from Ref. [13] are all at reasonably largevalues of h eff /L e , suggesting a possible explanation for the better data collapse observed inthat work. B. Spatial distribution of entanglement and conformational statistics
We now focus on the spatial inhomogeneities in the distribution of entanglement andconformational statistics throughout the films. This allows for tests of the spatially resolvedproperties of the entanglement network that (a) were predicted by the theory describedabove but (b) were not previously tested due to limits on the amount of data available [12].In this respect the much larger total number of particles simulated in this work is indis-pensible. Figure 4 shows the distribution of entanglement points (i.e. the location of thekinks determined by the Z1 algorithm) as a function of distance in the z -direction from thecenter of the film. Two points are of immediate interest. Qualitatively consistent with theexpectation that there are fewer chains to entangle with, near the film boundaries there isindeed a marked reduction in the entanglement density. Strikingly, though, there is not amonotonic increase in entanglement density towards the bulk; all three films plotted have amaximum in the entanglement density approximately 4 σ from the edge of the film (definingthe edge of the film by the fit value of ξ ). That distance is less than the characteristicdistance between kink points in the bulk, i.e., is less than the characteristic step size of thePP network. From the location of that maximum to the center of the film there is a decrease12 a) ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆◆◆◆◆◆▲ ▲▲▲▲▲▲▲▼▼▼▼▼▼▼▼ ○○○○○○ δ〈 Z 〉 / 〈 Z 〉 bulk F h eff L e ●■◆▲▼ h eff / L e F ( h e ff / L e ) (b) FIG. 3. (a) System-averaged reduction of entanglements as a function of the confinement scalerelative to both the chain dimension and L e . Data sets of constant h eff /L e (blue, orange, black,green, and purple points) correspond to the simulations of this work, and the remaining data (redpoints) correspond to the systems of N = 350 chains confined between hard amorphous wallsfrom Ref. [12]. (b) System-averaged reduction of entanglements normalized by an appoximatefunction of the dependence on h eff /L e , as described in the text. Inset. The value of the function F ( h eff /L e ) used to collapse the data sets in the main frame. Identical symbols in the main frameand inset correspond to the same data sets. in entanglement density that is qualitatively consistent with the theoretical predictions ofEq. 5 as applied to geometrically confined films [12]. Below we will explore a more quanti-tative comparison, but note that a qualitatively similar result is obtained by considering notthe density of entanglement points but the spatially resolved distance between successiveentanglement points along a chain, L e ( z ).To attempt to connect these measurements with the theoretical model described above,we have analyzed the distributions of angular orientations of both the chain end-to-endvector and the individual PP segments. Key assumptions of the theoretical model include(a) that the orientational correlations at the end-to-end vector scale created by geomet-13 | z - z cm | ρ k FIG. 4. Probability distribution of entanglements as a function distance in the z -direction fromplane corresponding to the center of mass of the film (solid curves) and the total density profileof the films (dashed curves). From bottom to top the curves correspond to N = 2000 chains infree-standing films with effective thickness h eff = 35 σ , h eff = 20 σ , and h eff = 13 . σ . - - < Cos ( θ )> P ( < C o s ( θ ) > ) FIG. 5. Film-averaged distributions of the angle with respect to the z -axis of the chain end-to-end vectors (dashed curves) and PP segments (solid curves). From bottom to top each set ofcurves corresponds to N = 2000 chains in free-standing films with effective thickness h eff = 35 σ , h eff = 20 σ , and h eff = 13 . σ . ric confinement are directly communicated to the PP scale, and (b) that these orientationalcorrelations at the PP scale alone allow one to predict changes in the properties of the entan-glement network. We will see that despite how well the theory predicts global entanglementproperties, neither of these assumptions is ultimately supported by the entanglement datagenerated by the Z1 algorithm.Figure 5 shows film-averaged data for the probability distributions of chain and PP14rientations with respect to the z axis. It is immediately clear that the orientations presenton the longest length scales are not communicated to the PP length scale, at least forthe primitive paths that are determined by the Z1 algorithm. This caveat is non-trivial,as the PP backbone detected by the Z1 algorithm may, in general, be different from theprimitive path steps detected by other algorithms. In particular, the way the chain shrinkingalgorithms affect the geometry of the melt may be crucial in the orientational statistics ofthe derived primitive paths, and it is entirely plausible that this change is in the directionof a more isotropic distribution of angles. Nevertheless, at least by this metric, the chainsare much more strongly orientationally ordered at the end-to-end scale than they are at theprimitive path scale. We note, though, that using the empirically observed distributionsof chain-scale orientational order leads to a more accurate prediction of the film-averagedentanglement density than using the simplified orientational distributions assumed by Ref.[12] and reported as the curve in Fig. 2: using the observed angular distributions in Eq:5 leads to predictions of (cid:104) Z (cid:105) / (cid:104) Z (cid:105) bulk = 0 . , . , and 0 .
76 for the three thinnest N =2000 films, for instance, compared with the measured values of (cid:104) Z (cid:105) / (cid:104) Z (cid:105) bulk = 0 . , . , and 0 .
77. In the context of the theory the success of the chain-scale predictions may beseen as a prediction of the “super-coarse-grained” needle mapping getting at some essentialphysics (albeit ignoring changes in the magnitude of the end-to-end vector of the chainsunder confinement, a severe approximation in the theory) [22], or it may be a reflectionthat the chain end-to-end orientation is a better representation of the entanglement-scaleorientational distribution than those reported by the Z1 algorithm.These angular distributions can be spatially resolved into angular distributions as a func-tion of the location of the end points of either the end-to-end vectors or the PP segments.This is shown in Fig. 6, where the data for each film is resolved into orientational distri-butions by layer in the film. The chain end-to-end vector orientational distributions, whilenoisy due to their more limited statistics, are slightly more sharply peaked around zero closeto the center of the film, as expected by chain-scale theories of conformations adopted inconfined geometries [11, 12]. In contrast, the PP-scale orientational distributions are muchmore strongly peaked at the film edges , and are essentially flat in the center of the film.While perhaps surprising that there is a qualitative disagreement here, one might expectthe chain-shrinking aspects of the Z1 algorithm to most strongly affect the angular distribu-tion of the PP network near the edges of the film. The length scale over which the primitive15 < Cos (θ)> PP P ( < C o s ( θ ) > PP ) (a) - < Cos (θ)> PP P ( < C o s ( θ ) > PP ) (b) - < Cos (θ)> PP P ( < C o s ( θ ) > PP ) (c) - - < Cos (θ)> R ee P ( < C o s ( θ ) > R ee ) (d) - - < Cos (θ)> R ee P ( < C o s ( θ ) > R ee ) (e) - - < Cos (θ)> R ee P ( < C o s ( θ ) > R ee ) (f) FIG. 6. Layer-resolved distributions of the angle with respect to the z -axis of PP segments (a-c)and the chain end-to-end vectors (d-f). Each set of three plots corresponds to data from N = 2000chains in free-standing films with effective thickness h eff = 35 σ , h eff = 20 σ , and h eff = 13 . σ .Within each graph, redder (more peaked) curves correspond to layers closer to the edge of the film,and bluer (flatter) curves correspond to layers near the film center of mass. paths adopt an essentially isotropic orientational correlations is ∼ σ . This is slightly longerthan the distance from the film edge of the entanglement density non-monotonicity shownin Fig. 4, and is shorter than L e .In principle the information in Fig. 6 is the direct ingredient in making spatially resolvedpredictions of the entanglement statistics with the theory in Eq. 5: one simply imaginestaking the integral over the angular distributions using the orientational distributions g ( (cid:126)u i )over whatever the local distribution is at that position in the film, rather than using theaverage over the entire ensemble. Figure 7 carries out this analysis, averaging over thelayer-resolved orientational distribution of both the primitive paths and the chain end-to-end vectors and comparing the normalized number of entanglement points in each layerwith the theoryetical expectation for the h eff = 35 σ and the h eff = 13 . σ film (whosebehaviors are representative of the films studied in this work). Other than the precipitousloss of entanglement near the film edge, the chain-scale theory predicts the h eff = 35 σ film properties very well. In contrast the chain-scale theory systematically over-predictsentanglement loss in the center of the thinner film and underpredicts it near the boundaries.The PP-scale theory is a qualitatively poor description of all films studied. The failure of16 -
10 0 10 200.00.20.40.60.81.0 z Z ( z ) Z bulk (a) - - Z ( z ) Z bulk (b) FIG. 7. Layer resolved, normalized distance between entanglement points (solid curve) for the N = 2000, h eff = 35 film (a) and the h eff = 13 . the chain-scale theoretical predictions are most acute in regimes where PP-scale informationis surely important, as indicated by the bulk properties discussed in the previous subseciton.The failure of the PP scale theory hinges on the connection between the geometry of thePP network revealed by the chain-shrinking algorithm and the true angular statistics of theentanglement strands in the melt. IV. DISCUSSION
In summary, we have studied both system-averaged and layer-resolved measures of en-tanglement statistics in a coarse-grained model of free-standing thin films made out of longlinear polymers. The simulations show that, in accord with films confined by hard, amor-phous walls, there is a pronounced decrease in the average number of entanglements per17hain with decreasing film thickness. The film-averaged values of the entanglement den-sity are reasonably well described by an existing theory [12], although this requires takingchain-scale orientational distributions rather than the PP-scale orientational distributionsthat are more natural in the context of the theory. Further, we find that existing proposalsthat collapse all entanglement results by plotting the entanglement density versus h eff /R ee are insufficient for the ultra-thin films simulated, and that an additional dependence on h eff /L e must be taken into account. While the separable-variable hypothesis in this paper, (cid:104) Z (cid:105) / (cid:104) Z (cid:105) bulk = F ( h eff /R ee ) F ( h eff /L e ), describes the data obtained in this paper fairlywell, formulating a theoretical description of the combined dependence of the entanglementdensity on both length scales is an outstanding challenge.While the film-averaged entanglement dilution is well-captured by the chain-scale ori-entational theory, the spatially resolved nature of this entanglement dilution is less welldescribed. The non-monotonicity of the entanglement density versus distance from the cen-ter of the film observed in Fig. 4 is a notable, qualitative prediction of the theory [12],but the quantitative details of the spatially-resolved nature of the entanglement dilution aremore poorly captured in the theory. We then detailed the orientational statistics observedat both the chain- and PP-scale, noting that two of the key assumptions of the theory do notseem to be supported by the simulation data. We note that in these thin-film systems thereis a weak segregation of chain ends near the interfacial layers of the film. While the chainends themselves do not correspond to entanglements, this nevertheless leads to a slightlyhigher number of entanglements at a characteristic distance (the combination of the typicalspacing between entanglements and the typical orientation of PP segments near the inter-face) from the free surface. However, this is a comparatively weak effect, affectly too smalla percentage of the total entanglements in the system to account for the effects seen in Figs.4 and 7.A key observation of this work, and a strong candidate for the failure of the theory tocorrectly predict the spatially resolved properties studied in this work, is that both thefilm-averaged and the layer-resolved PP orientational distribution are very different fromorientational distributions of the chain end-to-end vectors. It is qualitatively plausible thatin Fig. 7B the correct theoretical curve lies in between the one derived from the chainend-to-end vector orientational distribution and the one derived from the primitive pathinformation, in particular in the regime where the film thickness is not large compared to18 e . In polymer systems one can define a length-scale dependent orientational distributions,and this observation raises the possibility that it is not the orientational measures at theentanglement scale but some other length that dominates the problem.A strong caveat to this, though, is that we have taken the primitive paths to be deter-mined via the Z1 algorithm, and in particular the geometric properties of the primitive pathsare likely to be strongly affected by any of the algorithms that use chain-shrinking methodsto uncover the entanglement network. It is not surprising that the chain-shrinking algo-rithms tend to make the angular distributions of the primitive paths appear more isotropicthan the entanglement strands in the unperturbed confined melt, but in this work we seethat the algorithm may be completely erasing all informations about the geometry on thisscale. Whether alternate approaches to finding the entanglement strands, for instance analternate chain shrinking algorithm such as CReTA [39], or isoconifgurational time averag-ing of the chains [42], lead to strongly different geometric and orientational properties atthe entanglement strand level will be an object of future study.In general these results raise questions about the applicability of the theory, particularlywhen the primitive path orientational distribution does not match the chain-scale orienta-tional distribution. We propose that new simulations of geometrically confined polymers,both free-standing and with hard walls, be conducted in order to further clarify these ques-tions. We believe that simulations of reasonably long chains are crucial here, so that a cleandistinction can be made between lengths on the order of the tube diameter and lengths ofthe order of the chain radius of gyration. Of particular interest would be a more completesurface describing the dependence on film-averaged entanglement statistics on both h eff /R ee and h eff /L e . We also propose that new simulations and analyses be performed on entan-gled polymers under deformation; the entanglement dilution effect observed in steady-statecontinuous shear simulations [24] is very well captured by the theory [43], so in light ofthe present results a more critical evaluation of the theory, which would begin at the PPorientational scale, is called for. In closing, we reiterate that the connection between chainstructure away from bulk, equilibrium conformations and either the statistics of the entan-glement network or the dynamical properties of the system remains poorly understood (andnot captured by either the classic Doi-Edwards tube model nor any recent extensions ofit). We hope that the present results spurs the development of new theoretical descriptionselucidating these relationships. 19 CKNOWLEDGMENTS
This work was supported by the Advanced Materials Fellowship of the American Philo-sophical Society. The Tesla K40 used for this research was donated by the NVIDIA Cor-poration. The author would like to thank Rob Riggleman and Ken Schweizer for fruitfuldiscussions and comments on this manuscript. [1] J. L. Keddie, R. A. L. Jones, and R. A. Cory, Europhys. Lett. , 59 (1994).[2] P. A. O’Connell and G. B. McKenna, Science , 1760 (2005).[3] P. A. O’Connell and G. B. McKenna, Eur. Phys. J. E , 143 (2006).[4] J.-L. Masson and P. F. Green, Phys. Rev. E , 031806 (2002).[5] H. D. Rowland, W. P. King, J. B. Pethica, and G. L. W. Cross, Science , 720 (2008).[6] K. Shin, S. Obukhov, J.-T. Chen, J. Huh, Y. Hwang, S. Mok, P. Dobriyal, P. Thiyagarajan,and T. P. Russel, Nature Materials , 961 (2007).[7] F. Lange, P. Judeinstein, C. Franz, B. Hartmann-Azana, S. Ok, M. Steinhart, andK. Saalw¨achter, ACS Macro Lett. , 561 (2015).[8] A. Karantrantos, R. J. Composto, K. I. Winey, M. Kr¨oger, and N. Clarke, Macromolecules , 7274 (2012).[9] Y. Li, M. Kr¨oger, and W. K. Liu, Phys. Rev. Lett. , 118001 (2012).[10] M. P. Weir, D. W. Johnson, S. C. Boothroyd, R. C. Savage, R. L. Thompson, S. M. King,S. E. Rogers, K. S. Coleman, and N. Clarke, ACS Macro Lett. , 430 (2016).[11] A. Cavallo, M. M¨uller, J. P. Wittmer, A. Johner, and K. Binder, J. Phys.: Condens. Matter , S1697 (2005).[12] D. M. Sussman, W.-S. Tung, K. I. Winey, K. S. Schweizer, and R. A. Riggleman, Macro-molecules , 6462 (2014).[13] L. Si, M. V. Massa, K. Dalnoki-Veress, H. R. Brown, and R. A. L. Jones, Phys. Rev. Lett. , 127801 (2005).[14] Y. Liu, Y.-C. Chen, S. Hutchens, J. Lawrence, T. Emrick, and A. J. Crosby, Macromolecules , 6534 (2015).[15] A. Cavallo, M. M¨uller, and K. Binder, J. Phys. Chem. B , 6544 (2005).
16] H. Meyer, T. Kreer, A. Cavallo, J. P. Wittmer, and J. Baschnagel, J. Eur. Phys. J. – Spec.Top. , 167 (2007).[17] M. Vladkov and J.-L. Barrat, Macromolecules , 3797 (2007).[18] V. Sethuraman, D. Kipp, and V. Ganesan, Macromolecules , 6321 (2015).[19] M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press,Oxford, UK, 1986).[20] W.-S. Tung, R. J. Composto, R. A. Riggleman, and K. I. Winey, Macromolecules , 2324(2015).[21] H. R. Brown and T. P. Russell, Macromolecules , 798 (1996).[22] D. M. Sussman and K. S. Schweizer, Phys. Rev. Lett. , 168306 (2012).[23] D. M. Sussman and K. S. Schweizer, J. Chem. Phys. , 234904 (2013).[24] C. Baig, V. G. Mavrantzas, and M. Kr¨oger, Macromolecules , 6886 (2010).[25] M. H. Nafar Sefiddashti, B. J. Edwards, and B. Khomami, J. Rheol. , 1 (2014).[26] K. Kremer and G. S. Grest, J. Chem. Phys. , 5057 (1990).[27] “Hoomd-blue,” https://codeblue.umich.edu/hoomd-blue/ .[28] J. A. Anderson, C. D. Lorenz, and A. Travesset, J. Comput. Phys. , 5342 (2008).[29] R. Auhl, R. Everaers, G. S. Grest, K. Kremer, and S. J. Plimpton, J. Chem. Phys. ,12718 (2003).[30] “Bitbucket repository,” https://bitbucket.org/dmsussman/ .[31] F. P. Buff, A. Lovett, and F. H. Stillinger Jr., Phys. Rev. Lett. , 621 (1965).[32] T. Hapke, G. P¨atzold, and D. W. Heermann, J. Chem. Phys , 10075 (1998).[33] N. C. Karayiannis, V. G. Mavrantzas, and D. N. Theodorou, Phys. Rev. Let.. , 105503(2002).[34] R. S. Hoy and M. O. Robbins, Phys. Rev. E , 061802 (2005).[35] M. Kr¨oger, Comput. Phys. Commun. , 209 (2005).[36] N. C. Karayiannis and M. Kr¨oger, Int. J. Mol. Sci. , 5054 (2009).[37] R. S. Hoy, K. Foteinopoulou, and M. Kr¨oger, Phys. Rev. E , 031803 (2009).[38] S. Shanbhag and M. Kr¨oger, Macromolecules , 2897 (2007).[39] C. Tzoumanekas and D. N. Theodorou, Macromolecules , 2456 (2006).[40] G. Szamel, Phys. Rev. Lett. , 3744 (1993).[41] D. M. Sussman and K. S. Schweizer, Phys. Rev. Lett. , 078102 (2011).
42] W. Bisbee, J. Qin, and S. T. Milner, Macromolecules , 8972 (2011).[43] K. S. Schweizer and D. M. Sussman, arXiv:1604.03972 (2016)., 8972 (2011).[43] K. S. Schweizer and D. M. Sussman, arXiv:1604.03972 (2016).