Efficient equilibration of confined and free-standing films of highly entangled polymer melts

Equilibration of polymer melts containing highly entangled long polymer chains in confinement or with free surfaces is a challenge for computer simulations. We approach this problem by first studying polymer melts based on the soft-sphere coarse-grained model confined between two walls with periodic boundary conditions in two directions parallel to the walls. Then we apply backmapping to reinsert the microscopic details of the underlying bead-spring model. Tuning the strength of the wall potential, the monomer density of confined polymer melts in equilibrium is kept at the bulk density even near the walls. In a weak confining regime, we observe the same conformational properties of chains as in the bulk melt showing that our confined polymer melts have reached their equilibrated state. Our methodology provides an efficient way of equilibrating large polymer films with different thicknesses and is not confined to a specific underlying microscopic model. Switching off the wall potential in the direction perpendicular to the walls, enables to study free-standing highly entangled polymer films or polymer films with one supporting substrate.


I. INTRODUCTION
Polymer confinement plays an important role for many aspects of adhesion, wetting, lubrication, and friction of complex fluids from both theoretical and technological points of view. For more than two decades, both theoretical and experimental works [1][2][3][4][5][6][7] have shown that the dynamic and structural properties of polymers subject to confinement may deviate from that in the bulk as the interaction between polymers and confining surfaces becomes non-negligible. For example, the conformations of polymer chains in a melt near the wall in the direction perpendicular to the wall shrink remarkably compared to bulk chains while they only extend slightly parallel to the wall 3,8,9 . Long chain mobility in confined melts is also affected by entanglement effects while in turn the conformational deviations from bulk chains have influence on the distribution of entanglements [10][11][12][13][14] . Furthermore many studies have also focused on the dependency between the glass transition temperature T g and the nature of the confinement effects [15][16][17] . Therefore, it is important to understand the mechanical properties of confined polymer melts and how confinement impacts both viscous and elastic properties of amorphous polymer films with different surface substrates or even free surfaces.
Computer simulations provide a powerful method to mimic the behavior of polymers under well-defined external conditions covering the range from atomic to coarse-grained (CG) scales 18,19 . However, the cost of computing time rises dramatically as the size and complexity of systems increase. Therefore, applying appropriate coarse-grained models that keep the global thermodynamic properties and the local mechanical and chemical properties is still an important subject [20][21][22][23][24][25][26][27][28] . One of the successful monomer-based models, namely the bead-spring (BS) model 21,22 together with an additional bond-bending potential [29][30][31][32] , has been successfully employed to provide a better understanding of generic scaling properties of polymer melts in bulk. For such a model static and dynamic properties of highly entangled polymer melts in bulk have been extensively studied in our previous work 33,34 . We have verified the crossover scaling behavior of the mean square displacement of monomers between several characteristic time scales as predicted by the Rouse model 35 and reptation theory [36][37][38][39] over several orders of magnitude in the time. For weakly semiflexible polymer chains of sizes of up to N = 2000 monomers we have also confirmed that they behave as ideal chains to a very good approximation. For these chains the entanglement length N e = 28 in the melt, estimated through the primitive path analysis and confirmed by the plateau modulus 32,33,40 , resulting in polymer chains of N = 2000 ≈ 72N e . Thus, we here focus on this model and a related variant for the study of the equilibration of polymer melts confined between two repulsive walls as supporting films, and free-standing films after walls are removed. Each film contains n c = 1000 chains of N = 2000 monomers at the bulk melt density ρ 0 = 0.85σ −3 . Directly equilibrating such large and highly entangled chains in bulk or confinement is not feasible within reasonably accessible computing time.
A novel and very efficient methodology has recently been developed 27,41 for equilibrating large and highly entangled polymer melts in bulk described by the bead-spring model 21,22 . Through a hierarchical backmapping of CG chains described by the soft-sphere CG model 26,27 from low resolution to high resolution and a reinserting of microscopic details of bead-spring chains, finally, highly entangled polymer melts in bulk are equilibrated by molecular dynamics (MD) simulations using the package ESPResSO++ 42,43 . To first order, the required computing time depends only on the overall system size and becomes independent of chain length. Similar methodologies have also been used to equilibrate high-molecular-weight polymer blends 44 and polystyrene melts 45 . In this paper, we extend the application of the soft-sphere approach to confined polymer melts and subsequently free-standing films. As polymer chains are described at a lower resolution, the number of degrees of freedom becomes smaller. Here we adapt this hierarchical approach to equilibrate polymer melts confined between two walls in detail. Moreover, we apply our newly developed, related model 46,47 to prepare polymer films with one or two free surfaces. Differently from Refs. 26 and 41, we take the bending elasticity of bead-spring chains in a bulk melt into account for the parameterization of the soft-sphere CG model. Namely, the underlying microscopic bead-spring chains are weakly semiflexible (bending constant k θ = 1.5ε) instead of fully flexible (k θ = 0.0ε).
The outline of this paper is as follows: In Sec. II, we introduce the main features of the microscopic bead-spring model and soft-sphere coarse-grained model used for studying the confined polymer melts. The application of the soft-sphere CG model for confined coarse-graining melts, and conformational properties of fully equilibrated confined CG melts are addressed in Sec. III.
In Sec. IV, we reinsert the microscopic details of confined CG melts, and discuss the equilibration procedures. In Sec. V, we show how to prepare films with one or two free surfaces by switching to another variant of bead-spring model. Finally, our conclusion is given in Sec. VI.

