Influence of Chain Stiffness, Grafting Density and Normal Load on the Tribological and Structural Behavior of Polymer Brushes: A Nonequilibrium-Molecular-Dynamics Study

We have performed coarse-grained molecular-dynamics simulations on both flexible and semiflexible multi-bead-spring model polymer brushes in the presence of explicit solvent particles, to explore their tribological and structural behaviors. The effect of stiffness and tethering density on equilibrium-brush height is seen to be well reproduced within a Flory-type theory. After discussing the equilibrium behavior of the model brushes, we first study the shearing behavior of flexible chains at different grafting densities covering brush and mushroom regimes. Next, we focus on the effect of chain stiffness on the tribological behavior of polymer brushes. The tribological properties are interpreted by means of the simultaneously recorded density profiles. We find that the friction coefficient decreases with increasing persistence length, both in velocity and separation-dependency studies, over the stiffness range explored in this work.


Connectivity
A linear chain consisting of N beads represents each polymer chain within the polymer-brush system. Each bead in the chain may be representing more than a single monomer of a real polymer chain. Adjacent beads in each chain are bonded using a finite extendable nonlinear elastic (FENE) potential [15]. Each chain is attached by one of its ends to the substrate using a permanently tethered bead (red and blue beads in Figure 1), and the rest of the beads of each chain are free to move and interact with all other polymer and solvent beads via the ε pp and ε ps potentials specified below. Single spherical beads (cyan beads in Figure 1) represent the solvent particles. All of the simulations were performed for brush-against-brush systems in a planar setup. Periodic boundary conditions have been applied only along the axes parallel to the sliding surfaces.

