Adhesion and Separation of Nanoparticles on Polymer-Grafted Porous

Subscriber access provided by Monash University Library

Article

Adhesion and Separation of Nanoparticles on Polymer-Grafted Porous Substrates Kolattukudy P. Santo, Aleksey Vishnyakov, Yefim Brun, and Alexander V. Neimark Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.7b02914 • Publication Date (Web): 15 Sep 2017 Downloaded from http://pubs.acs.org on September 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Adhesion and Separation of Nanoparticles on Polymer-Grafted Porous Substrates Kolattukudy P. Santo1, Aleksey Vishnyakov1, Yefim Brun2 and Alexander V. Neimark1,* 1

Department of Chemical and Biochemical Engineering, Rutgers University, Piscataway, NJ 08854, USA; 2

DuPont Central Research & Development, Wilmington, DE 19803, USA. *- corresponding author, [email protected]

Abstract: This work explores interactions of functionalized nanoparticles (NP) with polymer brushes (PB) in a binary mixture of good and poor solvents. NP-PB systems are used in multiple applications, and we are particularly interested in the problem of chromatographic separation of NPs on polymer-grafted porous columns. This process involves NP flow through the pore channels with walls covered by PBs. NP-PB adhesion is governed by adsorption of polymer chains to NP surface and entropic repulsion caused by the polymer chain confinement between NP and the channel wall. Both factors depend on the solvent composition, variation of which causes contraction or expansion of PB. Using dissipative particle dynamics simulations in conjunction with the ghost tweezers free energy calculation technique, we examine the free energy landscapes of functionalized NPs within PB-grafted channels depending on the solvent composition at different PB grafting density and polymer-solvent affinity. The free energy landscape determines the probability of NP location at given distance to the surface, positions of equilibrium adhesion states and the Henry constant that characterizes adsorption equilibrium and NP partition between stationary phase of PB and mobile phase of flowing solvent. We analyze 1 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

NP transport through a polymer-grafted channel and calculate the mean velocity and the retention time of NP depending on the NP size and solvent composition. We find that with the increase of the bad (poor) solvent fraction and respective PB contraction, NP separation exhibits a transition from the hydrodynamic size exclusion regime with larger NPs having shorter retention time to the adsorption regime with smaller NPs having shorter retention time. The observed reversal of the sequence of elution is reminiscent to the critical condition in polymer chromatography, at which the retention time is molecular weight independent. This finding suggests a possibility of the existence of an analogous special regime in nanoparticle chromatography, at which NPs with like surface properties elute together disregarding of their size. The latter has important practical implications: NPs can be separated by surface chemistry rather than by their size employing the gradient mode of elution with controlled variation of solvent composition. KEYWORDS: nanoparticles, polymer brush, adhesion, nanoparticle chromatography, dissipative particle dynamics

2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

1. Introduction Nanoparticles (NPs) are being used extensively in modern-day technologies on versatile roles, including various applications in diagnostics and therapeutics1-2 as carriers for drug delivery, in imaging as contrast agents,3 in biosensing,4-5 in electrodes of supercapacitors,6 in light emitting diodes (LEDs),7-8 in solar cells,9-10 and in nanoelectronic devices.11 Many of these applications involve interaction of nanoparticles with polymers or polymer brushes (PBs) grafted to solid substrates. PBs doped by NPs are used in fabrication of nanocomposites, sensors and biomedical devices12, colloidal stabilization, lubrication13 and manipulation of nanoobjects on surfaces.14 Emerging applications of PB-grafted substrates in NP separation and purification15 attract special attention to NP – PB interactions. Traditionally, techniques for NP separation stem from those employed for characterization and separation of polymers, including size exclusion chromatography (SEC),16 hydrodynamic chromatography (HDC),17 field flow fractionation (FFF)18-19 and size selective precipitation20 (SSP). These techniques are designed to separate NPs mainly according to their size and shape. However, many applications require NP separation by surface chemistry (that may be manifested, for example, by the degree of hydrophobicity or surface heterogeneity), especially when NPs are modified by specific functional groups (“ligands”).21-23 In polymer chromatography, separation by chemistry, e.g. chemical composition or microstructure of polymer chains, is achieved at the so-called critical point of adsorption (CPA).38 At CPA, the repulsive entropic effect of polymer chain confinement in the pores of solids substrate is compensated by the attractive enthalpic (adsorption) interaction.24 The balance of entropic and enthalpic factors results in molecular weight-independent elution of polymers. The corresponding separation regime, often called liquid chromatography at critical conditions 3 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(LCCC), has been applied to separate polymers of different architecture24 such as diblock and triblock copolymers, star polymers and combs at the isocratic and gradient elution modes.25-27 CPA can be achieved by varying the solvent composition which controls the solvent quality.28 In a good solvent, attractive polymer-substrate interactions are suppressed and steric interactions dominate favoring separation by the hydrodynamic size (size exclusion chromatography), when large molecules elute faster than their smaller counterparts. A similar order of elution occurs in hydrodynamic chromatography in capillary channels or packed columns17. In the other limiting case of strong adsorption achieved with a poorer solvent, elution occurs in the reverse order with larger molecules retained longer. CPA is achieved at an intermediate composition when retention becomes independent of the molecular weight. Recently, we demonstrated that the critical conditions of polymer adsorption exist in the chromatographic columns packed with non-porous particles, in which the separation occurs in the interstitial volume between particles (flowthrough channels).29-30 The present work studies the interactions between NPs and PBs in a binary solvent by means of dissipative particle dynamics (DPD) simulations. We explore the dependence of the free energy landscape of NPs in the vicinity of PBs on the particle size and solvent composition and predict NP transport and retention in PB-grafted pore channels. Our motivation is to determine whether the reversal of the elution order and CPA, found in polymer chromatography upon varying the solvent composition, can also be achieved in NP chromatography, thus leading to size-independent NP elution. Although this particular task has (to the best of our knowledge) never been tackled in the literature, the NP-PB interactions in general have been actively explored both experimentally and theoretically.


Page 4 of 47

Page 5 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Of our particular interest are theoretical studies of the free energy of NP adhesion and sorption by PB and its influence on NP mobility. NP-polymer interactions are complex and involve multiple scales that requires theoretical approaches31-32 and computer simulations33-44 on a coarse-grained supramolecular level. Theoretically, NP-polymer interactions were explored using self-consistent field theory (SCFT),31, 33, 35-37, 39 Brownian dynamics (BD),40 and dissipative particle dynamics (DPD)38 simulations. Milchev et al.45 compared the interactions of NPs with PBs and polymeric solutions. They noted that free energy of NP rises steeply as NP penetrates deeper into an expanded PB, which could not be explained solely by the dependence of polymer density on the distance to the substrate.46-47 Halperin et al.31 explored the dependence of the disjoining pressure between NP and PB on solvent quality using SCFT. The reduction of the solvent quality, which leads to eventual collapse of PB, reduces the associated free energy barrier and favors NP penetration and absorption within PB. Merlitz et al. in their MD study 48 observed that chain polydispersity reduces the free energy of NP immersion into PB. Zhang et al.40 considered the effect of attraction between NP and polymer using BD simulations. For larger NPs, the free energy has a minimum with respect to NP location, thus indicating a preferential depth of NP immersion into PB. Our group49 explored NP interaction with PBs using DPD and focused at the interplay between attractive forces caused by polymer adsorption at NP surface and NP attraction to the substrate, and entropic repulsive forces at varying grafting density, NP size, and solvent composition. We described two basic mechanisms of NP-PB interactions: NP adhesion at the PB exterior and NP immersion into PB with adhesion to substrate, and calculated the free energies of these adhesion states and the energy barriers between them. Somewhat similar observations were made by Hua et al.50 who modeled the distribution of NPs of different sizes in the vicinity of a semi-flexible PB-grafted support. At certain conditions, NP


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