A. Generic microscopic bead-spring models
In the microscopic bead-spring model 21,22 , all monomers at a distance r interact via a shifted Lennard-Jones (LJ) potential U LJ (r), with where ε is the energy strength of the pairwise interaction, r cut is the cutoff in the minimum of the potential such that force and potential are zero at r cut = 2 1/6 σ . These LJ units also provide a natural time definition via τ = σ m/ε, m = 1 being the mass of the monomers. The temperature is set to T = 1.0ε/k B , k B being the Boltzmann factor, which is set to one. Along the backbone of the chains a finitely extensible nonlinear elastic (FENE) binding potential U FENE (r) is added, where k = 30ε/σ 2 is the force constant and R 0 = 1.5σ is the maximum value of bond length. For controlling the bending elasticity, i.e., chain stiffness, the standard bond-bending potential 29-32 is given by with the bond angle θ defined by where b j = r j+1 − r j is the bond vector between the ( j +1)th monomer and the jth monomer along the identical chain. The bending factor k θ is set to 1.5ε and the melt density is set to the widely used value of 0.85σ −3 throughout the whole paper.
For studying polymer melts under confinement, we first consider the simpler example of polymer melts that are confined between two planar, structureless repulsive walls. The walls placed at z = 0 and z = L z are described by the 10-4 Lennard-Jones planar wall potential 3,48 , with Here ε w is the interaction strength between monomers and the walls, z and L z − z are the vertical distances of a monomer from the two walls, respectively.
For preparing polymer films with one or two free surfaces, we have to stabilize the system at zero pressure, particularly, in the direction perpendicular to the walls. Only then we can switch off the wall potential and prevent system instability. Therefore, a short-range attractive potential to reduce the pressure to zero is added 46,47 with an additional shift term, such that U ATT (r) = U LJ (r) at r = r cut . Here α = 0.5145ε is the strength parameter, r a c = √ 2r cut ≈ 1.5874σ is the upper cut-off such that it has zero force at r = r cut and r = r a c . Note that this additional potential does not alter the characteristic conformations at T = 1ε/k B , so that we can switch between these different models as needed 21,22,46,47,49 .
Note that using U (old) BEND (θ ), bead-spring chains tend to stretch out as the temperature decreases. To avoid such an artificial chain stretching, Eq. (4) can be replaced by 46,47 if one is interested in studying polymer melts under cooling. The fitting parameters a θ = 4.5ε and b θ = 1.5 are determined so that the local conformations of chains remain essentially unchanged compared to those with k θ = 1.5ε using Eq. (4) at temperature T = 1.0ε/k B .