Excluded Volume
The excluded volume and attractive interaction between polymer and solvent beads is modeled following Soddemann [49] and Dimitrov et al. [41] in their previous works. Within this model, the Lennard-Jones (LJ/12-6) potential is smoothly truncated at a cut-off distance r c = 1.5. This is achieved by using the LJ/12-6 potential only up to its minimum and shifting to some desired depth (ε pp , ε ps ), continued from its minimum to zero with a potential having a cosine form. This provides a potential that is both continuous and also has a continuous derivative. Thus, at small distances, the model potential corresponds to the purely repulsive, so-called shifted Weeks-Chandler-Anderson (WCA) potential [50,51]: where the index (αβ) stands for the different types of pairs: polymer-polymer (pp), polymer-solvent (ps) or solvent-solvent (ss), and we choose the range parameters to be identical, σ pp = σ ps = σ ss . Scales for length and energy (and temperature) are chosen such that both σ pp = 1 and ε LJ = 1 (Boltzmann's constant k B = 1). At intermediate distances, the interaction is purely attractive and parameterized by a cosine potential of the form: while U LJcos (r) = 0 for larger distances. As explained in detail elsewhere [41], we have used a value of c 1 = 3.173 and c 2 = −0.856 to ensure that the potential (U LJcos ) and its derivative are continuous both at the potential minimum, r/σ αβ = 6 √ 2, and at the potential cut-off, r/σ αβ = 3/2. For our simulations, we have used the purely repulsive potential for polymer-polymer (ε pp = 0) and solvent-solvent interactions (ε ss = 0), while the polymer-solvent interaction is modeled using ε ps = 0.4. This combination ensures a good-solvent quality [41,52].

Temperature, Deformation and Wall
The temperature was maintained constant by rescaling peculiar velocities [53] of all, except tethered, beads. The macroscopic velocity due to shearing is thus taken into account in the definition of peculiar velocities. We have used a profile unbiased thermostatting (PUT) scheme where the center-of-mass velocity for each set of beads residing in layers parallel to the wall is computed and used to define the 'bias' velocity. This bias velocity is then subtracted from the velocities of individual beads to yield a thermal velocity for each bead. The thermal velocity is rescaled to the desired value, and subsequently, the bias velocity is added [53,54]. To ensure that polymer and solvent beads do not cross the grafting surface, a purely repulsive LJ/9-3 potential was used. The motivation for using a wall potential at the grafting surface was to reduce the computational cost. Another approach used in nonequilibrium-MD (NEMD) studies is to thermostat [21,22] the explicit wall, rather than the confined fluid. Bernardi et al. [55] have shown in their study that direct thermostatting of confined fluids under shear strongly influence their behavior. It is important to mention here that their study was carried out at a very large shear rate,γ = 1.0. In our study, we have applied shear rates orders of magnitude smaller with Weissenberg numbers well below unity (except at the maximum velocity, v = 0.05). Direct thermostatting of the confined fluid in NEMD studies of polymer brushes is fairly common.

Semiflexibility
For chain-stiffness modeling, we have used the following discrete version of a bending potential: where k b is the dimensionless bending stiffness and u i denotes the unit segment vector connecting the i-th with the (i + 1)-th bead along the polymer chain. The bending stiffness k b and persistence length L p are interrelated via L p = −b 0 / ln L(k b ) with a bond length b 0 ≈ 1, which can readily be estimated from the already stated intramolecular potential at the given temperature, and the Langevin function

Simulation Details
The dynamics of such a molecular model at constant temperature were explored by solving modified (thermostatted) Newton's equations with the explicit solvent. We have used LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [57], an open-source code, for all simulation work. The calculated chain trajectories provide complete information about chain alignment, density profiles and the stress tensor. All simulated quantities reported in the current study are given in terms of LJ units [54,58].
We have carried out simulations for the brush-against-brush model system described in Figure 1. Three types of beads: non-tethered, tethered and solvent, were used. The simulations were performed on randomly grafted polymer chains on flat surfaces. The system consists of M = 20 chains on each tethering surface, while each linear chain is composed of N = 30 beads. It is important to mention here that the critical grafting density (ρ * = 1/πR 2 g ) [20] for chains having N = 30 beads is ρ * ≈ 0.025 (in LJ units). We have considered four grafting densities: ρ = 0.015 (less than critical grafting density), ρ = 0.025 (equal to critical grafting density), ρ = 0.075 (three-times the critical grafting density) and ρ = 0.15 (six-times the critical grafting density). The simulations were carried out for different values of chain stiffness, k b = 0, 1 and 3. The total number of beads in the simulation box is such that the number density of all beads is approximately 0.8 at each separation between walls.
We equilibrated all systems carefully first with 5 × 10 5 integration steps with ∆t = 0.001, followed by 6 × 10 6 and 10 7 steps with ∆t = 0.002 for ρ = 0.075, 0.15 and ρ = 0.015, 0.025, respectively. Data were then extracted by additional 10 7 time steps with ∆t = 0.0025. A reduction of separation distance by one (LJ unit) was achieved as follows. An amount of solvent beads compatible with unchanged number density 0.8 at the new separation distance was randomly removed from the system. One grafting surface was kept fixed and the other was moved with a constant velocity, v = 0.01, for a duration of 100,000 steps at ∆t = 0.001. At each separation D between the polymer-chain-bearing surfaces, the polymer-brush system was allowed to equilibrate for 3 × 10 6 time steps (10 6 steps at ∆t = 0.001 followed by 2 × 10 6 steps at ∆t = 0.0025).
These equilibrated systems were used to run NEMD simulations. Steady shear was applied by moving the tethered beads on both sides with the same prescribed velocity in opposite directions, keeping the separation constant during each run of given shear velocity [59]. At each separation and velocity, the stress tensor was calculated using the Irving-Kirkwood expression [54,58,60]. The temperature was maintained constant at T = 1.2 using peculiar-velocity rescaling as discussed in the previous section. We have performed NEMD simulations keeping the separation between grafting surfaces constant. The Gattinoni et al. [61] study of confined fluids showed that a constant separation between walls may lead to results that differ from those obtained at constant load, especially at small wall separations. Experimental studies of the tribological behavior of polymer brushes have been performed subjected to both constant wall separation [7][8][9]27] and constant load [23,24,62]. For the present work, we have chosen to operate at constant separation and measured the load, rather than the opposite.
Two kinds of studies were carried out: (i) separation-dependent: different separations while maintaining a fixed velocity (Table 1) and (ii) velocity-dependent: different velocities while keeping the separation between tethering surfaces constant ( Table 2). For the velocity-dependency studies, the simulations were carried out over a range of shear velocities v = 0.00015-0.05 at fixed separations. At all velocities, data during the first 500,000 time steps (∆t = 0.001) belonging to the startup phase were ignored for the analysis of stationary quantities. Subsequently, the simulations were performed for different numbers of time steps at different velocities with ∆t = 0.0025. At the maximum velocity v = 0.05, simulations were run for 200,000 time steps following a startup period comprising half of this period of time. Similarly, at the slowest velocity v = 0.00015, simulations were run for 7 × 10 7 time steps following the startup phase of 3.5 × 10 7 time steps. Simulations for each separation (D) and velocity (v) were repeated using 10 different initial configurations of randomly grafted polymer chains. The mean values from these runs are reported with standard deviations in Section 3.  [20,59], with N = 30.

