Surface wave excitations and backflow effect over dense polymer brushes

Polymer brushes are being increasingly used to tailor surface physicochemistry for diverse applications such as wetting, adhesion of biological objects, implantable devices and much more. Here we perform Dissipative Particle Dynamics simulations to study the behaviour of dense polymer brushes under flow in a slit-pore channel. We discover that the system displays flow inversion at the brush interface for several disconnected ranges of the imposed flow. We associate such phenomenon to collective polymer dynamics: a wave propagating on the brush surface. The relation between the wavelength, the amplitude and the propagation speed of the flow-generated wave is consistent with the solution of the Stokes equations when an imposed traveling wave is assumed as the boundary condition (the famous Taylor’s swimmer).


DPD methodology
DPD is a coarse-grained Molecular Dynamics method. It was introduced by Hoogerbrugge and Koelman in 1992 [1] to simulate isothermal Navier-Stokes equations and it has been recently recognized that it can be successfully applied to any mesoscopic scale [2]. In the spirit of DPD, the simulated particle does not correspond to a single molecule, but rather to a significant large cluster of them. Each of the N point-like particles evolves in time according to Newton's equations that we solve by using the velocity Verlet algorithm [3]. Observables are then calculated as averaged all over time configurations. The force F ij of Eq. 1 on each particle has three contributions [4] resulting from the coarse-graining procedure: a conservative component F C ij , a dissipative component F D ij and a stochastic one F S ij . All these forces are pair-wise, to guarantee momentum conservation, and limited to a cut-off radius r c .
The conservative soft-core repulsive force has the following expression: with r ij = r i − r j is the vector distance between the i-th and j-th particle, r ij = | r ij | andr ij = r ij /r ij . Such repulsive potential models the averaged fast microscopic length and time scales and was first derived by coarse-graining particles interacting via a Lennard-Jones potential [5]. The indexes α, β indicates the particle type (solvent or polymer, in our case). The constant a α,β measures the force between two completely overlapping particles and it is casted from the compressibility of the modeled fluid. The two other forces account for the loss of details in the coarse-graining procedure. The dissipative force has the form: where the standard choice in literature for the "weight function"w D (r ij ) is: Equation (3) introduces a friction among particles proportional to the relative velocity v ij = v i − v j Thermal fluctuations are added via the stochastic force: where σ is a constant (related to temperature), w S (r ij ) is a weight function and θ ij is a random number extracted from a gaussian distribution with zero average and indipendent in time and among particle pairs: θ ij (t)θ lm (t ) = (δ il δ jm +δ im δ jl )δ(t−t ). Momentum conservation requires θ ij = θ ji .
Two additional constraints, w D (r ij ) = [w R (r ij )] 2 and γ = σ 2 2k B T , guarantee that the probability distribution respects the statistics of the NVT ensamble [6].
To mimic the brush chains we add a finite extensible nonlinear elastic potential (FENE) [7] F f ene which account for the neighbour particle connectivity. In Eq. (6), r eq is the equilibrium distance between neighbour monomers and R is the maximum allowed extensibility. Grafted points are randomly chosen from a uniform distribution and located on a flat surface. Note that polymers are non-ideal since they interact via an excluded volume potential.

DPD and physical units
We chose the DPD units such as r c = 1, m i = 1, t DP D = 1 and k B T = 1. Following [4], we fix the solvent-solvent interaction parameter a SS relating it to the adimensional compressibility of water κ −1 T . A comparison between the equation of state of a simulated DPD system and the experimental data (κ −1 T (water) = 16 at room temperature) suggests a SS = 75k B T /ρ = 25. We assume that the polymer-polymer interaction parameter has the same value (e.g. a SS = a P P ), while we select a smaller value for the solvent-polymer parameter a SP = 20 (good solvent conditions). The noise amplitude is fixed to σ = 3. For the FENE potential, finally, we use r eq = 0.86, R = 1 and k = 50.
In DPD the choice for the physical units depends on the level of coarse-graining is desired. Our reference system is the endothelial glycocalyx, therefore we fix the physical lengthscale l phys relating the spacing between different filaments of the glycocalyx network, d glyco = 20nm [8] with the average distance between anchor points of our brush, d graf t = 1/σ graf t = 0.82: We underline that with this choice also the brush thickness is in the range of endothelial glycocalyx (100nm ö1000nm [9]).
As for the physical mass scale m phys and physical time scale t phys , we exploit the comparison between viscosities and between energies. The viscosity of water η phys is η phys = 10 −3 P a s at 300K and the DPD viscosity η DP D in a bare channel, estimated from the slit-pore velocity profile relation η = ρA(L z /2) 2 /2v max , where v max is the maximum velocity (see Sect. 1.2 for definitions of the other parameters), corresponds in our case to η DP D = 0.84. The physical energy scale is defined as k B T phys , with T phys = 300K. Thus we can write down two relations from which we extract