B. Soft-sphere coarse-grained model
In the soft-sphere approach 26 , each bead-spring polymer chain in a melt is represented by a chain of N CG = N/N b fluctuating soft spheres in a CG view as shown in Fig. 1. The coordinate of the center and the radius of the sth sphere in the ith chain, r i (s) and σ i (s), are described by and respectively. Here r i,k is the coordinate of the kth monomer in chain i, R CM,i (s) and R g,i (s) are the center of mass (CM) and the radius of gyration of N b monomers in its subchains, respectively.
Since we extend the application of the soft-sphere CG model to semiflexible chains and follow the slightly different strategies of Refs. 27 and 41, we list the details of soft-sphere potentials as follows. The spheres are connected by a harmonic bond potential and an angular potential where d i (s) = r i (s + 1) − r i (s) is the bond vector, and b CG denotes the effective bond length.
The radius of the sth sphere in chain i and its fluctuation are controlled by the following two potentials 26,27 , and respectively. The non-bonded interactions between any two different spheres due to excluded volume interactions according to Flory's theory 20,36,50 is taken care of by the potential Here the parameters k bend , a 1 , a 2 , a 3 , b CG , and ε 1 depending on N b are determined by a numerical approach via curve fitting.
Similar as shown in Eq. (5), the two soft repulsive planar walls located at z = 0 and z = L z depending on the radius of each sphere is given by with V (CG) where ε The parameters  Since the excluded volume effect between N b monomers in each subchain is ignored in the parameterization of the soft-sphere approach and self-entanglements on this length scale is negligible (N b = 25 < N e = 28) subchains behave as ideal chains (alternative, one could include excluded volume for semi dilute solutions via a Flory term). In this case, we can simplify several steps of hierarchical backmapping 41 to only one step of fine-graining to introduce microscopic details of subchains once a CG melt reaches its equilibrated state.