concentration indicated an existence of a potential well near the brush support. CGMD simulation of Nasrabad et al.51 considered NP influence on PB structure. The authors modelled PBs grafted onto planar surfaces and in cylindrical pores. Strong adsorption of polymer chains onto NPs caused PB collapse. At a low grafting density, a strongly non-uniform PB outer boundary was observed with polymer chains clustered around the NPs that were located close to the surface. This phenomenon is qualitatively similar to the behavior of low density PBs in bad solvent. PB became more uniform when the polymer grafting density and NP concentration increased. Cao et al.52 modeled NP diffusion within a cylindrical channel with PB-grafted walls using DPD framework with Lennard-Jones quasiparticles. The authors varied the grafting density and the attraction between NPs and tethered chains and observed a transition from static NPs trapped by the chains to fast transport of NPs in the channel center. In the recent works from our group, we investigated the PB conformational transition and NP adhesion to PB in a binary solvent at varying solvent composition49, 53 using DPD simulations. The original ghost tweezers (GT) method was employed to calculate the free energy landscape of NP-PB interaction and analyze the equilibrium adhesion states and the energy barriers separating these states. The DPD model and the GT method proposed in refs.49, 53 constitute the methodological foundation of this work. The rest of the paper is structured as follows. Section 2 describes the methodological framework of our study: details of DPD simulations, parameterization of the coarse-grained model of polymer, solvent and nanoparticles, methodology of calculations of the free energy landscape of NPs near PBs by the ghost tweezers method, calculation of Henry constant of NP adhesion to PB in slit-like channels with PB-covered walls and modeling of NP transport through PB-grafted slit-shaped channels. In Section 3, we present the simulation results and discuss the 6 ACS Paragon Plus Environment

Page 6 of 47

Page 7 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

dependence of NP-PB interactions on the NP size and solvent composition, variation of the Henry constant and the retention time with the reduction in the solvent quality during NP flow through PB-grafted channels and demonstrate and quantify the observed transition from the sizeexclusion to the adsorption regime of NP separation which is accompanied with the reversal of the sequence of NP elution. Conclusions are summarized in Section 4.

2. Models and methods 2.1 Dissipative particle dynamics simulations In this work, we employed standard DPD implementation: molecules are dissected into soft quasiparticles (beads) each of which represents a group of atoms together. The beads interact with effective potentials parameterised to reproduce targeted structural and thermodynamic macroscopic properties of the material. The beads move according to Newton’s equations of motion while acted upon by a force, (C)

(B)

(D)

(R)

𝒇𝑖 = ∑(𝐅𝑖𝑗 + 𝐅𝑖𝑗 + 𝐅𝑖𝑗 + 𝐅𝑖𝑗 )

(1)

𝑖≠𝑗

where 𝒇𝑖 is the total force on the ith particle resulting from interactions with its neighbors. 𝒇𝑖 (C)

(B)

(D)

sums up the short-range conservative repulsion forces 𝑭𝑖𝑗 , bond forces 𝑭𝑖𝑗 , drag forces 𝑭𝑖𝑗 (R)

(C)

and random forces 𝑭𝑖𝑗 . 𝑭𝑖𝑗 is the harmonic repulsive potential between overlapping beads: (C)

𝑟

𝐅𝑖𝑗 = 𝑎𝑖𝑗 (1 − 𝑅𝑖𝑗) 𝒓̂𝑖𝑗 for 𝑟𝑖𝑗 ≤ 𝑅c and zero beyond the effective bead diameter 𝑅C . 𝑟𝑖𝑗 = c

|𝐫𝑖 − 𝐫𝑗 | and 𝑎𝑖𝑗 is known as the conservative repulsion parameter. Beads connected by bonds (B)

interact via harmonic bond forces 𝐅𝑖𝑗 = 𝐾B (𝑟𝑖𝑗 − 𝑟e )𝐫̂𝑖𝑗 where 𝐾B is the effective bond rigidity


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 47

and 𝑟e is the equilibrium bond length. The drag force and random force are implemented to maintain the Langevin thermostat, 54-55 (D)

(R)

𝐅𝑖𝑗 = −𝛾 𝑤 (D) (𝑟𝑖𝑗 )(𝐫̂𝑖𝑗 . 𝐯𝑖𝑗 )𝐫̂𝑖𝑗 ;

𝐅𝑖𝑗 = 𝜎𝑤 (R) (𝑟𝑖𝑗 )𝜃𝑖𝑗 𝐫̂𝑖𝑗

(2)

where 𝛾 is the friction coefficient, 𝜎 2 = 2𝛾𝑘𝑇 and 𝐯𝑖𝑗 = 𝐯𝑖 − 𝐯𝑗 is the relative velocity between ith and jth beads. The weight functions 𝑤 (D) and 𝑤 (R) are related as 𝑤 (D) = (𝑤 (R) )2 with 𝑟𝑖𝑗

𝑤 (R) (𝑟𝑖𝑗 ) = (1 − 𝑅 ) for 𝑟𝑖𝑗 ≤ 𝑅c and zero for 𝑟𝑖𝑗 > 𝑅c . 𝜃𝑖𝑗 is a random variable with Gaussian c

statistics. In this work, we use the conventional implementation of DPD: all beads have the same effective diameter Rc and the same intra-component repulsion parameter aII. The DPD model was parameterized in our previous works 49, 53 and effectively mimics polyisoprene natural rubber PB in benzene (good solvent) – acetone (bad solvent) binary solution. The bead size Rc is set to 0.71 nm. The bead density b equals 3𝑅c−3 and aII=42 kBT/Rc to effectively reproduce the overall compressibility of the system under consideration. 53, 56 DPD simulations were performed in reduced units with 𝑅c as the unit of length and 𝑘B 𝑇 as the unit of energy using LAMMPS57 software. Configuration snapshots were created using the Visual Molecular Dynamics (VMD) program.58 2.2 System setup In order to calculate the free energy landscape of NP-PB interaction, we use the ghost tweezers (GT) method with a simulation set-up similar to the one developed in our previous work 53 (Fig.1). The coordinate system is chosen so that the grafted substrate surface is parallel to XY plane and corresponds to z = 0. The PB is constructed by tethering linear polymer chains


Page 9 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