The system geometry
We simulate a parallelepiped box of sides L x = 30, L y = 5 and L z = 50. On the x-and y-axes we impose periodic boundary conditions, while at z = 0 and z = L z two impenetrable parallel walls of infinite mass are set. At the wall we use the so-called "bounce back reflection"conditions [10] in which all the velocity components are reversed: Such conditions assure no-slip boundaries. The integration time step is ∆t = 0.02. We set the density ρ = N/(L x L y L z ) = 3 [4], therefore N = 22500. We fix the number of monomers per chain n = 40 and the grafting density σ graf t = 1.5, defined as σ graf t = N ch /(L x L y ), with N ch the number of chains composing the brush. For the sake of ease and computational time, we attach polymers only at the bottom wall.
In order to produce a parabolic velocity profile inside the channel (e.g. along the z direction) a constant acceleration A = Ax is applied to all fluid particles: In turn, these fluid particle exchange momentum with the chains, dragging them. Different values of A allow us to probe different dynamic regimes.

Brush equilibrium properties
Polymer brushes immersed in a liquid at equilibrium (i.e. in the absence of flow) have been thoroughly investigated in the past [11,12,13,14]. Here we recall that the brush equilibrium conformation results from the balance between configurational entropy, that tends to make chains visit the whole available space, and excluded volume interactions, avoiding contact between monomers. The brush can be properly described by its profile ρ(z), indicating the probability distributions of finding a monomer at distance z from the grafting wall. We show in Fig. 1 ρ(z) at fixed σ graf t = 1.5 for several chain lengths (n = 20, 25, 30, 36, 40, 45). In line with previous theoretical results [12] and numerical studies [13,15], the higher the volume fraction, the more step-like the density profile is. A definition of the brush height h in term of the first moment < z > of the density profile depends on the shape of ρ(z). Assuming that each chain is completely elongated so that free ends are all at the maximum distance from the grafting wall, the profile ρ(z) is step-like (Alexander model ) and the brush height can be defined as h = 2 < z >. A less rigid theory (self-consistent field theory) hypothesizes that monomers of a same chain are distributed as a random walk, therefore the free end position ranges uniformely over the whole brush thickness. In such a case the profile is parabolic, ρ(z) = ρ(0)(1 − z 2 /h 2 ), and the brush height h is h = 8 3 < z >. For the simulations discussed in the manuscript we select n = 40, so that the evaluated ρ(z) is in between the step and the parabolic shapes (see Fig. 1). In the following, we define h following the definition of a step distribution, e.g.   Fig. 2 we observe that E is increasing with the compression. We now notice that flow produces a compression of the brush (see Fig. 3). Equating the compression imposed at equilibrium with the one produced by the applied flow (shaded areas in Fig. 2) we conclude that the brush elastic response is modified by a change in W i.

The Weissenberg number
Instead of using A as a measure of the flow intensity we prefer to introduce an adimensional quantity, the Weissenberg number W i := t brush /t f low , where t brush should be a characteristic structural time of the unperturbed brush and t f low should provide an estimate of the typical time scale associated to the flow. We evaluate t brush as the relaxation time obtained by the decay of the autocorrelation function, averaged over all distinct polymers, of the end-to-end vector amplitude R ee = | r 0 − r n |, where r n and r 0 are the position vectors of, respectively, the last particle and the anchor of one chain inside the brush. According to such definition, the relaxation time t brush assumes a dependency on the grafting density and on the degree of polymerization, thus representing a time scale of the whole brush rather than an isolated single-chain property.

Videos
We show in movies SI1 and SI2 the dynamics of a single chain inside the brush for two different W i, one giving flow inversion (movie SI1) and the other with no flow inversion (movie SI2): the recursive motion is constituted by an alternation of elongation-stretching-recoiling in both cases, thus it cannot be used as sufficient explanation for the backflow phenomenon. The motion of the whole brush under flow has been recorded in movies SI3 and SI4. In the first case no collective motion is detectable and no backward flow is observed in the velocity profile, while in the second case there is concurrence of a surface wave and of flow inversion.

Finite size effect
We have investigated the influence of finite size effect increasing the channel length L x at fixed acceleration A: the wave persists. For a larger channel, 30 < L x < 40, the wave increases in wavelength (λ = L x ) and decreases in frequency, as shown in Fig. 4, such that the propagation speed keeps constant. The wave amplitude decreases and, in accordance with the Taylor's relation (see Letter, Eq. (1)), v min decreases too.
Another test has been attempted increasing, up till the double, both L x and L z (L x = 60  and L z = 100). For such a system we found again, for W i = 464, a surface wave, with double amplitude and double wavelength (see movie SI5). However, investigating channels up to L x = 360, L y = 20 and L z = 150, we found that the amplitude of the surface wave does not keep on growing upon increasing L x : in fact, its value seems to reach a plateau. We interpret this non-linear behaviour as a consequence of the finite extensibility of the polymers, e.g. due to the higher and higher energy required to significantly compress and stretch the brush.