Simulations of Thin Polymer Films on Flat and Curved Substrates

Properties of the bulk systems are denoted by the subscript $0$ and are shown in Table S1. For the flat substrates, we characterize the local polymer film properties depending on the distance $z$ perpendicular to the substrate, such as the monomer density $\rho(z)$, the polymer center-of-mass (CM) density $\rho_{\text{CM}}(z)$, and the mean square radius of gyration of the chains $\langle R_{g}^{2}\rangle(z)$ ^{48, 49}. For the curved substrates, we work both with Cartesian 2D maps in the $xz$-plane, and with plots in curvilinear coordinates as described in Fig. 2. Here, we explore offset curves ⁵⁰ as a function of arc length and conversely analyze fixed positions on the curve as function of normal distance. Each bin is a rectangular cuboid of sides $\Delta n=0.2\,\sigma$, $\Delta s=0.5\,\sigma$ and $L_{y}$.

For reporting chain properties, we use the center-of-mass position of the chain to perform the spatial binning. Due to the nature of the curvilinear coordinates, normal lines on positive (negative) curvature regions will converge (diverge). Thus, a monomer/polymer may contribute to more than one bin in the arc length direction. We averaged the results over $500$ snapshots, separated by $1000\,\tau$, and used the mirror symmetry of the substrate to improve statistics.

The stability of a polymer film depends crucially on the wetting conditions between the monomers and the substrate. The basic equation that governs capillary phenomena at interfaces is Young’s equations, which relates the interfacial tensions $\gamma_{ij}$ between the substrate (s), liquid (l) and vapor (v) phases to the contact angle $\theta$ that the liquid-vapor interface forms with the solid substrate

where $S=\gamma_{\text{sv}}-(\gamma_{\text{sl}}+\gamma_{\text{lv}})$ is the so-called spreading coefficient. For a droplet that fully wets the substrate, the contact angle at equilibrium is $\theta=0^{\circ}$ such that $S\geq 0$. To form stable thin films, the adhesive substrate-monomer forces must be strong enough to outweigh the cohesive forces in the polymer film. Otherwise, the polymers will form a droplet with a contact angle of $\theta>0^{\circ}$ that might even desorb completely ($\theta=180^{\circ}$).

In our simulation model, the adhesion strength is directly controlled through the parameter $\varepsilon_{\text{w}}$ of the LJ potential acting between monomers and substrate beads, while the strength of the cohesive forces within the film are tuned via the monomer-monomer interaction strength $\varepsilon$. In the following, we will present a simple continuum model to estimate the contact angle of a (polymer) droplet on a flat substrate as a function of $\varepsilon_{\rm w}$. We assume that the polymers form a spherical cap with radius $R$ and contact angle $\theta$ [Fig. 3(a)]. We approximate the potential energy of the system (neglecting contributions from bonded interactions, which are identical in all cases) as:

where $V_{\text{p}}$ and $A_{\text{p}}$ are the volume and surface area of the spherical cap, and $A_{\rm srf}=\pi c^{2}$ is the surface area of the disk in contact with the substrate. Note that $A_{\text{p}}$ includes the area of the droplet-substrate interface, $A_{\rm srf}$ [indicated by blue region in Fig. 3(a)]. Further, we assume a homogeneous monomer number density within the droplet with value identical to the bulk, i.e., $\rho_{0}\approx 0.91\,\sigma^{-3}$ for our model. From the bulk simulations, we further determined the potential energy per particle due to the LJ potential as $e_{0}\equiv E_{0}/N\approx-3.6\,\varepsilon$ (Table S1). In deriving Eq. (6), we have also assumed that the monomers in the outermost droplet layer with thickness $h=\sigma$ have half the number of neighbors, and that the potential energy due to the wall-monomer interaction can be described by a step-like potential with (effective) coordination number $Z=16.88$, that characterizes the number of contacts between an adsorbed monomer and substrate beads. Since the polymer droplet is incompressible, the droplet volume $V_{\text{p}}$ must be preserved making $R\equiv R(\theta)$. With this constraint, the potential energy written in Eq. (6) becomes a function of the contact angle $\theta$. We then identify the preferred contact angle $\theta$ of our model by finding the minimum of Eq. (6).

Figure 3(c) shows that the contact angle $\theta\to 180^{\circ}$ in the limit of $\varepsilon_{\text{w}}\to 0$, as expected. With increasing wall attraction $\varepsilon_{\text{w}}$, $\theta$ monotonically decreases until it saturates at $\theta=0^{\circ}$ for $\varepsilon_{\text{w}}\gtrsim 0.2\,\varepsilon$. To test these estimates of the contact angle based on simple energetic considerations, we performed MD simulations in which we initially placed a hemispherical polymer droplet of radius $R=30\,\sigma$ (consisting of $N=51456$ monomers) on a flat substrate. We then periodically increased the wall attraction in steps of $\Delta\varepsilon_{\text{w}}=0.05\,\varepsilon$ and let the system evolve for $1000\,\tau$. We verified that this time is sufficient to let the polymer droplet relax into its equilibrium shape, and determined the contact angle $\theta$ of the droplet with the planar substrate [Fig.3(b)]. The resulting contact angles $\theta$ from the MD simulations and our continuum model are in excellent agreement, which provides an a posteriori justification for the simplifications of our model calculation.