to the solid substrate at given grafting density. Each chain consists of 90 beads of type P that corresponds to polyisoprene molecular weight of 6120 Da. The substrate is formed by an array of immobile S beads arranged in a cubic order at a very high density of 19Rc-3, which effectively makes the substrates strongly repulsive towards all polymer or solvent beads. We consider two grafting densities, 𝛤 = 0.6 and 𝛤 = 0.36 chains per nm2. PB is immersed in the binary mixture of good (G) and bad (B) solvent components. Solvent composition (expressed as volume fraction of good component G in the solvent xG = NG/(NG + NB), where N is the number of beads of particular type) determines the PB conformation: the PB expands in the good solvent and contracts upon addition of the bad solvent. The parameters for polymer – solvent and solvent – solvent interactions are taken from refs.49, 53 to effectively represent PINR in acetone (bad solvent, ΔaBP= aBP - aII =7.0 kBT/Rc) and benzene (good solvent, ΔaGP = aGP-aII =1.5 kBT/Rc). The solvents are miscible with each other, ΔaGB=1.5 kBT/Rc. We also consider another system with the bad solvent component denoted as B2, which interacts with the polymer more favorably than B, ΔaB2P=4.0 kBT/Rc. NPs of spherical shape are composed of beads arranged on hexagonal simple closepacked lattice and linked by strong harmonic bonds. NP is formed by two types of beads. The two outer layers of NP consist of beads labeled N which favorably interact with polymer. The inner layers of NP are formed by the “core” beads of type C that strongly repel all other beads in the system to ensure that NPs are impenetrable to solvent and polymer. Because of high bond rigidity, NP maintained an approximately spherical shape. To study the effect of NP size we consider NPs of 3 different radii, 4𝑅c = 2.84 nm, 6𝑅c = 4.26 nm and 8𝑅c = 5.68 nm, each consisting of M = 740, 2496 and 5917 beads respectively. Short ligand chains are tethered to the NP surface to mimic NP functionalization. Each ligand chain is composed of six beads of type L. 9 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 47

Sidechains are distributed uniformly over the NP surface. In all simulations reported in this work, the ligand grafting density is set to 0.39 nm-2 which corresponds to 40 ligands for 𝑅NP = 4𝑅c , 90 ligands for 𝑅NP = 6𝑅c and 158 ligands for 𝑅NP = 8𝑅c . The polymer-polymer (P-P), NP-NP (N-N, N-C and C-C), ligand-NP (L-N) and ligandligand (L-L) bonds have a force constant 𝐾𝐵 = 120 𝑘B 𝑇/𝑅c2 and an equilibrium bond length 𝑟e = 0.8𝑅c (0.57 nm). All parameters are listed in Table 1. The boxes for simulations of 8 𝑅c , 12 𝑅c and 16 𝑅c NPs have dimensions 40×40×50 𝑅c3 , 40×40×52 𝑅c3 and 40×40×56𝑅c3 respectively. The extension in Z dimension provides a sufficient space for NP and PB to secure the uniformity of the bulk solvent outside PB. Table 1.Parameters for short-range repulsion and harmonic bond potentials. The repulsion parameter 𝑎𝑖𝑗 (in units 𝑘𝐵 𝑇/𝑅𝑐 ) between the beads and the bond parameters, the force constants KB and the equilibrium bond length re used in the simulations. Bead types are denoted as: C-core beads; N- NP surface beads; L- ligand beads; A- ghost particle beads; P-polymer beads; G-good solvent; B-bad solvent1, B2-bad solvent 2, S-substrate.

Repulsion parameters aIJ, kBT/Rc I\J C N L A P G B/B2 S

C 42.0

N 60.0 42.0

L 60. 42.0 42.0

A 0.0 0.0 0.0 0.0

P 60.0 38.5 42.0 0.0 42.0

G 60.0 43.5 43.5 0.0 43.5 42.0

B/B2 60.0 43.5 49.0/46.0 0.0 49.0/46.0 43.5 42.0

S 60.0 42.0 42.0 0.0 42.0 42.0 42.0 42.0

Harmonic bond parameters Type P-P L-L

KB, (kBT/Rc2) re/Rc 120 0.8 120 0.8

Type C-C GT

KB, (kBT/Rc2) 120 0.01

re/Rc 0.8 0

To analyze the dependence of the NP-PB adhesion on the NP size and solvent composition, we explore three characteristic systems differed by PB density and solvent-polymer interaction. In Systems 1, the PB consists of 484 polymer chains uniformly grafted on the substrate with the 10 ACS Paragon Plus Environment

Page 11 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 1. (a) Schematics of free energy calculations by the ghost tweezers method. Polymers are attached to a substrate to form a brush immersed in a binary solvent. NP is tied to its immobile ‘twin ghost tweezer’ particle via harmonic bonds. The NP-GT displacment Δz gives the force excerted by PB on NP. (b) Representative snapshot of PB-NP system used for free energy simulations. Colors: white beads- polymer, yellow-substrate, blue-ligands and green-NP. The ghost particle beads are shown as transperent red beads. Solvent beads are not shown for clarity. (b) Representative snapshot of NP in a polymer-grafted channel under flow. The solvent velocity v(z) is approximated by a parabolic profile indicated by the arrows..

grafting density Γ = 0.6 nm-2. In Systems 2, the grafting density is reduced to Γ = 0.36 nm-2 which amounts to 289 polymers uniformly grafted on the substrate. It is worth noting that we considered a limited range of grafting densities (0.36-0.6 nn-2) at which PB is homogeneous. At lower grafting densities, one may expect additional effects related to the inhomogeneous distribution of grafted chains as the solvent quality worsens. These effects deserve a special investigation. In system 3, we improve quality of the bad solvent to make bad solvent-polymer interactions more favorable; the grafting density is the same as in Systems 1, Γ= 0.6 nm-2. The bad solvent in Systems 3 is denoted as B2. In these systems, we obtain adhesion free energy landscapes for NPs of three sizes at five (Systems 1 and 2) and six (Systems 3) different solvent compositions. For each system, we calculate the free energy landscape using the GT method


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

explained below. Calculation of the free energy of NP at given distance z from the substrate requires an independent simulation of 1 million time steps or more. All together, we report the results of more than 590 independent simulations of at least 1 million time steps. 2.4 Calculation of NP-PB interaction free energy by the ghost tweezers method. To measure the free energy of interactions between NP and PB, we apply the ghost tweezers (GT) method53, which in-silico mimics lab experiments with optical tweezers.59 The GT method is implemented by introducing a ghost particle that is an “identical twin” of the real NP (Fig. 1a). The ghost twin (GT) particle is composed of beads of type A and of the same number of beads as the real NP kept in the undisturbed hexagonal order. GT is pinned at given position 𝐑 GT ; it is immobile and does not interact with any system components except for the NP. The NP is bound to GT by M harmonic springs of rigidity kGT connecting each respective pair of real and GT particle beads (i=1,..,M). Effectively, the NP is tethered to the given point 𝐑 𝐺𝑇 by an effective force exerted by the GT, 𝐅GT = ∑𝑀 𝑖=0 𝑘GT (𝐑 𝑖 − 𝐑 GT,𝑖 ), where 𝐑 𝑖 is the fluctuating position of the i-th bead of NP and 𝐑 𝐺𝑇,𝑖 is the fixed position of its counterpart bead of the GT particle. This force counter-balances the sum of the forces between NP and polymer and solvent beads at given NP position. Provided the system equilibration in the course of DPD simulation with the GT placed at 𝐑 GT , we calculate the average distance 𝑧̅̅ between the NP center of mass and substrate surface and the average GT force, 𝐹GT =< 𝐅GT >, which due to the lateral uniformity of the system may have only one non-zero normal component, which depends of the deviation, z = 𝑧̅̅ − 𝑧GT , of the NP position 𝑧̅̅ from the GT position 𝑧GT , 𝐹GT = 𝐾GT (𝑧̅̅ − 𝑧GT ). Here, 𝐾GT = 𝑀𝑘GT is the cumulative spring constant of the effective harmonic force acting between the real and ghost NPs. 𝐾GT is chosen to optimize the efficiency of calculations. At the beginning of a simulation the GT is pinned to a certain location z0 far enough from PB to ensure negligible NP12 ACS Paragon Plus Environment

Page 12 of 47

Page 13 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