III. EQUILIBRATION OF SOFT-SPHERE CHAINS IN A CONFINED CG MELT
In this section, we extend the application of the soft-sphere CG model for polymer melts in bulk to polymer melts confined between two repulsive walls {Eq. (16)}, first focusing on the CG melt containing n c = 1000 chains of N CG = 80 spheres at the CG bulk melt density ρ For the comparison to our reference systems, we set the distance between two walls compatible with the bulk melt. Thus, we locate two walls at z = 0σ , and z = L z ≈ 133σ while keeping the periodic boundary conditions along the xand ydirections with the lateral linear dimensions L x = L y = 133σ . Of course, one can adjust L z and extend/reduce L x , L y for keeping the bulk melt density as needed.
The initial configurations of soft-sphere chains in terms of {σ i (s), d i (s), θ i (s)} are randomly generated according to their corresponding Boltzmann weights exp[−βU is computationally more efficient to perform Monte Carlo simulations to equilibrate confined CG melts. Similar to Ref. 27, our simulation algorithm including three types of MC moves at each step is as follows: (i) For a local move, one of n c N CG spheres is randomly selected, e.g.
the sth sphere in the ith chain, the sphere of radius where η is a random number and η ∈ [0, 1). (ii) For a snake-slithering move, one end of n c chains is randomly selected, and σ  26 since the contributions for r > 20σ are negligible. Nevertheless, there is no influence on measurements of any physical observable while it speeds up the simulations by a factor of four.
Applying a linked-cell algorithm with the cell size L c = 2.66σ (L x /L c = L y /L c is very close to an integer), smaller than the cut-off value r (CG) cut , speeds up the simulation even more by an additional factor of 2.5, i.e. all together it speeds up by a factor of ten. It takes about 10 hours CPU time on an Intel 3.60GHz PC for a confined CG melt to reach its equilibrated state (after 2 × 10 7 MC steps are performed). The acceptance ratio is about 73% for a trial change of sphere size, 45% for a local move, and 41% for a snake-slithering move.
Choosing the wall strength ε , and P(θ ) comparing to that for an equilibrated CG melt in bulk as shown in Fig. 2. Since soft spheres are allowed to penetrate each other in CG melts, we observe that the distributions P(σ , N b ) and P(d, N b ) are slightly narrower for both equilibrated CG melts in bulk and in confinement than that for the reference systems while the average values of radius and bond length remain the same. We have also compared the estimates of R 2 (s, N b ) and g(r, N b ) for an equilibrated confined CG melt to an equilibrated CG melt in bulk and the reference data in Fig. 3. The curve of R 2 (s, N b ) taken the average over all n c = 1000 chains for a confined CG melt deviates from the bulk behavior for s > 10 due to the confinement effect while for s < 6, is a bit smaller compared to the bulk value. It is due to the artifact of the soft-sphere CG model, where the excluded volume effect within the size of spheres is not considered. After monomers are reinserted into soft-sphere chains, local excluded volume and the corresponding correlation hole effect automatically correct for these deviations. However, the discrepancy for s < 6 is still within fluctuations observed in bulk. When monomers are reinserted into soft-sphere chains, the estimate of g(r, N b ) starts to increase at r ≈ 3σ , and then decrease at r ≈ 5σ . It indicates that near the walls, the distance between any two spheres decreases due to the confinement effect.
To investigate the confinement effect on packing and conformations of a polymer melt, we determine the soft-sphere density profile between two walls as follows, since P(σ , N b ) has its maximum at z = σ i (s) ≈ 3.06σ (see Fig. 2). The two components of the mean square radius of gyration depending on the z-component of CMs of chains are defined as follows, where r i,α=|| (s) = r i,x (s) + r i,y (s) and r i,α=⊥ (s) = r i,z (s).  (20) and where k cm and k g determine the coupling strength. The forces derived from these two potentials can drive the center of mass and the radius of gyration of each subchain to the center and the radius of the corresponding soft sphere, respectively. Namely, each bead-spring chain is then sitting on top of its corresponding soft-sphere CG chain (see Fig. 1 in a microscopic representation, we set L z = 134σ instead of L z = 133σ taking the repulsive potential of the wall, which is a steep but smooth function, into account. The reinsertion MD time is about 80τ and the CPU time is about 2.8 hours on a single processor in an Intel 3.60GHz PC.

B. Equilibration procedure
In the next step the full excluded volume interaction as listed in Section II A has to be introduced. To avoid the "explosion" of the system due to the overlap of monomers, we have to switch on the excluded volume interactions between the non-bonded pairs of monomers in a quasi-static way (slow push-off) 40,49 . Therefore, the shifted LJ potential for each non-bonded pair of monomers at a distance r, {Eq. (1)}, is first replaced by a force-capped LJ potential where r fc is an adjustable cut-off distance in this warm-up procedure. r fc = r o fc decreases monotonically at each cycle from 2 1/6 σ to 0.8σ in the warm-up procedure for the non-bonded pairs except the next-nearest neighboring (nn) pairs along identical chains. For the nn pairs, r fc = r nn fc is set to 1.1σ initially, and tuned according to the following cost function 40 (23) at the end of each cycle in the warm-up procedure with the restriction that 0.8σ ≤ r nn fc ≤ 1.1σ . Master curve refers to fully equilibrated polymer melts in bulk (the reference systems). In this process the nn-excluded volume is adjusted according to the value of C, namely reduced (increase creasing the cut-off value of r o fc for non-bonded monomer pairs except nn pairs are required for the confined polymer melt system due to the competition between the excluded volume effect and the confinement effect. In the first 100 cycles, r nn fc is updated at each cycle associated with the cost function. In the rest 20 cycles, r nn fc decreases by 0.01σ before reaching the minimum value 0.8σ (see Fig. 6a). (b) In the relaxation procedure, the shifted LJ potential U LJ (r) is restored. We first perform 10 5 MD steps with Γ = 0.5τ −1 and ∆t = 0.001τ, and then another 2 × 10 6 MD steps with ∆t = 0.005τ to ensure that the confined melt reaches its equilibrated state. Afterwards, we can set the time step to its standard value, ∆t = 0.01τ for the further study of the confined polymer melt in equilibrium.
During warming up the confined polymer melt, the local monomer density profile ρ(z) near the walls varies with the interaction strength ε w between monomers and walls as shown in Fig. 6a.
For ε w = 0.005ε, the bulk melt density ρ 0 is conserved all the way up to the wall, where the bin size is set to 0.5σ . At the same time, after the confined polymer melt is warmed up, the curve of R 2 (s) coincides with the master curve for s < 50 as shown in Fig. 6b. A typical variation of cut-off distances, r nn fc and r o fc for the non-bonded pairs of monomers in the warm-up procedure are shown in Fig. 6c. At the end of the warm-up procedure, U FC−LJ (r) approaches to U LJ (r). Thus, there is no problem to switch back to U LJ (r) for the further relaxation of the confined polymer melt.
First we examine the difference between the equilibrated confined melt and the reference bulk melt as shown in Fig. 7. We compare the whole system in terms of the mean square internal distance R 2 (s) , the single chain structure factor S c (q), and the collective structure S || (q) in the direction parallel to the walls. We see that the curve of R 2 (s) estimated only for those chains having their CMs satisfying 0.3L z ≤ r In (b), the master curve for the average behavior of the reference systems is also shown for comparison.
any anisotropy of chain conformations under confinement, we distinguish S c,⊥ (q = q z ) where the wave vector q is oriented in the z-direction perpendicular to the wall, and S c,|| (q = (q 2 x + q 2 y ) 1/2 ). S c,|| (q) follows the same chain structure as in the bulk melt as shown in Fig. 7b. The shift of S c,⊥ (q) toward to a slightly larger value of q from the bulk curve also indicates that the estimate of R 2 g,⊥ for all 1000 chains is smaller. Nevertheless, chains still behave as ideal chains. Comparison of the collective structure factor of the whole melt between the confined polymer and the bulk melt in the parallel direction to the walls, we see that there is no difference as the distance between two walls is compatible to the bulk melt (see Fig. 7c). The local packing of monomers characterized by the pair distribution g(r) for both inter and intra pairs of monomers is also in perfect agreement with the bulk melt as shown in Fig. 7d.
Finally, we go back to the outset and analyze the conformational properties of the confined equilibrated melt in a CG representation by mapping bead-spring chains to soft-sphere chains using Eqs. (9), (10). As shown in Fig. 3 and Fig. 4b keeping the short range repulsion from the walls. Fig. 8a shows that the three diagonal terms of the pressure tensor P αβ (t) drop from 4.9ε/σ 3 (5.0ε/σ 3 for bulk melts) to 0.1ε/σ 3 in a very short time about 20τ. We further relax the confined polymer film for 30000τ ≈ 13τ e , τ e ≈ 2266τ being the entanglement time 33 , by performing MD simulations in the NPT ensemble (Hoover barostat with Langevin thermostat 52,53 implemented in ESPResSo++ 42,43 ) at temperature T = 1.0ε/k B and pressure P = (P xx +P yy +P zz )/3 = 0.0ε/σ 3 to finally adjust the pressure from 0.1ε/σ 3 to 0.0ε/σ 3 .
Under this circumstance, an equilibrated free-standing film is generated after removing two walls by turning off the wall potential at z = 0σ and z = L z = 134σ . If we only remove one of the walls, we get a polymer film with one supporting substrate, where one of course can introduce appropriate adhesion interactions.
The overall conformations of all chains and inner chains (0.3L z ≤ r (i) CM,z ≤ 0.7L z ) as characterized by R 2 (s) between two walls for a confined polymer melt based on two variants of beadspring model, and after relaxing for 30000τ (NPT MD simulations) are preserved within fluctuation as shown in Fig. 8b,c. The monomer density profile ρ(z) in the direction perpendicular to the wall compared to the density in the interior of confined melt is also preserved as shown in Fig. 8d.
The thickness of films, D, is determined according to the concept of Gibbs dividing surface that has been applied to identify the interface between two different phases 54-56 , e.g. liquid and vapor, polymer and vacuum, etc. based on the density profile in the direction perpendicular to the inter-    Following the same strategy, one can easily equilibrate highly entangled polymer melts confined between two walls at distances ranging from thick films to thin films (in which the distance between two walls is smaller than the radius of gyration of chains) within easily manageable computing time. Our work opens ample possibilities to study static and dynamic properties of highly entangled polymer chains in large polymer films, including e.g. entanglement distributions. Varying the interaction potential between walls and monomers, or even replacing the wall potential by other potentials, only requires local short relaxation runs starting from a fully equilibrated polymer melt confined between two repulsive walls. Switching to our recently developed coarse-grained model for studying polymer melts under cooling 46,47 , both fully equilibrated confined and free-standing films at the temperature T = 1.0ε and pressure P = 0.0ε/σ 3 are also obtained in this work. This provides a direct route to further investigate the relation between the glass transition temperature and the thickness of films of highly entangled polymer chains at the zero pressure. Beyond that, it is also interesting to analyze the rheological properties and local morphology of deformed films by stretching or shearing.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.