Equilibrium
The equilibrium-brush height h of a single brush layer we define by the normal distance from the grafting surface where the brush begins to interact with the opposing brush, as seen from the time-averaged density profile. To be precise, we define h where the density profile falls below 1% of its maximum value close to the tethering surface. Using this criterion, the opposing brushes basically do not interact for D ≥ 2h. Other definitions using moments of the density profile are certainly possible and common, as well. Figure 2 shows the effect of grafting density on the equilibrium-brush height (h) of flexible polymers as calculated from the density profiles for each different grafting density at a large separation between the brush-bearing surfaces, while Figure 3 shows the density profiles for both flexible and semiflexible polymer-brush systems when the chains grafted on opposite surfaces were not in contact with each other. These density profiles were used to calculate and study the effect of chain stiffness (k b ) on the equilibrium-brush height (h).  Equilibrium-brush height (h) against chain stiffness (k b ) is plotted for brushes at grafting densities ρ = 0.075 and ρ = 0.15 in Figure 4. The results reveal an increase in the equilibrium-brush height with increasing chain stiffness for different grafting densities in agreement with the experiments of polymer brushes in solvent mixtures [63]. Following [64], for a planar brush under good-solvent conditions with L L p , one expects for the equilibrium-brush height within a Flory-type theory: where L = (N − 1)b 0 is the length of the fully extended polymer chain, b 0 the mean bond length, ρ the tethering density and τ a dimensionless solvent-quality parameter, i.e., a positive k b -independent constant of order unity under good-solvent conditions, and R 0 is the undisturbed equilibrium size of the free, semiflexible polymer. In the absence of excluded-volume interactions, our multi-bead chain with dimensionless bending stiffness k b should exhibit the statistical properties of a classical worm-like chain (WLC), characterized by mean square end-to-end distance R 2 ee = 2L p L{1 − L p (1 − exp(−L/L p )]/L} and squared radius of gyration R 2 where the persistence length L p is related to the bending stiffness k b via the before-mentioned relationship L p = −b 0 / ln L(k b ). For the comparison with the MD data, it does not qualitatively matter if we use R ee or R g as the undisturbed equilibrium size R 0 in Equation (4). We use the latter in Figure 4 with τ = 1 and an overall prefactor of order unity (1.68) to fit the data. While the Flory predictions are in excellent agreement with our MD data when the calculated equilibrium size of the polymer is used, the WLC naturally fails to capture the effect of excluded volume on the size and brush height in the limit of flexible chains (k b = 0), because the flexible FENE chain model has already a finite L p for k b = 0. The Flory-type result in Equation (4) follows from minimization of the free energy of a single chain, F/k B T = h 2 /R 2 0 + τNφ, with respect to h, where φ = Nb 3 0 ρ/h denotes the h-dependent volume fraction within the planar brush. The limitations of the Flory approach were discussed in detail in [64]. Within the mushroom regime, where polymers do not overlap, the height h of the polymer layer should become independent of ρ and scale like the radius of a polymer, for which Flory theory predicts h (τb 0 L 2 R 2 0 ) 1/5 . This corresponds to employing φ = Nb 3 0 /h 3 in the free energy expression. The equilibrium behavior of the model brush in the absence of shear is characterized by the normal stress as a function of separation between opposing tethering surfaces. In Figure 5, we highlight both the significant effect of surface density on the normal stress and the apparent insensitivity of the normal stress on the intrinsic stiffness of chains. Note that tethering of chains does not prefer a particular orientation of the anchored segments. The equilibrium behavior is of relevance to rate the nonequilibrium behavior to be presented in the subsequent sections. Figure 5. Normal stress versus gap size (system with M = 20, N = 30 subject to good-solvent conditions, constant bead number density) for different bending stiffness and surface densities, k b = 0 (red), k b = 1 (green) and k b = 3 (blue), for a brush system at ρ = 0.075 and k b = 0 (yellow) for a more densely grafted system at ρ = 0.15. The results are from equilibrium MD simulations. The normal stress is basically unaffected by bending stiffness over the shown range of distances.