PB interaction so that at equilibrium NP fluctuates around z0,z = 0 within the accuracy of the simulation and respectively, the measured GT force, 𝐹GT = 0. Then the GT is moved along the Z-axis in small increments toward the substrate. At each GT position 𝑧𝐺𝑇 , the system is equilibrated allowing for the effective force acting between the NP and PB to be counterbalanced by the GT force, which is determined as a function of the NP position, 𝐹GT (𝑧), from the average deviation of the NP center of mass from the tether point as described above. 𝐹GT (𝑧) > 0 corresponds to an effective repulsion between PB and NP, while 𝐹GT (𝑧) < 0 corresponds to effective attraction. At the equilibrium positions, 𝐹GT (𝑧) = 0. With the simulation length we afforded in this work, an average error for the force data points is about of 0.5 kT/Rc that is sufficient to determine the free energy landscapes that display different qualitative behavior with reasonable accuracy. The change in the Helmholtz free energy, 𝐴(𝑧) is determined as the work done by the GT in the course of quasi-equilibrium pulling of NP from the bulk solvent (position z0) to the position z at the substrate, 𝑧

𝐴(𝑧) = ∫ 𝐹GT (𝑧)𝑑𝑧

(3)

𝑧0

The NP location z0 far from the substrate serves as a reference point where A = 0. The free energy landscape, 𝐴(𝑧) allows one to estimate the equilibrium distribution of NP near the substrate, determine the excess Gibbs adsorption of NPs and calculate the Henry constant that quantifies the adhesion equilibrium. Let us consider a slit-like pore of half-width w with walls grafted with PB. The probability density 𝑃(𝑧) of the NP location can be calculated as

𝑃(𝑧) = 𝑒

−

2𝑤 𝐴(𝑧) 𝐴(𝑧) − 𝑘B 𝑇 ⁄ ∫ 𝑒 𝑘B 𝑇

𝑑𝑧 .

0


(4)

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 47

The partition coefficient 𝐾of the nanoparticle between the polymer brush and the bulk solvent can be obtained using the Gibbs excess adsorption theory and calculating the Henry constant, as detailed in the previous works.29-30, 60 The Henry constant 𝐾𝐻 is defined from the excess adsorption of NPs per unit pore surface area, which represents the difference between the total amount of NPs in the channel 𝑁tot and the amount of NPs in a reference volume of bulk fluid with given concentration of NPs denoted as 𝑐B . This concentration is considered as low so that interactions between NP are ignored. Since our ultimate goal is to understand the specifics of NP separation in pore channels with PB-grafted pore walls, it is reasonable to define the reference bulk volume as the volume available for the solvent flow. As shown below, the solvent flow does not penetrate into the PB and the velocity profile is approximately parabolic (Fig.1c) that would be observed in a channel with solid walls of half-width 𝑤H = 𝑤 − 𝑤PB . We define 𝑤H as the hydrodynamic width of the PB grafted channel and 𝑤PB as the hydrodynamic thickness of the PB. The reference bulk volume for the definition of excess adsorption Nex represents the volume of the mobile phase outside of the PB

𝑁ex =

𝑤 𝐴(𝑧) 1 − 𝑁tot − 𝑐B 𝑤H = cB ∫ 𝑒 𝑘B 𝑇 𝑑𝑧 − 𝑐B 𝑤H 2 0

(5)

Here Nex is the total amount of NPs per unit substrate surface area. Accordingly, the Henry constant is given by

𝐾H =

𝑁ex 𝑐B

𝑤

= ∫0 𝑒

−

𝐴(𝑧) 𝑘B 𝑇

𝑤

𝑑𝑧 − (𝑤 − 𝑤PB ) = 𝑤PB + ∫0 (𝑒

−

𝐴(𝑧) 𝑘B 𝑇

− 1) 𝑑𝑧.

(6)

As the reference technique, we consider the hydrodynamic chromatography employed for size exclusion separation of polymers and NPs alike17, 61-62. In HDC, the separation is achieved by the particle size due to the inhomogeneity of the solvent flow in the pore channels with solid walls. It 14 ACS Paragon Plus Environment

Page 15 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

is worth noting that in this case the Henry constant defined by Eq. 6 in the absence of PB (wPB=0) is negative and equal to the NP radius 𝑅NP , 𝐾H = −𝑅NP, that reflects NP depletion (negative excess adsorption) at the solid surface (A(z) = ∞ at z < 𝑅NP ), that is characteristic for HDC. 2.4 Modeling NP transport In order to model NP transport through a chromatographic channel, we simulate the solvent flow through a slit-like channel with PB-grafted walls. The system is periodic in the lateral directions. We consider 16 systems, which differ by the solvent composition and PB grafting density corresponding to Systems 1, 2, and 3. The simulation boxes in solvent flow modeling have the size of 40×40×70 𝑅𝑐3 (28.4×28.4×49.7 nm3). Selected simulations of NP flow are performed in larger simulation box of size 40×40×90 𝑅𝑐3 (28.4×28.4×63.9 nm3). We first equilibrate the systems at given solvent composition, then apply a constant force in X-direction on the solvent particles and run the simulation until the steady flow is reached. The velocity profile 𝑣s (𝑧) of the solvent across the channel is obtained by averaging the solvent particle velocities over a long time (100000 steps). For all systems considered, the velocity of solvent particles inside PB is zero while in the channel outside PB, the velocity profile tends to be parabolic similar to the Poiseuille flow in the channel with solid walls, see Fig. 1c. This observation allows us to introduce the hydrodynamic thickness wPB of PB as the cut-off of the parabolic approximation assuming that solvent velocity is zero at z= 0 ⁄𝑤 ∫0 𝜌s (𝑧) 𝑑𝑧

(7)

Eq. 7 accounts for an inhomogeneous distribution of the solvent in the channel due to its depletion within PB and the fact that solvent located within PB is effectively immobile at z . The normalized velocity profile 𝑣n (𝑧) = 𝑣s (𝑧)/ < 𝑣s > is found to be independent of the applied force. In the parabolic approximation, the mean velocity is given by, 𝑤

< 𝑣𝑠𝑃𝐴 > =

𝑣𝑚 ∫0 𝜌s (𝑧) (1 −

(𝑤 − 𝑧)2 ) 𝑑𝑧 (𝑤 − 𝑤PB )2 ⁄𝑤 = 𝑣m 𝜒 ∫0 𝜌s (𝑧) 𝑑𝑧

(8)

𝜒 =< 𝑣sPA >/𝑣m , the mean solvent velocity normalized by the maximum solvent velocity, defined by Eq. 8. The normalized velocity profile in the parabolic approximation is then

𝑣nPA (𝑧) = 𝜒 −1 (1 −

(𝑤 − 𝑧)2 ) 𝑎𝑡 𝑤PB ≤ 𝑧 ≤ 𝑤 (9) (𝑤 − 𝑤PB )2 16

ACS Paragon Plus Environment

Page 17 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

and zero for other values of z. The solvent density is constant outside the brush, 𝜌𝑠 (𝑧) = 𝜌b ( 𝜌b = 3 𝑅𝑐−3 ) while inside the brush, 𝜌𝑠 (𝑧) = 𝜌b −𝜌PB (𝑧), where 𝜌PB (𝑧) is the PB density. In estimating the mean NP velocity, we followed the assumption61, 63-64 accepted in the conventional HDC with solid wall channels that the velocity of particles positioned at distance z from the wall equals the solvent velocity at this distance. The mean velocity of NPs is determined by the distribution of solvent velocity across the channel. The fluid velocity in HDC, 𝑣nHDC (𝑧) is assumed to be Poiseuillean and is given by Eq. 9 with 𝑤PB = 0 and 𝜒 =2/3. Provided that the particle center cannot be closer to the wall than its radius 𝑅NP , the mean particle velocity reduced to the mean solvent velocity equals