At this point, it should be noted that our continuum model does not consider chain connectivity effects explicitly, but only implicitly through the bulk monomer number density $\rho_{0}$, which depends on the degree of polymerization $M$.⁵¹ To test the effect of this approximation, we performed additional MD simulations of shorter dimer chains and determined the resulting contact angles $\theta$. As can be seen in Fig. 3(c), there are only minor differences in $\theta$ for the two different chain lengths: At a given $\varepsilon_{\text{w}}$, the contact angle $\theta$ was systematically lower for the dimers, and the transition to the fully wetted regime occurred at a weaker wall attraction. These findings are in agreement with recent MD simulations and density functional theory calculations by Midya et al.²⁹ who observed that the temperature for complete wetting ($\theta=0^{\circ}$) increased with increasing chain length $M$. The differences between the short and long polymers fully disappear for attraction strengths well above the wetting transition, $\varepsilon_{\text{w}}\gg 0.2\,\varepsilon$. In practice, however, one can not use an arbitrarily large interaction strength $\varepsilon_{\text{w}}$, because a too strong substrate-monomer attraction will result in an increase of the local monomer density near the substrate, eventually crossing the phase boundary to a crystalline solid. Indeed, for $\varepsilon_{\text{w}}\gtrsim\varepsilon$ we observed that the monomers near the substrate formed a thin crystalline layer with hexagonal in-plane order, which greatly reduced the mobility of the adsorbed monomers (see SI).

To establish a reference for polymer films on curved substrates, we first simulated a thin polymer on a planar substrate (see SI for details on the film preparation). We set the substrate-monomer interaction strength to $\varepsilon_{\text{w}}=0.35\,\varepsilon$, so that the polymers fully wet the substrate. The nominal film thickness is set to $H\approx 14.8\,\sigma\approx 6.2\,R_{\text{g},0}$, with $R_{\text{g},0}\equiv\left\langle R_{\text{g},0}^{2}\right\rangle^{1/2}\approx 2.4\,\sigma$ being the root mean square radius of gyration of the chains in the bulk melt. Figure 4(a) shows the normalized monomer number density in the direction perpendicular to the substrate, $\rho(z)/\rho_{0}$. Near the substrate, the density profiles show a distinct layering, which is typical for short molecules near hard walls,⁵² with the first peak occurring near the minimum of the integrated substrate-monomer interaction potential, $z\approx\sigma$ (integrating the 12-6 potential over the discrete lattice results in a 10-4 potential⁴⁹ which has the minimum at $z=\sigma$ instead of $z=2^{1/6}\,\sigma$). The density oscillations diminish with distance from the substrate, and the local monomer density in the thin film approaches the bulk value at about $5$ monomer layers away from the substrate. At the liquid-vapor interface, the monomer concentration decreases rapidly to zero, with $\rho(z)$ following the characteristic $\tanh$ shape. If we consider the location of the polymers based on their center-of-mass positions, we find two local maxima near the substrate-liquid and liquid-vapor interfaces [Fig. 4(b)]. The surplus of chains near the interfaces can be understood by considering the conformation of chains in those regions, which are flattened in the directions tangential to the substrate and shrunk in the direction perpendicular to the substrate [Fig. 4(c)].^{21, 22}.

Polymer films on the curved substrated are created by placing the end-configurations of the flat film simulations such that the bottom film surface lies $\sigma$ above the top of the Gaussian bump. This distance is within the attraction range of the substrate-monomer interaction (see Sec. 2), so that the film slowly adsorbs to the curved substrate. We then let the systems relax for $5.5\times 10^{5}\,\tau$, which is typically long enough so that the films do not evolve anymore in time. Snapshots of the deposition process are shown in Fig. S2. Density distributions of fully deposited films are shown in Fig. 5 as colormaps in the Cartesian $xz$-plane, for films of nominal thickness $H=7.4\,\sigma$ adsorbed on substrates of height $A_{\text{g}}=20\,\sigma$ and widths $S_{\text{g}}=10\,\sigma$ and $S_{\text{g}}=2\,\sigma$, respectively. Similar to the flat films, the monomer density in these films is identical to the bulk density in regions sufficiently far away from the interfaces (yellow areas in Fig. 5). In a region within $\approx 5\,\sigma$ of the substrate, we observe the typical density oscillations (layering) ^{48, 49}. As is the case for flat films at the liquid-vapor interface, we also observe the characteristic decaying profile of the density in the direction perpendicular to the interface. In this region, we define the liquid-vapor interface as the region $\rho/\rho_{0}=0.1$ (magenta lines in Fig. 5).