Effect of Shear Velocity at Various Separations, but Fixed Tethering Density
In Figure 6, the shear and normal stresses are plotted for flexible chains (k b = 0) against shear velocity at different separations D = 40 (uncompressed brushes), D = 30 (marginally compressed brushes) and D = 20 (highly compressed brushes). The normal stress is found to remain constant at all velocities at each separation, but increases with increasing brush compression. We observed an increase in friction (shear) forces with increasing compression, i.e., increasing normal load at all investigated velocities. Similar trends were observed in previous experimental studies [10,24,25]. The coefficients of friction (ratio between shear and normal stress) also increase with increasing velocity and with increasing compression. Figure 6 further shows the density profiles at the corresponding separations for uncompressed brushes, marginally compressed brushes and highly compressed brushes, where interpenetration is observed. At a separation of D = 40, the grafted opposite polymer chains are not interacting with each other, and there is a layer of solvent in between the polymer outer surfaces, leading to a tribological behavior that is dominated by the friction behavior of the solvent, implying a very low or immeasurable friction coefficient [23]. As the surfaces are brought closer, the opposite chains interact with each other, and the coefficient of friction increases. With further compression, the penetration between opposing chains increases, resulting in a further increase in friction.

Effect of Grafting Density and Separation at Fixed Velocity
Here, we have run the simulations at different separations D between polymer-bearing surfaces at a constant velocity, v = 0.001, applied in opposing directions on tethered beads on both sides. Brush-against-brush systems were studied at various grafting densities, ρ = 0.015, 0.025, 0.075 and 0.15 (LJ units) for a fixed M = 20 and N = 30. As discussed, the critical grafting density for linear chains having N = 30 beads is 0.025. We hence probed the mushroom regime (ρ = 0.015), the system at critical grafting density, ρ = 0.025, and the brush regime (ρ = 0.075 and 0.15). The separation D was varied for each density, such that h < D < 2h for the respective case, i.e., over a similar compression.
The stresses calculated at different grafting densities are quantitatively very different. Figure 7 summarizes the results for for shear stress σ xz against normal stress σ zz calculated for different separations for systems with surface densities ρ = 0.015, 0.025, 0.075 and 0.15, upon applying a shear velocity v = 0.001 in opposing directions on tethered beads on each wall. The shear stress against normal stress data exhibit a linear dependency. Fitting a straight line (considering the error in each value) [65] to each curve, the coefficient of friction was calculated from the slope of the curve, as was also done in experiments [23]. The coefficient of friction thus obtained from the slope of the curves is plotted against the corresponding grafting densities in Figure 8. We observe a decrease in the coefficient of friction with increasing grafting density, while this coefficient tends to become insensitive to grafting density within the brush regime. There is a stark decrease of the coefficient of friction when we enter from the mushroom regime (ρ = 0.015) into the brush regime (ρ = 0.075 and ρ = 0.15). This can be explained as follows. With the increase in grafting density (ρ), chains are more stretched out due to the excluded volume effect to support higher normal load and ensure a thin solvent lower at the interface, even at higher compression. This thin solvent layer ensures lower friction at higher compression. The decrease in the friction coefficient with increasing ρ was observed in experiments [26], where at a moderate load/pressure, the system with highest grafting density exhibits the lowest friction. As the applied load increases, a "transition load" was observed at each grafting density in the experimental work [26]. Above the transition load, the friction was found to be increasing with the increase in grafting density. Over the range of separations D we studied, we did not see any "transition load" in our simulation work.  Figure 7 (so-called Approach A) and not for a particular velocity (so-called Approach B).