< 𝑣pHDC >=

𝑤 1 1 ∫ 𝑣𝑛HDC (𝑧)𝑑𝑧 = 1 + 𝛿 − 𝛿 2 (10) 𝑤 − 𝑅NP 𝑅NP 3

where  = 𝑅NP /𝑤. The reciprocal of < 𝑣pHDC > represents the ratio of the particle retention time 𝜏pHDC to the solvent retention time 𝜏𝑠 , 𝜏pHDC 1 = ≈ 1 − 𝛿. 1 𝜏s 1 + 𝛿 − 3 𝛿2

(11)

The latter approximate equality holds at 𝛿 ≪ 1. In the case of NP flow in a PB-grafted channel, we assume that NP velocity is equal to the solvent velocity given by Eq. 9 at 𝑤PB + 𝑅𝑁𝑃 ≤ 𝑧 ≤ 𝑤 and equals 0 at 𝑧 < 𝑤PB + 𝑅NP , i.e. when NP is partially or completely immersed into PB, it is assumed to be immobile. Since the probability of NP location at distance z is proportional to the Boltzmann factor 𝑒 −𝐴(𝑧)/𝑘B 𝑇 , the mean NP velocity (reduced to the mean solvent velocity) is calculated as 17 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

< 𝑣𝑝 >=

𝐴(𝑧) − 𝑤 ∫𝑤 +𝑅 𝑒 𝑘B 𝑇 𝑣n (𝑧)𝑑𝑧 PB 𝑁𝑃 ⁄

Page 18 of 47

𝑤

∫0 𝑒 −𝐴(𝑧)/𝑘B 𝑇 𝑑𝑧

. (12)

The ratio of the particle retention time 𝜏𝑝 to the solvent retention time 𝜏𝑠 is respectively calculated as 𝜏p ⁄𝜏s = 1⁄< 𝑣 > . (13) p In the case of poor adsorption of NPs to PB, one may expect that Eq. 13 converts to the HDC Eq. 11 with 𝛿 = 𝑅NP /(𝑤 − 𝑤PB ). The parabolic approximation (Eq. 9) with the hydrodynamic thickness 𝑤PB of PB determined from the solvent flow simulation in one channel of particular width allows one to estimate the mean NP velocity and respectively, the retention time for the channels of any width. This is particularly useful considering the fact that simulations of NP flow through channels suffers from serious computational difficulties. Firstly, the real flow velocities in chromatographic columns are very low compared to the thermal velocities of particles in DPD simulations. Thus, the force applied to create such a flow is extremely weak, and it is hardly possible to obtain necessary statistics on NP diffusion across the channel cross-section even in a reasonably long simulation. Secondly, simulations of large channels are computationally unfeasible. In fact, we have shown with selected simulations of NP flow that the assumption that the NP velocity can be approximated by the solvent velocity in the absence of NP does not cause any significant error (see Supporting Information (SI), Section 1). It is worth noting that the parabolic approximation for NP velocity is a simplified model; it does not consider the effects of possible partial NP penetration into PB, which would introduce additional frictional forces to NP flow. It also ignores the NP rotational motion that is pronounced especially at the edge of PB. 18 ACS Paragon Plus Environment

Page 19 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 2. PB density (monomers per unit volume) profiles at different good solvent fractions in the simulations. Red symbols corresponds to system 1 simulations with polymer-grafting density 0.6 nm-2 at xG=0.99 (circles), 0.90 (asterisks) and 0.80 (squares). Green symbols represents densities in system 2 with grafting density 0.36 nm-2 where circles correspond to xG=0.99, asterisks to 0.89 and squares to 0.80. Blue symbols corresponds to system 3 with bad solvent 2 and grafting density 0.6 nm-2. Here densities at xG=0.99 (circles), 0.72 (asterisks) and 0.60 (squares) are shown.

3. Results and Discussion 3.1 Free energy of NP adhesion to the polymer brush We study the dependence of the adhesion free energy on the solvent quality by varying the good solvent fraction xG from 0.99 to lower values up to 0.6. At xG = 0.99, PB is fully expanded; as solvent quality worsens PB gradually transforms to a collapsed state. This is shown in Fig. 2 that displays density profiles obtained in different simulation systems: as xG decreases, the brush becomes denser and its effective width decreases. For example, in Systems 1, PB at xG = 0.8 is effectively about 15 % thinner and respectively, denser than PB at xG = 0.99. Solvent composition strongly affects interaction between the brush and NP. Figure 3 displays the dependence of the effective NP-PB force and the free energy landscape as functions of the NP coordinate (distance of NP center from substrate surface) for Systems 1 (G-B solvent mixture, Γ


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. The force (upper panel) and free energy (lower panel) profiles (right y-axis) in system 1 at different good solvent fractions xG ranging from 0.80 to 0.99 for different NP diameters 8Rc (blue), 12Rc (green) and 16 Rc (red). z- coordinate corresponds to the position of the NP center of mass. Points show the actual values of the effective force obtained with the GT method while the curves display the smoothed data (moving average). Free energy curves are obtained by integrating the smoothed data. Note the different scales of the right y-axis chosen to show the magnitude of the respective values of force and energy. Regions of negative force correspond to the NP-PB attraction, and zero force values indicate free energy minima correspond to the equilibrium adhesion positions. Polymer grafting density is 0.6 nm-2. Black squares show the local density profiles of PB beads (left y-axis).

= 0.6nm-2). At xG = 0.99, we observe no adhesion of NP disregarding of its size: the effective force is zero in the channel center beyond PB and becomes repulsive as NP is getting closer to PB (Fig. 4a). Repulsion raises steeply as NP is immersed into PB as it was observed in previously published works.50, 53 As solvent quality is reduced (xG decreases), an enthalpic attraction between NP and PB emerges, since interaction of NP surface (or ligands) with polymer is more favorable than with the bad solvent component. In the region of attraction, force is negative – GT prevents NP from penetration into PB. The position of zero force corresponds to the adhesion equilibrium: at smaller distances entropic repulsion overpowers the attractive forces, which prevents the NP from further immersion into the PB. Already at 𝑥G = 0.92 , Systems 1 exhibits slightly attractive force (negative values) at PB surface but repulsion increases steeply as NP is submerged into PB. Attraction becomes progressively more pronounced as the bad solvent fraction increases (xG decreases) and PB contracts. The location 20 ACS Paragon Plus Environment

Page 20 of 47

Page 21 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

of the minimum of the force is also shifted closer to substrate which results not only from the enthalpic attraction itself, but also from the decrease in the brush thickness (Fig. 2).

Figure 4. Characteristic equilibrium positions of NPs in size exclusion (upper panels) and adsorption (lower panels) regimes: (a) Systems 1, (b) Systems 2 and (c) Systems 3. The NP diameter is 16 𝑅𝑐 . The colors: yellow-substrate, pink-polymer, green-NP and blue-ligands.