To study the adsorbed polymer layer near the substrate in more detail, we plotted in Fig. 6 offset curves of $\rho$ at distances $n=1$, $2$, and $3\,\sigma$ as a function of arc length coordinate $s$,⁵⁰ as defined in Fig. 2. These curves demonstrate that lines of constant distance are also lines of (almost) iso-density. However, for the narrowest substrate, we observe a peak in density near $s=20\,\sigma$, where the substrate has the highest curvature $\kappa=0.47\,\sigma^{-1}$ (lowest positive radius of curvature $R\equiv\kappa^{-1}\approx 2.12\,\sigma$). Thus, the tendency for monomers to wet the substrate induces a local increase in monomer density near the region of highest positive substrate curvature. This effect can be understood by considering the local distribution of substrate monomers in the vicinity of positive curvature ($\kappa>0$): A monomer in those regions interacts with more surface beads within the cutoff distance $r_{\text{c}}$ of the attractive substrate-monomer interaction potential, which thus results in a local increase of the monomer density. Similarly, regions of negative curvature ($\kappa<0$) will experience the inverse effect with a decreased interaction when compared to a flat substrate.

The result for the radius of gyration colormaps are shown in Fig. 7 for the same parameters as in Fig.  5. We observe that in a region of $n\approx\sigma$ away from the substrate, the chains are stretched along the surface, resulting in a greater $R_{\text{g}}$ (thin red line in the colormaps). As for the flat films, we also observe that immediately after this layer we encounter a layer of reduced chain size (blue region). Again, we encounter a region of bulk properties between the substrate-liquid and liquid-vapor interfaces (yellow region). Consistent with the result presented above for the narrowest bump, we observe a dip in the size of the chains. Because we assign the chain property to its center-of-mass position, some bins inside the substrate are filled when chains are stretched along the top of the bump [red curved region in Fig. 7(b)]. From these data, it appears that regions of high positive (negative) curvature induce shrinking (stretching) of the adsorbed chains.

To support this finding, we plot in Fig. 8 the normalized radius of gyration as function of normal distance $n$ at some selected positions $s$. Here, we focus on a film of thickness $H=22.2\,\sigma$, such that we have full coverage on the top of the bump. Again, we find that the chains are stretched in regions $n\approx\sigma$ away from the substrate, whereas chains in regions of positive curvature (orange circles in Fig. 8) are significantly compressed when compared to chains near the flat regions of the surface (purple crosses in Fig. 8). Likewise, near the top of the substrate, where the curvature is negative (green squares in Fig. 8), we can find stretched chains even for regions of $n<\sigma$, since the the center of mass of these chains can fall in regions within the substrate.

For all substrates, we found that the liquid-vapor interface becomes flatter with increasing film thickness $H$, which indicates that polymers far away from the substrate do not experience the curvature anymore. To better understand (and possibly predict) the shape of the adsorbed polymer films, we constructed a continuum model, where the polymer film is treated as a homogeneous incompressible medium that fully wets the substrate. Within this framework, the free energy of the system is given by:

where $z(x)$ denotes the position of the liquid-vapor interface. This expression corresponds to the surface energy of the modulation $z(x)$ relative to a flat plane which is a smooth function (no terraces) of $x$. Thus, the liquid-vapor interface $z(x)$ of the deposited film can be computed by solving the variational problem of minimum surface (arc length in 1D) for the film with constraints (i) $L_{y}\int\left[z(x)-G(x)\right]\text{d}x=V$ to preserve the total film volume (incompressibility), and (ii) $z(x)\geq G(x)+h_{\text{min}}$ (full wetting of the film). These conditions are imposed by penalty factors in the free energy functional:

and

where $\Theta$ is the Heavisde function, and $\mu,\nu$ are the penalization parameters for the fixed volume and minimum thickness conditions, respectively. Since the underlying substrate $G(x)$ is a symmetric function of $x$, the solution will be symmetric also.

We convert the problem to its Weak Form, which allows greater flexibility in handling complex boundary conditions and constraints and solve it on an one-dimensional grid of spacing $\Delta x=0.1\,\sigma$ (roughly $\approx 1000$ points) with no-flux boundary conditions at $x=-L/2$ and $x=L/2$. Table S2 specifies the parameters used in the numerical solution. Figure 9 shows the comparison of the liquid-vapor interfaces obtained from MD simulations (blue dashed lines) and from solving the variational problem (red solid lines). We compare values of equal film thickness, assuming a constant monomer number density for the MD simulations $H^{*}=N/(\rho_{0}L_{x}L_{y}$). Note that local variations in the monomer density $\rho(z)$ [Fig. 5(a)] lead to a small mismatch between $H^{*}$ from simulations and our continuum calculations. Further, it is important to mention that polymer films with $H^{*}=5.9\,\sigma$ were unstable in the MD simulations with a narrow bump $S_{\text{g}}=2\,\sigma$, as the polymers dewetted the region of high curvature near the bump. In a similar manner, for the continuum model, enforcing the minimum thickness condition was not possible for this film thickness, resulting in a liquid-vapor interface that crossed the solid substrate. For the rest of the cases, the parameter $\nu$ was chosen (see Table S3) such that the solution converged and the condition was met. The resulting shapes of the liquid-vapor interface from MD simulations and the continuum model are in qualitative agreement: As the film thickness increases, the information from the surface modulation is increasingly lost, until the liquid-vapor interface becomes completely flat for sufficiently large thicknesses.