Effect of Shear Velocity for Two Grafting Densities at Fixed Separation
Here, normal and shear stresses were calculated for various velocities v at a fixed separation D = 30. The simulations were performed for systems with different grafting densities ρ = 0.075 and ρ = 0.15 (Figure 9), but for the same number of grafted chains M = 20 on each side having N = 30 beads. We observed an increase in normal and shear forces with increasing grafting density at all shear velocities, consistent with an expected greater osmotic pressure, and with a more significant contribution of chain-chain interactions, respectively. In the shear stress versus velocity curve (Figure 9b), there is a slight increase of the shear stress with increasing velocity for brush-brush systems at relatively low velocities. Similar trends were observed in experimental studies [24,25] and are designated as a "boundary-regime"-like behavior. For the corresponding grafting densities and separation (D = 30), we have plotted density profiles of solvent and polymer chains. Here, we observe that at lower grafting density (ρ = 0.075), the polymer chains are marginally compressed or uncompressed, but at higher grafting density (ρ = 0.15) for the same separation (D = 30), the brushes are more compressed; and interpenetration between opposite chains is observed, due to larger equilibrium-brush height with increasing grafting density ( Figure 2). This higher interpenetration between chains causes higher friction that we observe at higher grafting density. These results are in agreement with Galuschko et al. [20], where an increase in friction with an increase in grafting density at constant separation distance was observed for a brush-against-brush system with a dissipative-particle-dynamics (DPD) thermostat.
To summarize, we explored two different approaches to study the effect of grafting density on the tribological behavior of polymer brushes. In the first approach (A), we kept the separation constant and varied the shear velocity for systems with different grafting densities. In that case, the coefficient of friction was calculated from the ratio between shear and normal stresses at each shear velocity for the different systems. Here, we observed an increase in the coefficient of friction with increasing grafting density at each shear velocity. This increase can be attributed to the fact that more polymers enter the penetration (interaction) zone for a system with higher grafting density. In the second approach (B), we maintained the shear velocity constant, but varied the separation D between opposing grafting surfaces for systems with different grafting densities over a range, such that h < D < 2h holds. Normal and shear stresses were calculated at each separation. The coefficients of friction were calculated for each system from the slope of the shear stress against normal stress curve. Following this approach, we observed a decrease in friction with increasing grafting density. Here, we would like to make an attempt to rationalize the two different results obtained via A and B for the coefficient of friction. From our simulations, it seems that the normal stress is independent of v and increases with decreasing separation distance D, i.e., σ zz (v, ρ, D) = σ bulk zz + d 2 (ρ)/D, where σ bulk zz would correspond to the bulk value (D → ∞), and the force density d 2 (ρ) depends on ρ alone. The shear stress σ xz seems to increase linearly with v, σ xz (v, ρ, D) = [c bulk xz + d 1 (ρ)/D]v, parameterized by coefficients c bulk xz and d 1 (ρ). From Figure 9, we infer that both d 1 (ρ) and d 2 (ρ) tend to increase with increasing ρ, and Figure 10 seems to suggest that c bulk xz is extremely small or even vanishes. For approach A, the definition of the coefficient of friction is: while for B, we look for a linear relation between σ xz and σ zz as we vary D. Inverting the relation σ zz (v, ρ, D) for fixed ρ and v to find D(σ zz , v, ρ) and inserting into the corresponding expression for σ xz (v, ρ, D), we obtain such a linear relationship, σ xz = const. + CoF B σ zz with: First of all, CoF A and CoF B cannot be expected to be identical. Our results for approach B have shown that d 1 (ρ)/d 2 (ρ) tends to decrease with increasing ρ, while A revealed an increasing CoF A with increasing ρ. Assuming c bulk xz = 0, we can rationalize our results if we tentatively approximate d 2 (ρ) ≈ ρd 1 (ρ) and assume small enough values of ρ/σ bulk xz D. For sufficiently strong compression, the bulk contributions become irrelevant, and CoF A and CoF B lead to identical results.