The lower panels in Fig. 3 display the dependence of the Helmholtz free energy calculated by integrating the measured force using Eq. 3, on the position of the NP center of mass. For integration, simulated force data points are smoothened using moving averages (solid lines of the force profiles). At 𝑥G = 0.99, free energy landscapes are monotonic and do not indicate any particular equilibrium state (since NP-PB interactions are repulsive). This behavior is reminiscent to that in HDC. As the bad solvent fraction increases, free energy landscapes exhibit well-defined minima that corresponds to the adhesion equilibrium with NP partially immersed into PB. The corresponding snapshots are shown in Fig. 4a. The location of equilibrium and the corresponding equilibrium adhesion energy depends on the NP size in a nontrivial fashion. When the good solvent fraction is high, that is, when PB is fully expanded, repulsion is stronger for larger NP: the larger the NP, the more pronounced is the restriction on 21 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. The force (upper panel) and free energy (lower panel) profiles of the NPs in Systems 2 at different good solvent fractions. z is distance from the substrate surface expresed in Rc.The polymer grafting density of the brush is 0.36 nm-2 Colors: black – brush density. Blue, green, and red are data with NPs of diameters 8R c, 12Rc and 16 Rc respectively. The representation scheme is the same as in Fig.3.

the conformational freedom of the polymer chains. When the good solvent fraction is low and PB is denser, the adsorption interactions prevail and larger NPs have larger adsorption energy due to larger NP surface area on which polymer segments can get adsorbed. Thus, there is a cross-over, upon the solvent quality variation, in the adsorption free energy dependence on the NP size, which is demonstrated by the free energy profiles (Fig. 3). At sufficiently high good solvent fraction, attraction is more pronounced for smaller NPs: for example, at 𝑥G = 0.92 and 𝑥G = 0.90, smaller particles have wider regions of attraction where adsorption interactions prevail due to increased NP penetration into PB because of their smaller size. As xG decreases, (𝑥G = 0.87 and 𝑥G = 0.80) the adsorption energy increases with the NP size. The intermediate solvent compositions of xG= 0.9-0.92 correspond to a region of transition from the entropydominated to the enthalpy-dominated NP-PB interaction, or from the size exclusion mode to the adsorption mode of NP separation. As discussed in the introduction, a similar behavior is typical for adsorption of polymers, where there exists a compensation point, known as the critical point of adsorption (CPA), where the enthalpic attraction is counter-balanced by the entropic repulsion 22 ACS Paragon Plus Environment

Page 22 of 47

Page 23 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

disregarding of the polymer size. For polymers, CPA is characterized by the size independence of the Henry constant and retention time. We could not identify such a point for NP adhesion, however the observed transition from repulsive to adsorption regimes occurs in a narrow range of solvent compositions for all the systems studied. Figure 5 shows the effective forces and free energy landscapes for Systems 2, which is similar to Systems 1 except for a lower grafting density (Γ = 0.36 nm-2). The PB in Systems 2 is effectively about 4Rc (~21%) thinner than in Systems 1 (Fig. 2). Reduced PB density leads to less efficient entropic repulsion and therefore the interval of z values corresponding to attractive NP-PB interaction widens. This causes significantly stronger adhesion: free energy minima at lower xG in Systems 2 are much lower than those in Systems 1 at the same compositions. However, the transition from size-exclusion to adsorption modes occurs in the same range of compositions as in Systems 1, at xG = 0.92 - 0.89. Essentially, the cross-over behavior is similar to Systems 1. Fig. 4b shows snapshots of the equilibrium positions of the NP of diameter 6 Rc at different solvent compositions. In Systems 3, the polymer grafting density is the same as in Systems 1 ( 𝛤 = 0.6 nm-2) but the bad component interacts with polymers more favorably (𝑎P,B2 = 46.0 𝑘B 𝑇/𝑅c instead of 49.0 𝑘B 𝑇/𝑅c in Systems 1 and 2). The resulting force and free energy profiles in Systems 3 are shown in Fig.6. The snapshots of NP equilibrium states in size exclusion and adsorption modes are given in Fig. 4c. A better interaction between bad solvent and polymer radically weakens the enthalpic attraction between NP and PB because the latter results from the contrast between NPPB and solvent-PB interactions. In Systems 3, we observe the same cross-over between the size exclusion mode, where repulsion is stronger for larger NPs, and the adsorption mode with the opposite trend of stronger attraction for larger NPs. However, due to generally weaker attractive 23 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 47

Figure 6. The force (upper panel) free energy (lower panel) profiles of NPs in Systems 3. Polymer grafting density is 0.6 nm-2, same as in Systems 1 but the bad solvent component is B2 that interacts more favourably with the polymer than the one in Systems 1. Color and represeantion schemes are the same as that of Fig.3.

forces in Systems 3, the cross-over region shifts to lower xG values compared to Systems 1 and 2: in Systems 3 it occurs at xG ≈ 0.72 (Fig. 6). Overall, we conclude that the transition between the size exclusion and adsorption modes is observed for all systems studied, and the conditions at which the cross-over occurs depend on the PB density and specifics of solvent-polymer interactions. Adhesion of NP to PB results from the preferential interaction of the polymer with NP over the bad solvent. Bad solvent is excluded not only from the exterior of PB, but also from the gap between PB and NP that causes deeper penetration of NP into PB and increased effective NP-PB attraction. The strength of NP-PB attraction therefore depends on the strength of repulsion between polymer and bad solvent. Fig.7a and 7b show snapshots of NP-PB configurations at xG=0.99 and xG=0.80 for PB of density 0.6 nm-2. Visual comparison of these two snapshots shows the specifics of NP-PB adhesion in binary solvents. In the absence of bad solvent (Fig. 7a), PB prefers to be surrounded by good solvent effectively repelling NP. Addition


Page 25 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 7. (a) Snapshot of NP (RNP = 6Rc) at the closest proximity (z = WPB + RNP) to PB in the size-exclusion regime, at xG =0.99. Bad solvent fraction is negligible and NP is repelled from the brush. (b) Snapshot of adsorbed NP on PB of density 0.6 nm-2at xG=0.80. Bad solvent beads are shown in yellow, NP in green and ligands in blue. Polymer beads are shown as transparent pink beads. Bad solvent is excluded from PB causing its contraction and enhancing NP-PB interaction. (c) Solvent (red-bad solvent, green-good solvent) and PB (black) density profiles at xG=0.80 (squares) and at xG=0.99 (circles)

of bad solvent (Fig.7b), causes PB contraction and solvent exclusion between NP and PB that enhances NP penetration into PB. The differences in the polymer and solvent density profiles shown in Fig.7c quantify the change in the composition of the PB interface upon addition of bad solvent. 3.2 Solvent flow in polymer-grafted channel Chromatographic separation of NPs depends not only on the adhesion interactions but also on the specifics of the solvent flow through polymer-grafted channels that are also controlled by the solvent quality. Fig. 8 presents the results of simulation of solvent flow through slit–like channels between PB-grafted walls as shown in Fig.1c. We consider the same PB-binary solvent Systems 1, 2, and 3 as in the free energy calculations. The simulation boxes are open and periodic in the XY plane. A constant force of 0.002 kBT/Rc is applied to each solvent particle in X-direction and the simulation for several hundred thousand steps is performed until a steady flow is reached with minor fluctuations of velocities around the averaged velocity profile vs(z). Fig. 8a-8c show the profiles of the solvent velocity along with PB and solvent densities across the channel cross-section at different solvent composition. The computed solvent velocity vs(z) 25 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8. The solvent velocity and density profiles across PB-grafted channel of width 66.75 Rc at different good solvent fractions xG in Systems 1 (a), 2 (b), and 3 (c) . Computed solvent velocity vs(z) (dashed lines), normalized solvent velocity vn(z) (solid lines), PB density (circles) and binary solvent density (crosses) are shown. Red, green, blue, cyan and magenta colors represent, respectively, xG= 0.80, 0.87, 0.90, 0.92 and 0.99 in (a) and (b) and 0.60, 0.70, 0.72,0.80 and 0.99 in (c).

shown by dashed lines depends on the magnitude of the applied force. The normalized velocity, vn(z), reduced to the mean solvent velocity (shown by the solid lines) does not depend on the magnitude of the applied force (Eq. 7), as shown in the SI (Section 2). The use of the normalized velocities is practical for comparison of the role of solvent composition on the fluid flow in the chromatographic columns at a given volumetric flow rate. The normalization of vn(z) implies that the mean reduced velocity =1 ensuring the same volumetric flux at different solvent compositions. It is worth noting that PB conformation at given solvent composition and respective PB density profile are not affected by the flow (see SI Section 2); this feature was also shown experimentally.65-66 At the same time, the fluid velocity is sensitive to the change of the solvent composition because the latter affects the brush extension and respectively, the hydrodynamic width of the channel. Fig. 8a presents the solvent velocity profiles at good solvent fractions 0.99, 0.92, 0.90, 0.87 and 0.80 in Systems 1 with PB density 𝛤 = 0.6 nm-2. The channel width, which is the distance between the substrate surfaces in Z-direction is 66.75 𝑅c (47.4 nm). The profiles show that the solvent velocity is zero inside PB, but within the space between the brushes the solvent 26 ACS Paragon Plus Environment

Page 26 of 47

Page 27 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

velocity profiles tend to be parabolic as in a Poiseuille’s flow between solid walls. The density distributions of the solvent and the polymer are correlated with the velocities in the outer region of the brush, where the velocity profile changes from zero to the parabolic type. As shown in Fig.8, the solvent density ρs(z)=ρb= 3Rc-3 in the bulk region (outside the brush) while inside the brush ρs(z)= ρb  ρPB(z). Thus, there is a considerable amount of solvent inside PB and this solvent is practically immobile. In HDC channels without PB, Eq. 8 integrates to 2/3 vm (see Section 2.4) and thus the normalised velocities vn(z) have a peak of 3/2. The presence of PB decreases the mean velocity because a part of the solvent that is inside PB has zero velocity. Integrating Eq. 8, taking into account that the solvent density is constant in the bulk (z>wPB) and is approximately given by ρb

 in the interior of the brush, we obtain

≈3 (𝜌b −< 𝜌PB >)𝑤PB + (𝑤 − 𝑤PB )𝜌b , (14)

from which we get back to Eq. 8 when, = ρb (means no penetration of the solvent into the brush) and also when wPB=0. When approaches the limiting value 2/3 vm(wwPB)/w. Thus, as the mean velocity is smaller in the presence of PB, vn(z) becomes larger than 3/2. The mean solvent velocity of the PB-grated channel is determined approximately by the relative hydrodynamic width ((wwPB)/w = wH/w) of the channel that determines the fraction of the solvent that flows. In Fig. 8a, the normalized velocities become as high as close to 3 for the narrowest channel. Fig. 8b presents the density and velocity distributions in a channel with PB of density 0.36 nm-2 for different solvent compositions. The hydrodynamic width of the channel in


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 47

this case is increased for all xG because of the low grafting density. The magnitude of the solvent velocities also increases because of the lower PB density causing less resistance to the flow and the increased volume fraction of the solvent. On the other hand, vn(z) is lowered because of the larger hydrodynamic width of the channel. The use of a less repelling bad solvent B2 also causes changes in the velocity distributions (Fig. 8c). In this case, Systems 1 and 3 at xG = 0.99 are nearly the same, since the bad solvent fraction is negligible and as a result, the velocity and density profiles are similar. However, as xG decreases, the change of the PB conformation and the flow velocity is more moderate than in Systems 1. The hydrodynamic width of PB, 𝑤𝑃𝐵 is found by fitting vs(z) obtained from the simulations with parabolic approximation vsPA(z) at a given mean velocity, 𝑤

∫

𝑤

𝑣𝑠PA (𝑧) 𝑑𝑧

= ∫ 𝑣𝑠 (𝑧) 𝑑𝑧 .

𝑤PB

(15)

0

Table 2. The hydrodynamic PB thickness (wPB) at different good solvent fractions for Systems 1, 2 and 3

xG 0.99 0.92 0.90 0.89 0.87 0.80 0.72 0.70 0.65 0.60

Systems 1 21.1 19.8 19.2 18.8 17.5

WPB (RC) Systems 2 16.5 15.0 14.4 14.2 12.9

Systems 3 21.2

20.1 19.5 19.4 19.0 18.3

wPB can be iteratively found by solving Eq. 14 (see section 3 in the SI for details). The values of wPB obtained for various systems are shown in Table 2. 28 ACS Paragon Plus Environment

Page 29 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The hydrodynamic thickness wPB is the appropriate definition of the effective PB thickness when our purpose is to describe the flow through PB-grafted channels. It is worth noting that the

Figure 9. (a) Comparison of parabolic approximation (solid lines) with simulation results (circles, squares). Red, green and blue represent solvent velocities vs(z) for xG=80 in Systems 1, 2 and 3 respectively. Black, cyan and magenta are the corresponding normalized profiles vn(z). (b) Normalized velocity profiles (solid lines) at xG=0.99 in Systems 1 in channels of different widths w, obtained using the parabolic approximation. Black, w = 33.4Rc; red, w = 50Rc; green, w = 100Rc; blue, w = 200Rc; magenta, w = 500Rc. The hydrodynamic thickness 21.12 Rc is used. The solvent densities (dashes lines) which are taken to be constant in the bulk, are also shown.

hydrodynamic PB thickness is different from the effective geometrical thickness wG defined by approximating the PB density profile with the rectangular one, assuming a uniform equilibrium PB density at z < wG, see Section 3 of SI for details. Fig. 9a compares vs(z) and vn(z) obtained from the simulations with the corresponding parabolic approximations at xG=0.80. Both solvent (circles) and normalized velocities (squares) showed


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 47

Figure 10. Solvent composition dependence of the Henry constant in (a) Systems 1, (b) Systems 2 and (c) Systems 3. Y-axis represents the decimal logarithm of KH shifted by 10 nm to avoid negative values under logarithm in the size exclusion regime. Insets present KH - xG dependence to show the transition from the size exclusion regime with negative KH to the adsorption regime with positive KH.

very good agreement with their respective parabolic approximations for all systems (solid lines). Fig. 9b presents the predicted velocity distributions using parabolic approximations for channels of different widths ranging from the width of the simulated channel w0=33.4 Rc to 500 Rc containing PB (Systems 1) at xG=0.99. The solvent density ρs(z) in the channels wider than w0 was assumed equal to the computed solvent density in the channel of width w0 at z < w0 and set equal to the bulk density ρb at w0 < z < w. The solvent densities in Fig. 9b are shown by dashed lines. The normalized velocities vnPA(z) are shown in Fig. 9b by solid lines. Since they are constructed to provide the reduced flux equaled 1, the maximum normalized velocity in the channel center gets smaller as the channel size increases. At w =500 Rc, the maximum value of vsPA(z) =1.54, close to the value of 3/2 that would be observed in the channels with solid walls without PB. 3.3 Henry Constant of NP Adhesion We quantify NP-PB adhesion by the Henry constant calculated according to Eq. 6 from the free energy landscapes determined in the simulations with the GT method. The Henry constant is calculated by Eq. 6 as the excess of adsorption per unit substrate surface area using the reference 30 ACS Paragon Plus Environment

Page 31 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