Effect of Chain Stiffness and Separation at Fixed Velocity and Surface Density
The simulations were performed for a brush-against-brush system having specifications, M = 20, N = 30 and ρ = 0.075. Such a model brush system was studied for different chain-stiffness values k b = 0, 1 and 3. A constant shear velocity v = 0.005 was applied on tethered beads on both sides, but moving in opposite directions. The simulations were run over different separations D, and stresses were calculated at each D. Figures 10a-b shows normal and shear stresses against separation D curves obtained for chains with different persistence lengths. We observe a lower shear stress at a given normal stress for systems with chains of higher stiffness. The normal stress versus distance plot does not show any effect of chain stiffness, while the shear stress versus distance graph (Figure 10b) highlights a decrease of the shear stress with increasing chain stiffness at identical separations. This can be explained by the fact that even though at the same separation the system with higher stiffness shows higher interpenetration between opposite chains, brushes having stiffer chains align themselves along the shear direction easily, so the friction decreases rather than increases.

Effect of Velocity and Chain Stiffness at Fixed Separation and Surface Density
Here, normal and shear stresses are plotted against shear velocity at a separation of D = 30 and surface density ρ = 0.075 for two different values of chain stiffness k b = 0 and k b = 3 in Figure 11. We observe an increase in normal stress with the increase in chain stiffness, whereas the shear stress decreases with increasing chain stiffness. These two observations lead to the conclusion that the coefficient of friction decreases with increasing persistence length over the range of stiffness studied.
Here, again, we have observed boundary-regime-like behavior [24,25] at lower velocities in the shear stress versus shear velocity representation (Figure 11b). The corresponding density-profile plots show that flexible (k b = 0) polymer chains are marginally compressed at separation D = 30, but for stiffer chains (k b = 3), interpenetration is observed for the same separation D = 30.

Comparing Flexible and Semiflexible Brushes of Identical Height
In Figure 4, we have observed that different combinations of grafting density (ρ) and chain stiffness (k b ) can result in a similar equilibrium-brush height (h). This is also clear from Equation (4). For example, two different combinations, first ρ = 0.075 and k b = 3 and second ρ = 0.15 and k b = 0, result in a basically identical equilibrium-brush height, h ≈ 18. In Figure 12a-c, the normal and shear stresses are plotted against shear velocities for these two different cases discussed earlier at the same separation, D = 36. The corresponding concentration profiles in Figure 12d-e confirm the similar equilibrium-brush heights. We have observed higher normal and shear stresses at all shear velocities for the combination where the chains are more flexible and the grafting density higher (though the equilibrium-brush height is very comparable in both cases). The height h and dimensionless ratio h/D alone do therefore not uniquely determine the tribological properties, but the chain stiffness and grafting density serve as additional design parameters.

Conclusions
NEMD simulations were performed on a brush-against-brush system with explicit solvent beads. The effect of grafting density and chain stiffness on the tribological behavior of polymer brushes was studied. When the separation between polymer-bearing surfaces is maintained constant and the friction behavior of brushes monitored over different shear velocities, the friction is found to increase with increasing grafting density. However, when the velocity is kept constant and the friction behavior was studied over different separations, the coefficient of friction is found to decrease with increasing grafting density. This can be explained in terms of the load-bearing capacity of polymer brushes with higher grafting density, and we provided some rationale to interpret the differences at weak compression, while we expect agreement between both measures at sufficiently strong compression. The friction coefficient furthermore depends clearly on chain stiffness and is found to decrease with increasing persistence length, both in the velocity-and separation-dependency studies, over the stiffness range explored.
This work was undertaken to help establishing design rules for polymer-brush assisted lubrication and to shed some light on the effect of tuning parameters, such as grafting density and persistence length on the tribological properties of brushes sheared against each other. Semiflexible polymer brushes are particularly relevant in the context of DNA-coated nanoparticle crystallization. A natural progression from this work includes exploring the effect of molecular weight and solvent quality to better understand the interplay of these additional tuning parameters on the tribological behavior of the model polymer brush. The solvent quality can be varied by changing the attractive part of interactions between solvent and polymer beads. Furthermore, constant-load simulations would be helpful to address the experimentally observed transition load.