bulk volume based on the PB hydrodynamic thickness 𝑤PB defined from the Poiseuille approximation of the solvent flow at given solvent composition. The calculated Henry constant 𝐾H as a function of the good solvent fraction xG is presented in Fig. 10. In the repulsive size exclusion regime, when the good solvent fraction is sufficiently high, the free energy is either positive or close to zero and the integrand in Eq. (10), 𝑒 −

𝐴(𝑧) 𝑘𝑇

≤ 1.

As a result, 𝐾H is negative. In a limiting case of very dense brush and purely repulsive interactions, 𝐾H ≈ −𝑅NP , as in the HDC case, since the distances z< wPB + RNP from the substrate are effectively inaccessible to NP. In HDC, 61 this size exclusion leads to shorter elution times for larger NPs: they are located closer to the channel center and move faster within the Poiseuille flow. Thus, we call the repulsive region of NP-PB interaction at high good solvent fractions as the hydrodynamic, or size exclusion mode where 𝐾H is negative and decreases as the NP size increases. As the NP free energy profiles at lower good solvent fractions have regions of attraction with negative energy values, the exponential factor in integrand of Eq. (10) becomes large and positive and consequently, the Henry constants with the addition of bad solvent become positive and grow exponentially. This situation is referred as the adsorption regime of chromatography, in which the Henry constant 𝐾𝐻 is positive and increases with the NP size. As 𝐾𝐻 values vary with the good solvent fraction in an exponential manner, it is convenient to represent them in logarithmic scale. To avoid negative values under logarithm in the hydrodynamic regime, we shift 𝐾𝐻 by 10 nm and plot log (𝐾𝐻 + 10 𝑛𝑚) against the good solvent fraction xG in Fig. 10. The Henry constant values without such modification near the cross-over are displayed in the


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

insets to demonstrate the transition from the size exclusion regime with negative KH to the adsorption regime with positive KH.

Figure 11. The probability distributions P(z) of nanoparticles in the channel with PB grafting density 0.6 nm-2 (Systems 1) at good solvent fractions 0.99 (a),0.90 (b) and 0.80 (c). The solid red, green and blue curves respectively represent P(z) of NPs with diameter 16Rc, 12Rc and 8Rc while the dashed vertical lines correspondingly marks the mobility boundary wPB+RNP for each NP on PBs of both sides. The PB relative density ρPB’(z) (black) and the normalized velocity profiles vn(z) (magenta) are also shown.

The transition from the size exclusion regime with negative KH in good solvents to the adsorption regime with positive KH at poorer solvent is observed in all systems modelled in this work. In each system, this transition happens in a relatively narrow range of solvent compositions: decreasing xG by approximately 0.02 reverses the order, the largest NPs that are most strongly repelled becomes most strongly adsorbed. From our results, we cannot pinpoint the exact value of xG where the Henry constant becomes NP size independent (that xG would correspond to the critical point of NP adsorption at PB). Nevertheless, the narrowness of the concentration range corresponding to the cross-over in each particular system and relatively strong dependence of the location of the cross-over region on the interactions indicate that NP separation by surface chemistry can be achieved at least for select classes of NP mixtures. 3.4. Mean velocity and retention time of NPs in PB-grafted channels


Page 32 of 47

Page 33 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

To analyze NP separation in polymer-grafted channels under flow, we follow the approximate approach adopted in HDC61, 63-64 and assume that the NP velocity equals the solvent velocity at the position of the NP center of mass.17 Further, we employ the parabolic approximation of the flow velocity, introduce the hydrodynamic thickness of PB wPB and assume that NPs immersed (at least partially) in PB (z wPB + RNP) move with the solvent velocity vs(z). As such, the mean NP velocity and retention time are calculated by convolution of the solvent velocity in parabolic approximation and the NP probability distribution according to Eq.12 and 13. In the absence of PB, this approximation reduces to the conventional HDC equations 10 and 11. The probability density distribution P(z) of NP location in the channel is determined by the NP-PB free energy landscape for particular NP size and solvent composition according to Eq. 4. Fig.11a shows P(z) at xG=0.99 in Systems 1 along with the PB relative density, ρPB’(z)= ρPB(z)/ ρPBmax and normalized velocity vn(z) obtained from simulations. The assumed boundary of NP mobility, wPB+RNP, is shown as dashed lines. Since NPs are completely repelled from PB in such a good solvent, they tend to distribute uniformly in the bulk flow where the free energy is zero and the probability density is constant over z. It is clear from the plot that the NP probability density at z< wPB+RNP is nearly zero while at z> wPB+RNP it is constant. The hydrodynamic width of PB, wPB serves as a wall preventing penetration of NPs beyond z=wPB+RNP. In addition, the size exclusion constrains the larger particle towards the channel center where it moves faster with the Poiseuille-type flow while the smaller particles sample a wider space and wider range of velocities and therefore, move slower in average. This picture is similar to that in HDC. Fig.11b shows probability distributions for NP location at xG=0.90. Here,


Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 12. Transition from the size exclusion to adsorption regimes upon the decrease of solvent quality in Systems 1(I), 2( II) and 3( III). The channel half-width w= 70.43 Rc = 50 nm. (a) Mean nanoparticle velocity . (b) NP retention times reduced to solvent retention times τp/τs.. τp/τs 1 in the adsorption regime and τp/τs ≈1 in the transition regime.

adsorption causes penetration of NPs in the interior of PB beyond the mobility boundary. However, the chances of finding the NP in the solvent bulk are still nonzero as the NP adsorption is weak. The overall velocity of NP in this case is determined by the fraction of NPs retained within PB where the NP velocity is zero. At xG=0.80, as shown in Fig. 11c, the NP probability density is concentrated within PB at z< wPB+RNP. At these conditions NPs are mainly retained in PB in adsorbed state with a small probability to desorb. In this adsorption mode, the probability of desorption is determined by the depth of the adsorption well and decreases with the increase in NP size. Respectively, the mean velocity of NP decreases with their size. Fig. 12a depicts the mean NP velocity calculated using Eq. 9 as a function of the solvent composition in Systems 1(I), 2(II), and 3(III). is calculated using the parabolic approximation assuming a channel haft-width, w = 70.43Rc which is equivalent to a real channel 34 ACS Paragon Plus Environment

Page 34 of 47

Page 35 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

width of 100 nm. In each case , decreases with the bad solvent fraction as the adsorption interactions gradually strengthen and attract NPs towards the PB interior. This effectively makes the probability of NP location in the mobile bulk fluid lower and therefore restricts its mobility. In all three systems, the dependence of < vp> on the NP size undergoes a transition from size exclusion to adsorption modes. At higher xG, the mean velocities are larger for larger NPs due to the size exclusion of the NPs from the effectively repulsive PB as shown in Fig. 11. At low xG, the NP probability density is determined by the adsorption energy, the larger particle with larger free energy of adsorption has a smaller probability to be located in the mobile solvent bulk at any given moment. Thus, the mean velocities are smaller for larger particles. When the adsorption is very strong the NP velocity is nearly zero (Fig.12a I &II) but still decreases with the NP size. At the intermediate values of xG, the adsorption is weak and the size dependence of arises from both size exclusion and adsorption effects which approximately counter balance each other. Fig. 12b shows the logarithmic values of NP retention times reduced by the solvent retention times (τp/τs) which is the reciprocal of (Eq. 13). In all cases, the retention times in the adsorption regime indicate that NPs are retained longer than the solvent (τp /τs >1) while in the size exclusion regime the solvent is retained longer ( τp/τs

Adhesion and Separation of Nanoparticles on Polymer-Grafted Porous

Recommend Documents