Field-Theoretic Simulations of the Distribution of Nanorods in Diblock

Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Field-Theoretic Simulations of the Distribution of Nanorods in Diblock Copolymer Thin Films Huikuan Chao, Benjamin J. Lindsay, and Robert A. Riggleman* Chemical and Biomolecular Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104, United States ABSTRACT: Using block copolymer microphases to guide the self-assembly of nanorods in thin films can give rise to polymeric materials with unique optical, thermal, and mechanical properties beyond those found in neat block copolymers. Often the design and manufacture of these materials require exquisite control of the nanorod distribution, which is experimentally challenging due to the large parameter space spanned by this class of materials. Simulation approaches, on the other hand, can access the thermodynamics that contribute to the nanorod distribution and hence offer valuable guidance toward the design and manufacture of the materials. In this work, we employ complex Langevin field-theoretic simulations to examine the thermodynamic forces that govern the assembly of nanorods in thin films of block copolymers with a particular focus on vertically oriented cylindrical and lamellar domains. Our simulations show that the nanorod geometry, the substrate selectivity for the distinct blocks of the copolymer, and the film thickness all play important roles in engineering both the nanorod orientation and spatial distribution in diblock copolymer thin films. In addition, we employ thermodynamic integration to examine how the nanorods alter the stability of vertical and horizontal domains in thin films, where we find that the tendency of the nanorods to stabilize a vertical orientation depends on both the film thickness and the nanorod concentration.

■

INTRODUCTION Assembly of nanoparticles in inhomogeneous polymers has in recent years received increasing academic and industrial attention due to the potential for the development of new types of functional polymer nanocomposites (PNCs). The combination of polymer inhomogeneity and nanoscale features of the nanoparticles confer these types of materials with advanced mechanical, thermal, and optical properties.1−6 Block copolymers (BCPs) containing spherical nanoparticles have been some of the most extensively investigated of the PNC materials.7−9 Experimental and computational studies have shown that the spherical nanoparticles are found to distribute to defects in BCP morphologies, relieving the strain imposed by the defects.10,11 The self-assembly of systems with nonspherical particles may be different, however, because it has been shown that BCP morphologies depend on the physical and chemical properties of nanoparticles added to the polymer.12 Although the principles controlling the distribution of spherical particles have been well established, other nanoparticle shapes, such as nanorods (NRs) or nanoplates, have not been investigated as thoroughly.9 One of the key challenges in the NR assembly problem is how to simultaneously manipulate the spatial distribution and orientation of NRs within BCP microdomains, where a number of factors can complicate the assembly, including orientationdependent NR−NR interactions and NR−BCP interactions. Xu and co-workers13,14 introduced small ligands that can form hydrogen bonds with PS-b-P4VP BCP to form a BCP complex. © XXXX American Chemical Society

By adjusting the design of the BCP complex and geometry of the NRs, they discovered that the NR−NR interactions and NR− BCP interactions can be cooperatively manipulated. For example, they were able to achieve NR assemblies with end-toend NR alignment parallel to cylindrical domains. For microscopically confined BCP systems, such as BCP thin films, the NR orientation can be further manipulated to align vertically or horizontally with respect to the film surface. Experimental work from the Composto group15 synthesized PEG-decorated gold NRs. By implementing solvent annealing, they were able to manipulate the NRs within poly(methyl methacrylate) lamellae with an orientation in the plane of the lamellae. More recently, the Composto group has demonstrated that defects in block copolymers can lead to segregation of nanorods to the defect regions16 to relieve chain stretching.11,17−19 Due to the computational challenges associated with equilibrating block copolymers, and nanoparticle-loaded block copolymers in particular, computational studies on the assembly of NRs in BCP are quite limited.20,21 Efficient theoretical tools and alternative simulation approaches, such as field-theoretic simulations, could be a boon to the study of anisotropic particles in BCPs. The development of polymer field theory (PFT) over the past few decades has led to some great successes in studying the equilibrium behavior of inhomogeneous polymers,22,23 Received: August 7, 2017 Revised: November 13, 2017 Published: November 14, 2017 A

DOI: 10.1021/acs.jpcb.7b07862 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B including the prediction of the BCP phase diagram.24 In fieldtheoretic polymer models, polymers are typically modeled as Gaussian chains consisting of segments harmonically bonded together, and nonbonded interactions between segments are decoupled through a particle-to-field transformation. In the end, each chain interacts with the chemical potential field generated by the other chains, which opens the door to analytic treatment and in some limits can significantly reduce the computational expense.22,23 In particular, the mean-field approximation is often applied to further simplify the calculation, leading to the so-called self-consistent field theory (SCFT). Incorporation of NRs into polymer field theory is still an ongoing topic in the polymer field theory community. Early work by Balazs and co-workers has established the self-consistent field theory (SCFT)/density functional theory (DFT) method, where DFT approximation of the free energy of a hard sphere fluid is used to augment the SCFT free energy to capture the equilibrium behavior of PNCs.7 The SCFT/DFT framework has been extended to investigate the equilibrium configurations of NRs in BCPs.25 However, the orientation of the NRs was required to be predefined to render the calculations tractable. Sides and coworkers developed a hybrid particle/field theory (HPF), where the degrees of freedom associated with nanoparticles are left out during the particle-to-field transformation12 and retained as explicit degrees of freedom in the simulation. However, the numerical realization of HPF can be time consuming due to the requirement that converged SCFT simulations must be obtained to update the nanoparticles’ configuration. Ma et al. used HPF to study how the load and geometry of NRs impact bulk BCP morphologies in two-dimensional simulations.26 More recently, Ting et al.27 used HPF to tune the grafting density of polymer ligands on the ends of NRs relative to the faces to drive NRs into equilibrium side-by-side or end-to-end assemblies. Recently, our group has developed the polymer nanocomposite field theory (PNC-FT) method to incorporate nanoparticles into the PFT framework.28−30 In this approach, we begin with models very similar to the HPF approach, but we have shown that the nanoparticles can be directly integrated into the field theory, giving a unified theoretical treatment of both the polymer matrix and the particles. Importantly, PNC-FT is fully compatible with complex Langevin field-theoretic simulations, which can capture the fully fluctuating field theory. Furthermore, there are no formal restrictions on the shape of the nanoparticles incorporated into the PNC-FT, which enables us to study anisotropic particles, such as nanorods. In a recent study,16 we compared the results of our PNC-FT approach to experiments on PS-b-P2VP block copolymers containing P2VP-functionalized gold nanorods (P2VP-AuNRs). In a composition regime where the P2VP domains formed cylinders, it was found that when solvent annealing was used to drive the P2VP cylinders into a vertical orientation, the P2VP-AuNRs segregated to the base of the cylinders. The PNC-FT calculations suggested that this was due in part to the wetting conditions at the bottom substrate and that the segregation of the P2VP-AuNRs could be partially relieved by switching to a substrate with neutral wetting. In this study, we significantly extend the PNC-FT calculations from this prior study to understand what controls the distribution of NRs in block copolymer thin films. In addition to studying the effects of substrate wetting, we also investigate the role of film thickness, nanorod dimensions, and the phase of the block copolymer. Our results are interpreted in terms of the entropy of the block copolymer matrix chains, and generally speaking, when we observe a transition from horizontal to

vertical NR orientation, the transition is driven by the block copolymer entropy. In addition, we also investigate the role of the nanorods on the free-energy difference of the vertically oriented and horizontally oriented cylinder domains. We find that the NRs can help stabilize vertically oriented cylinder phases, but similar to neat diblocks, this depends on the film thickness.

■

METHODS Theory. The field-theoretic model of BCP/NR composites employed here closely follows the PNC-FT framework developed in our previous work.28 Briefly, the diblock copolymer is modeled as a discrete Gaussian chain, where N = NA + NB interacting segments are connected with the harmonic potential 3 |r|2 ub(r) = 2 , where b is the statistical segment length. In the 2b absence of nonbonded interactions, a single Gaussian chain with this setup would have an unperturbed radius of gyration, Rg =

(N − 1)b2 6

. The mass of a polymer monomer is distributed

over a unit Gaussian h(r), which has the benefit of giving correlations to the segments28 and regularizes the theory against UV divergences.31,32 A density distribution function with both spatial and orientational dependences, Γ(r, u), is used to distribute the density of a single nanorod about its center ⎛ (|r × u| − RP) ⎞ ⎛ (|r·u| − L P /2) ⎞ erfc⎜ ⎟ erfc⎜ ⎟ ξP ξP 4 ⎝ ⎠ ⎠ ⎝

ρo

Γ(r, u) =

(1)

where ρo is the bulk segment density and u is the unit vector pointing along the long axis of the nanorod. RP and LP are the radius and length of the nanorod, respectively, and ξ P characterizes the interfacial width, over which the nanorod density decays smoothly from the bulk value (ρo) at the center to zero far from the interface. The total microscopic densities of the polymer species and nanorods are given by nD

ρ̂K (r) =

NK

∑ ∑ h(r − ri ,j) i

j

(2)

and nNR

ρ̂NR (r, u) =

∑ Γ(r − ri, ui) i

(3)

respectively, where K represents the A and B segment densities arising from the BCPs. The interaction potential between unlike species I and J is given by χI,J UIJ = dr ρÎ (r) ρĴ (r) ρo (4)

∫

where ρ̂K(r) is the total density distribution of type K. In this study, the nanorods are regarded as chemically identical to the A blocks such that χA,B = χB,NR and χA,NR = 0. The density fluctuations away from the bulk value are penalized by a quadratic potential33 κ Uκ = 2ρo

∫

⎛ ⎞2 ⎜ dr⎜∑ ρK̂ (r) − ρo ⎟⎟ ⎝K ⎠

(5)

where the index K goes over all of the species plus an additional density due to the confining wall, and κ controls the strength of B


Article

The Journal of Physical Chemistry B the penalty for density fluctuations. We define the wall density (ρ̂W) using a masking function, which is commonly used in fieldtheoretic studies34

the type of the jth segment in the diblock chain. The segment center density operators for species A and B can be obtained as ρÃ (r) = −nD

⎛ | r | − L W ⎞⎤ ρ⎡ ̂ (r) = o ⎢1 − tanh⎜ z ρW ⎟⎥ 2 ⎢⎣ ⎝ ξW ⎠⎥⎦

∫ +w+∫ +w− e−/[w ]

N

ρ̃B (r) =

ρK̆ (r) = (h*ρK̃ )(r)

Q NR =

1 4πV

∫ dr∫ du e−w

NR (r, u)

(14)

The operator that yields the density of NRs with orientation u and their centers located at r can be shown to take the form

∫ dr w+

∫ dr qD(r, j = N )

(13)

The single NR partition function (QNR) takes the form

ρNR ̃ (r, u) = −nNR

Here, the spatial dependence of the potential fields are dropped out for the sake of simplicity. The potential fields wA and wB are given by wA(r) = h*(iw+ − w−)(r) and wB(r) = h*(iw+ + w−)(r), respectively, where * indicates a convolution. wNR is the potential field conjugated with the nanorods given by wNR(r, u) = [Γ*(iw+ − w−)](r, u). The prefactor z0 in the partition function is independent of the fields and only contains constant self-energy, factorials correcting for indistinguishability, and the normalization constants that result from the particle-to-field transition; these constants do not play a role in the calculations presented below. We note that the PNC-FT framework implemented here to build the model is not restricted to two chemical species and is able to describe multispecies systems with nanoparticles. Because the majority of the study only involves two chemical species, we use the two-species exchange model for the sake of brevity and refer readers interested in the multispecies field-theoretic models to either our previous work28 or the article by Delaney and Fredrickson.35 As a final point, we note that our primary interest in this work is to sample the fluctuating fields appearing in eq 7 without approximations, such as the mean-field approximation, which leads to the so-called self-consistent field theory (SCFT). As detailed below, our method of choice for sampling the fluctuating fields is the complex Langevin (CL) method. In the above model, the partition function of a single BCP chain, QD, can be calculated as 1 V

(12)

respectively. The total segment density is obtained by convolving the center density operator with h(r)

(8)

QD =

nD ∑ q (r, j) ewB(r)qD(r, N − j) VQ D j = N + 1 D A

(7)

∫

(11)

and

where /[w±] is the effective Hamiltonian of the model, a functional depending on the fields w+(r) and w−(r), which has the form

∫

nD A ∑ q (r, j) ewA(r)qD(r, N − j) VQ D j = 1 D

=

±

/[w±] = − nD ln Q D[wA ; wB] − nNR ln Q NR [wNR ] ρ ρo + o dr w−2 + dr w+2 − iρo χA,B χA,B + 2κ

δwA N

(6)

where LW is the wall thickness, ξW characterizes the interfacial width of the wall, and |rz| is the distance from the point r to the z = 0 plane in the model. This wall is centered at z = 0 and crosses the periodic boundary of the simulation box along the z-dimension, imposing a substrate on both the top and bottom of the thin film. With the model thus defined, we next use standard Hubbard− Stratonovich transformations22 to develop the field-theoretic form of our partition function A = z0

δ ln Q D

=

δ ln Q NR δwNR (r, u)

nNR e−wNR (r, u) 4πVQ NR

(15)

The NR density including the finite size of the rod can be obtained by convolving the center density operator with Γ(r, u) ̆ (r, u) = (Γ*ρNR ρNR ̃ )(r, u)

(16)

where we emphasize that the convolution is only performed over the spatial (r) coordinates. The total local NR density due to nanorods of all orientations is obtained by integrating over u, ρ̆NR(r) = ∫ du ρ̆NR(r, u). The local ordering of the NRs can be quantified using the nematic order tensor Q(r) =

̃ (r , u ) ∫ du(uu − 13 I)ρNR ̃ (r , u ) ∫ duρNR

(17)

where uu represents an outer product and I is the identity tensor. To quantify the level of ordering, we calculate the largest eigenvalue of the tensor and average it over the field fluctuations. The selectivity of the substrate between the A and B components is controlled by modifying the field wB(r) to wB(r) + s(r), where ̂ (r) s(r) = λW ρW,l

(18)

Here, λW controls the selectivity strength and the index l denotes that the selectivity is imposed only on the lower surface; the top surface remains neutral under all calculations. Finally, some simulations were performed with a soft interface at the top of the film to test the role of our boundary conditions. Similar to previous SCFT simulations,36 a soft interface was imposed by adding a 2-Rg-thick homopolymer layer between the top of the diblock film and the upper portion of the masking function. The homopolymer chains were composed of NC = 40 C-type segments, which had no enthalpic preference for any other component such that χA,C = χB,C = χNR,C, and our initial field

(9)

where qD(r, j) is the diblock chain propagator, which is calculated by iterating a Chapman−Kolmogorov equation

∫ dr Φ(r − r′)qD(r′, j − 1)

qD(r, j) = exp[−wK j(r)]

(10)

given the initial condition qD(r, j = 1) = exp[−wA(r)] and the normalized bond transition probability Φ(r). Here, Kj denotes C


Article

The Journal of Physical Chemistry B configurations were biased so that all of the C homopolymers were near the top surface. Hybrid Particle Field Theory Simulation. In several examples below, we are interested in the free-energy difference between two systems with a single nanorod that has a particular position rNR and orientation u. We use the standard HPF method12 for these calculations, and in this approach, the chemical potential fields w± are evaluated at their mean-field configurations. We emphasize that the underlying particle model in these calculations is identical to the pure field-theoretic model described above, only here we have retained a single explicit nanoparticle out of the particle-to-field transformation. The resultant HPF Hamiltonian has the form

−

(19)

(20)

/λ = λ /[w±] + (1 − λ)/ref [w±]

When λ = 1, the composite system is identical to the fully fluctuating field-theoretic model, whereas if λ = 0, then the composite system regresses to the reference system. Finally, the

∂ ln A ∂T z ∂ e−/HPF e−/HPF ∂z 0 + k bT 2 0 A ∂T A ∂T

∫ dr χ

A,B

(24)

thermodynamic integration is performed on

⎡ 3k bT NnD − ⎢ 2 ⎣⎢ +

(23)

α is a spring constant that controls the stiffness of the harmonic fluctuations around the mean-field values. We next constructed a composite Hamiltonian of /ref and /[w±]

F = Fref − (21)

where

1

dλ

∂/λ ∂λ

over λ

(25)

λ

is an ensemble averaging evaluated using complex

M

ρ

α 3 dx ∑ [(w+, i − w+*, i)2 + (w−, i − w−*, i)2 ] 2 i

∫ dr χ o (w−*)2 A,B

⎤ (w+*)2 ⎥ ⎥⎦ + 2κ

∂/λ ∂λ λ

∫0

∂/λ ∂λ

Langevin simulation in the composite system with the Hamiltonian, /λ . The free-energy of the reference system, Fref, can be evaluated from /ref , where the contribution from /* can be numerically evaluated from the mean-field solution w*± . The contribution from the harmonic potential on the other hand can be understood coming from uncoupled harmonic springs on the collocation grid in the simulation box fluctuating about the meanfield values

where the second form of the equation is evaluated under the mean-field approximation. The first term on the RHS of eq 21 accounts for the ideal gas contribution, whereas the second term comes from the field-dependent contributions. To evaluate the second term, we assume that both χA,B and κ scale inversely with T (e.g., χA,B ∼ 1/T, κ ∼ 1/T). The mean-field internal energy thus takes the following form U (rNR , u) =

(w+*)2

∫

where F(rNR, u) is the free energy of the system and can be evaluated under the mean-field approximation (to within a constant) as F(rNR , u) = /HPF[rNR , u ; w±*]. Here, and in the remainder of this work, the asterisk (*) superscript denotes a mean-field solution. U(rNR, u) is the internal energy of the system and can be evaluated from

= k bT 2

+ 2κ

To investigate how the orientation of NRs changes the BCP entropy, we set up HPFT simulations, where a single NR is positioned in contact with the substrate, in either a horizontal or vertical orientation, within a cylinder domain. The per-chain entropy difference from the BCP chain between the vertical and horizontal cases is then calculated according to ΔSvh = S(rNR,v, uv) − S(rNR,h, uh). Thermodynamic Integration. In some of our calculations below, we are interested in the free-energy difference between a film with vertical and horizontal domains using our fully fluctuating model. Because we are sampling the field fluctuations, we do not have direct access to the free energy of our system, and as a result, we have adopted the thermodynamic integration proposed by Lennon et al.37 for complex Langevin field-theoretic simulations. Briefly, this method first calculates the free-energy associated with the mean-field solutions w*± , where the freeenergy is simply the Hamiltonian evaluated at w*± , /*. A reference system is constructed on the basis of the mean-field result with a reference Hamiltonian α /ref = /* + dr[(w+ − w+*)2 + (w− − w−*)2 ] 2

where ρ̂NR(r; rNR, u) is the total density of the fixed NR with a configuration (rNR, u). After converging the fields to their meanfield configuration, we then have direct access to the free energy and other thermodynamic properties of the system. For example, the entropy can be estimated as the difference between the free energy and the internal energy

U (rNR , u) = k bT 2

∫ dr χ

2ρo

A,B

̂ (r; rNR , u)]w+(r) + ρNR ̂ (r; rNR , u)w−(r)} ∫ dr{i[ρo − ρNR ρo ρ + ∫ dr w−2 + χ +o 2κ ∫ dr w+2 χ

S(rNR , u)T = −F(rNR , u) + U (rNR , u)

2ρ

∫ dr χ o (w−*)2 A,B

−

A,B

̂ (r; rNR , u)]w+* ∫ dr{i[ρo − ρNR

̂ (r; rNR , u)w−*} + nD ln Q D[wA*; wB*] − + ρNR

/HPF[rNR , u ; w±] = − nD ln Q D[wA ; wB]

A,B

3k bT NnD + 2

S(rNR , u)T =

where dx is the spatial resolution in the grid and M is the total number of collocation points. The free-energy of the 2M uncoupled harmonic springs can be analytically evaluated as Gaussian integrals for each w±,i degree of freedom. We employed 20 values of λ between 0 and 1 to discretize the thermodynamic integration.

ρo

(22)

This gives rise to the equation that we used to calculate the average entropy of the diblock system D


Article

The Journal of Physical Chemistry B

Figure 1. (a) Fraction of vertically oriented nanorods f V as a function of the strength of the bottom substrate selectivity λW in films with hard confinement (solid) or soft confinement (dashed). Soft confinement was simulated with a layer of immiscible polymer C on the top surface of the film, where NC = ND and χC,K = χA,B for all K ∈ {A, B, NR}. (b) Density field maps containing isosurfaces that enclose regions where the A-block (minority component) local volume fraction is given by the colors on the colorbar with λW = 0.0 (top) and λW = 6.0 (bottom). The solid black surfaces enclose regions where the NR volume fraction is greater than 5ϕNR. (c) Density isosurfaces of NRs from the systems in (b) normalized by ϕNR are given by the colors on the colorbar with λW = 0.0 (top) and λW = 6.0 (bottom).

Numerical Methods and Parameters. We adopted the complex Langevin (CL) scheme2,28,30 to sample the fully fluctuating chemical fields using the equations ∂w±(r) ∂t

= −λ±

δ/ + η±(r, t ) δw±(r)

experimental systems), and unless otherwise noted, the volume fraction of NRs, ϕNR, is set to 0.5%. The length scales in all of the results presented below are scaled by Rg unless otherwise specified. The x and y dimensions of the simulation box were optimized under the mean-field condition to accommodate two unit cells of the bulk diblock chain in the cylindrical or lamellar phase. We do not expect this to be a significant approximation, as recent studies have shown that the fluctuation corrections to the domain spacing in neat block copolymers are expected to be small35,39 at our chosen value of C = 51; similar trends have also been found in polymer nanocomposites.40 The spatial resolution dx is set to 0.26, whereas the orientation of NRs is resolved in spherical coordinates with 24 intervals in the azimuthal angle ϕ and 12 intervals in the polar angle θ. Integrations of orientation u are performed using the Gaussian−Legendre integration scheme. Given the NR concentration is very low, we assume that there is no macrophase separation into particle-rich and particle-poor regions at these concentrations; explicit tests of this assumption would require either work in a grand-canonical ensemble or the Gibbs ensemble41−43 and is the subject of ongoing investigations. In other words, we assume that the mixture is at least metastable at all concentrations studied, and we do not observe evidence for phase separation.

(26)

where η(r, t) is a Gaussian white noise with the statistics

⟨η±(r, t )⟩ = 0

(27)

⟨η±(r, t )η±(r′, t ′)⟩ = 2λ± δ(r − r′) δ(t − t ′)

(28)

The CL results are averaged over 50 000 iterations with step size λ± δt = 10−5. We also note that when the Gaussian white noise is turned off, the above CL equations drive the system to the meanfield solution of the model. Our implementation of both the CL simulations and the SCFT calculations is achieved via a pseudospectral method, and all convolutions (e.g., those involved in iterating the Chapman−Kolmogorov equation) are evaluated in Fourier space. The field updates in eq 26 are performed using a first-order semi-implicit algorithm.38 The calculations in the study use the parameters χA,BN = 40, κN = 60, fA = 0.3, and ξ = 0.26Rg. The chain density C=

ρo R g3 N

= 51 (which we note is larger than common E


Article


Figure 2. (a) Fraction of vertically oriented nanorods f V as a function of the NR radius RP scaled by the cylinder domain diameter DA. (b) Isosurfaces enclose regions where the A-block (minority component) local volume fraction is given by the colors on the colorbar with RP = 0.61 (top) and RP = 0.39 (bottom) and fixed LP = 3.06, whereas the solid black surfaces enclose regions where the NR volume fraction is greater than 5ϕNR. (c) Density isosurfaces of NRs from the systems in (b) normalized by ϕNR are given by the colors on the colorbar with RP = 0.61 (top) and RP = 0.39 (bottom).

To map the above-mentioned dimensional quantities in our model to parameters in real diblock thin films, one can consider the PS-b-P2VP thin films in our recent work16 as an example. The PS-b-P2VP (Mn,PS = 180 kg/mol-b-Mn,P2VP = 77 kg/mol) there has a statistical segment length bexpt ≈ 0.67 nm44 and segments Nexpt ≈ 2460. Mapping this polymer to our diblock chain model ( N = 4 0 ) g i v es t h e s t a t is t ic a l s eg m e n t le n g t h , b≈

2 Nexptbexpt

N

from either surface was considered horizontal. All other NRs were considered vertical, and f V was then simply the ratio of the number of vertical NRs to the total number of NRs. As a result, f V has the following expression ⎡

fV =

= 5.25 nm and Rg ≈ 13.4 nm. In this scenario,

(

̃ (r, u) × ⎢⎣H uz − ∫ dr∫ du ρNR

2 2

) + H( − u

z

−

2 2

)⎤⎥⎦

̃ (r, u)⟩ ∫ dr∫ du⟨ρNR

(29)

the NRs considered in the article can be mapped to real NRs with lengths (LP) and radii (RP) in the ranges of 37−41 and 5−8 nm, respectively.

where H(x) is the Heaviside step function with unit value for x > 0 and zero otherwise and is used to select for NRs with a vertical orientation. Figure 1a shows how the vertical fraction changes as a function of the wetting parameter λW for three different NR lengths and for both rigid and soft confinements at the top surface. When the substrate is neutral to both the blocks (λW = 0), a finite value of f V is observed for NRs with a length comparable to a half of the film thickness, LP = 2.84. As λW increases and the substrate increasingly wets the cylinder domain, f V decreases and reaches zero at λW = 6.0, indicating all of the NRs are oriented parallel with the substrate. The distortion of the cylindrical domains near the substrate appears to attract the NRs, similar to the attraction of additives to grain boundaries in BCPs. The density maps of the A-block and the NRs at the two terminal λW values are commensurate with the change in the NR orientation (Figure

■

RESULTS Substrate Selectivity. We first study how substrate selectivity on the two blocks of the BCP affects the NR distribution and orientation in the thin films, and we focus on the NR orientation in the cylinder domain. To quantify the NR orientation, we calculated the fraction of NRs that are vertically oriented, f V, as functions of the wetting parameter, λW. The only locations in our films where we observe a finite population of horizontal NRs is near the substrates in our system. As such, we chose an angle of 45° as the cutoff between horizontal and vertical, and any NR whose center was within a distance of LP / 2 2

F


Article


Figure 3. (a) Local density of the NR centers in the thin film for the systems shown in Figure 2b. (b) ΔSvh as a function of the NR radius RP scaled by the cylinder domain diameter DA.

Figure 4. (a) Fraction of vertically oriented nanorods f V as a function of film thickness h. (b) Local density profiles of the NR centers in the thin film for selected cases. (c) Angle ⟨θ⟩ between the long axis of NRs that are not horizontally oriented and the norm of the substrate is calculated as a function of h. RP is fixed at 0.54 for all results presented in this figure.

The length of the NR also affects the NR orientation as the wetting conditions are changed, and longer nanorods are more likely to be vertically oriented. When the NR length LP reaches about 80% of the film thickness, the fraction f V becomes saturated and is invariant against λW for the range shown in Figure 1a. For stronger wetting with λW > 6, the vertical orientation of the cylinders becomes unstable for all NR lengths,

1b,c). The ordering of NRs near the substrate is quantified by calculating the largest eigenvalue of the nematic order tensor, Q(r), for the horizontally oriented NRs. For both the neutral (λW = 0) and A-block-selective (λW = 6.0) cases, the average of the largest eigenvalue is less than 0.1. This small eigenvalue indicates that the level of ordering is trivially small, as expected for the low NR concentration in our systems. G


Article


Figure 5. (a) Schematics show two sets of typical density isosurface figures for the thin films with vertically (left column) and horizontally (right column) oriented cylindrical phases. The figures in the top row show scaled A-block density isosurfaces, whereas the lower row includes the density isosurfaces from NRs scaled by ϕNR. RP = 0.61 and LP = 2.98 are fixed for NRs here. (b) Per-chain free-energy difference, ΔFvh, between vertically and horizontally oriented cylindrical phases as a function of the NR fraction ϕNR for varied film thickness h. All of the data are shifted by subtracting the value for pure polymer films, and the original values are shown in the inset. (c) Per-chain free-energy difference, ΔFvh from pure polymer cases (ϕNR = 0) as a function of h.

thin films. For two values of LP, Figure 2a−c shows how the vertical fraction f V changes with NR radius. We observe that although the nanorod lengths LP are quite different, both f V’s increase from approximately 15% to nearly 90% as RP increases from 0.34 to 0.61Rg. To further elucidate how the NRs are distributed in the cylinder domain, we plot the local density profiles of the NR centers across the film in Figure 3a, and two peaks are observed in the profiles. The peaks approximately RP away from the supporting substrates correspond to NRs oriented in the plane of the film, whereas the other peaks located LP/2 away from the substrates can be attributed to vertical NRs with one end adjacent to the substrates. The latter type shows a pronounced intensity for the larger RP, which indicates that the

and we observe a dense A layer near the bottom substrate with horizontal cylinders through the center of the film. Finally, we have tested whether confinement with a soft surface, which we model as an immiscible homopolymer, would lead to changes in f V as a function of the wetting conditions. Presumably, this soft confinement is more akin to a polymer−air interface than the rigid confinement typically used in field-theoretic simulations, and similar approaches have been used in other studies of block copolymer thin films.36 We find that a smaller fraction of the nanorods adopts a vertical orientation for the two shorter NR lengths (LP = 2.84 and 3.41), but there are no substantial qualitative changes in behavior. NR Geometry. We next explored how the NR radius RP and length LP would affect the NR distribution and orientation in the H


Article


Figure 6. (a) Fraction of vertically oriented NRs f V as a function of the nanorod radius RP. (b) Isosurfaces enclose regions, where the A-block (minority component) local volume fraction is given by the colors on the colorbar with RP = 0.58 (top) and RP = 0.39 (bottom) and fixed LP = 2.84, whereas the solid black surfaces enclose regions where the NR volume fraction is greater than 5ϕNR. (c) Density isosurfaces of NRs from the systems in (b) normalized by ϕNR are given by the colors on the colorbar with RP = 0.58 (top) and RP = 0.39 (bottom).

thickness, we first calculate f V as a function of h, as shown in Figure 4a for two NR lengths. The results show that overall f V does not have a strong dependence on the film thickness and remains in the range of 1.0 and 0.6, which indicates that the majority of NRs are vertically orientated. Because both the NR lengths are larger than DA (LP/DA = 1.9 and 2.1), the shorter NR has a smaller enthalpic penalty in the horizontal orientation than the longer NR. This might contribute to the finite decrease in f V for the shorter NR when the film thickness becomes large enough to hold two end-to-end-aligned NRs. Other than the finite change, the majority of vertical orientation can be confirmed by looking at the density profiles of NR center in Figure 4b, which feature a single peak in the center of the thickness when h = 2.5. When h increases to 5.6, two peaks associated with vertically orientated NRs with their ends placed near the boundaries of the film are observed in the profiles; in addition, there are less intense peaks corresponding to horizontal NRs laying near the substrates. It is surprising to see that f V remains near 1 even when the film thickness is slightly smaller than the NR length LP. To investigate this, the averaged angle ⟨θ⟩ between the long axis of NRs that are not horizontally oriented

NRs with larger radius are more likely to adopt a vertical orientation. To investigate the thermodynamics behind the effect of varying RP, we adopt the HPFT model described above and extract the entropy of the diblock copolymers. We calculate the per-chain entropy when we fix a single, explicit nanorod at the most probable location, where the NRs have a vertical orientation Sv, and from this, we subtract the per-chain entropy when a single NR is fixed at the most likely horizontal configuration Sh. We find that ΔSvh monotonically increases with RP for both choices of LP, indicating that the vertical orientation becomes entropically favored as for larger NRs (Figure 3b). We note here that our calculations assume at least metastable miscibility between the NRs and block copolymers, but according to a study on nanoparticles in bulk diblock copolymers,45 NRs might phase-separate from the block copolymer if they are sufficiently large in radius compared to the cylinder domain diameter, DA, and have a volume fraction above ∼5%. Film Thickness. In addition to the NR geometry, the thin film thickness h may also play a role in determining the NR distribution and orientation. To explore the role of the film I


Article


Figure 7. (a) Local density profiles of the NR centers in the thin film for the systems shown in Figure 6b. (b) Difference in the diblock chain entropy as a function of the NR radius RP.

and the norm of the substrate is calculated as a function of h using the following expression ⎡

⟨θ ⟩ =

(

) cos (u ) + H(−u − )(180° − cos ⎡ (r, u) × ⎢⎣H(u − ) + H(−u − )⎤⎥⎦

̃ (r, u) × ⎢⎣H uz − ∫ dr∫ du ρNR ̃ ∫ dr∫ du ρNR

2 2

−1

z

z

where H(x) is the Heaviside step function with unit value for x > 0 and zero otherwise. The results in Figure 4c show that as h decreases toward LP, the value of ⟨θ⟩ increases. When the value h is smaller than LP (e.g., h = 2.5), NRs with the two lengths exhibit θ with values larger than 20°. This result plus the density profiles indicate that the NRs still reside in the center of the film but becomes tilted toward the surfaces when the thickness begins to decrease below the NR length. NR Influence on Cylinder Orientation. In addition to understanding the parameters that control the NR distribution and orientation, it is also of interest to investigate how the NRs affect the thin-film morphology. Specifically, we studied how the inclusion of NRs influences the cylinder domain orientation in the thin films. We calculate the free-energy difference (ΔFvh) between two cylinder domain orientations, vertically and horizontally oriented relative to the substrate, to quantify NRs’ influence on the two morphologies (see Figure 5a). In Figure 5b, ΔFvh is plotted as a function of NR volume fraction, ϕNR, at different film thickness levels. When we shift ΔFvh by the values from pure polymer cases, we observe that overall ΔFvh decreases as ϕNR grows and the reduction in ΔFvh depends on the film thickness in a nonmonotonic manner. Our results indicate that there is an optimal film thickness that maximizes the reduction in ΔFvh due to the inclusion of NRs. A previous study on BCP thin films by Panagiotopoulos et al.46 proposed that one of the governing factors in the cylinder orientation is the interfacial energy between the unlike domains. When the film thickness exceeds the equilibrium spacing, Dhex, needed to accommodate a layer of cylinder domains, the vertical orientation becomes energetically preferred over the horizontal orientation. This view is compatible with the decreasing of unshifted ΔFvh in the NR free cases, as shown in Figure 5c, where the ratio h/Dhex increases from 1.41 to 1.61 as h increases from 4.65 to 6.04. The inset of Figure 5b shows that the vertical orientation is preferred for thicker films and higher NR concentrations, although the

z

2 2

2 2

z

2 2

⎤ (uz))⎥⎦

−1

(30)

preference for vertically oriented cylinders is nonmonotonic in the film thickness, similar to the results for neat diblocks.46 When the substrate becomes selective, equilibrium domain orientation in neat diblock thin films becomes even more complex and beyond the scope of the current study. We refer readers interested in this topic to previous studies.47,48 Lamellar Thin Films. Finally, we briefly show that trends similar to those found in the cylinder phase are observed when we vary the NR geometry in lamellar thin films. Figure 6a plots f V as a function of RP at three values of LP. The results show that f V increases sharply from approximately 0 to above 0.6 as RP increases from 0.34 to 0.6, and the effect is more pronounced for longer NRs. Two cases from Figure 6a with distinct RP values are chosen to highlight the difference in the NR orientation and are shown in Figure 6b,c. The NRs with larger RP are found primarily near the two substrates with a vertical orientation, whereas for smaller RP, the NRs are found horizontally aligned along the substrates. The change in the density profiles of the NR centers (Figure 7a) is commensurate with the change in the NR distribution, and the larger RP NRs have a pronounced peak corresponding to the vertically oriented NRs. As in the cylindrical phase, the thermodynamic reason behind the change can be traced back to the entropy of the diblock chains. Figure 7b plots ΔSvh as a function of RP and shows that ΔSvh increases with RP for both NR lengths considered here.

■

CONCLUSIONS In this article, we have used complex Langevin field-theoretic simulations and hybrid particle field theory calculations to examine the distribution of nanorods in block copolymer thin films as a function of NR dimensions, film thickness, substrate selectivity for the two blocks, and the phase of the diblock copolymer. We find that preferential wetting of the substrate for the cylinder-forming block can lead to segregation of the NRs to the substrate, and the NR dimensions play an important role in J


Article

The Journal of Physical Chemistry B determining the preferential location of the NRs in the film. In addition, we also performed thermodynamic integration to understand how the presence of NRs affected the relative stability of vertically oriented cylindrical domains compared to horizontally oriented domains. In general, we find that the entropy of the block copolymers plays a prominent role in determining the location of the NRs in thin films, and we expect that our results will be important for understanding and controlling the equilibrium distribution of anisotropic particles in block copolymer films.

■

(13) Thorkelsson, K.; Mastroianni, A. J.; Ercius, P.; Xu, T. Direct Nanorod Assembly using Block Copolymer-Based Supramolecules. Nano Lett. 2012, 12, 498−504. (14) Thorkelsson, K.; Nelson, J. H.; Alivisatos, A. P.; Xu, T. End-toEnd Alignment of Nanorods in Thin Films. Nano Lett. 2013, 13, 4908− 4913. (15) Deshmukh, R. D.; Liu, Y.; Composto, R. J. Two-Dimensional Confinement of Nanorods in Block Copolymer Domains. Nano Lett. 2007, 7, 3662−3668. (16) Rasin, B.; Chao, H.; Jiang, G.; Wang, D.; Riggleman, R. A.; Composto, R. J. Dispersion and Alignment of Nanorods in Cylindrical Block Copolymer Thin Films. Soft Matter 2016, 12, 2177−2185. (17) Kim, Y.; Chen, H.; Alexander-Katz, A. Free Energy Landscape and Localization of Nanoparticles at Block Copolymer Model Defects. Soft Matter 2014, 10, 3284−3291. (18) Listak, J.; Bockstaller, M. R. Stabilization of Grain Boundary Morphologies in Lamellar Block Copolymer/Nanoparticle Blends. Macromolecules 2006, 39, 5820−5825. (19) Ryu, H. J.; Sun, J.; Avgeropoulos, A.; Bockstaller, M. R. Retardation of Grain Growth and Grain Boundary Pinning in Athermal Block Copolymer Blend Systems. Macromolecules 2014, 47, 1419−1427. (20) Kumar, S. K.; Ganesan, V.; Riggleman, R. A. Perspective: Outstanding Theoretical Questions in Polymer-Nanoparticle Hybrids. J. Chem. Phys. 2017, 147, No. 020901. (21) Ganesan, V.; Jayaraman, A. Theory and Simulation Studies of Effective Interactions, Phase Behavior and Morphology in Polymer Nanocomposites. Soft Matter 2014, 10, 13−38. (22) Fredrickson, G. H. The Equilibrium Theory of Inhomogeneous Polymers; Oxford University Press: New York, 2009. (23) Matsen, M. W. The Standard Gaussian Model for Block Copolymer Melts. J. Phys.: Condens. Matter 2001, 14, R21−R47. (24) Leibler, L. Theory of Microphase Separation in Block Copolymers. Macromolecules 1980, 13, 1602−1617. (25) Shou, Z.; Buxton, G. A.; Balazs, A. C. Predicting the SelfAssembled Morphology and Mechanical Properties of Mixtures of Diblocks and Rod-Like Nanoparticles. Compos. Interfaces 2003, 10, 343−368. (26) Tang, Q.-Y.; Ma, Y.-Q. Self-Assembly of Rod-Shaped Particles in Diblock-Copolymer Templates. J. Phys. Chem. B 2009, 113, 10117− 10120. (27) Ting, C. L.; Composto, R. J.; Frischknecht, A. L. Orientational Control of Polymer Grafted Nanorods. Macromolecules 2016, 49, 1111− 1119. (28) Koski, J.; Chao, H.; Riggleman, R. A. Field Theoretic Simulations of Polymer Nanocomposites. J. Chem. Phys. 2013, 139, No. 244911. (29) Chao, H.; Hagberg, B. A.; Riggleman, R. A. The Distribution of Homogeneously Grafted Nanoparticles in Polymer Thin Films and Blends. Soft Matter 2014, 10, 8083−8094. (30) Koski, J.; Chao, H.; Riggleman, R. A. Predicting the Structure and Interfacial Activity of Diblock Brush, Mixed Brush, and Janus-Grafted Nanoparticles. Chem. Commun. 2015, 51, 5440−5443. (31) Wang, Z.-G. Fluctuation in Electrolyte Solutions: The Self Energy. Phys. Rev. E 2010, 81, No. 021501. (32) Villet, M. C.; Fredrickson, G. H. Efficient Field-Theoretic Simulation of Polymer Solutions. J. Chem. Phys. 2014, 141, No. 224115. (33) Helfand, E. Theory of Inhomogeneous Polymers: Fundamentals of the Gaussian Random-Walk Model. J. Chem. Phys. 1975, 62, 999. (34) Matsen, M. W. Thin Films of Block Copolymer. J. Chem. Phys. 1997, 106, 7781−7791. (35) Delaney, K. T.; Fredrickson, G. H. Recent Developments in Fully Fluctuating Field-Theoretic Simulations of Polymer Melts and Solutions. J. Phys. Chem. B 2016, 120, 7615−7634. (36) Stasiak, P.; McGraw, J. D.; Dalnoki-Veress, K.; Matsen, M. W. Step Edges in Thin Films of Lamellar-Forming Diblock Copolymer. Macromolecules 2012, 45, 9531−9538. (37) Lennon, E. M.; Katsov, K.; Fredrickson, G. Free Energy Evaluation in Field-Theoretic Polymer Simulations. Phys. Rev. Lett. 2008, 101, No. 138302.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Robert A. Riggleman: 0000-0002-5434-4787 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge support from the National Science Foundation through award DMR-1410246 (H.C. and R.A.R.) and the donors to the American Chemical Society Petroleum Research Fund for support through award 55662ND7 (B.J.L. and R.A.R.). Computational resources were also made available at the Texas Advanced Computing Center through XSEDE award TG-DMR150034.

■

REFERENCES

(1) Moll, J. F.; Akcora, P.; Rungta, A.; Gong, S.; Colby, R. H.; Benicewicz, B. C.; Kumar, S. K. Mechanical Reinforcement in Polymer Melts Filled with Polymer Grafted Nanoparticles. Macromolecules 2011, 44, 7473−7477. (2) Chao, H.; Riggleman, R. A. Effect of Particle Size and Grafting Density on the Mechanical Properties of Nanocomposites. Polymer 2013, 54, 5222−5229. (3) Cheng, Y.; Yu, A. S.; Li, X.; Oh, T.-S.; Vohs, J. M.; Gorte, R. J. Preparation of SOFC Cathodes by Infiltration into LSF-YSZ Composite Scaffolds. J. Electrochem. Soc. 2016, 163, F54−F58. (4) Yamamoto, R.; Yoshida, Y.; Yoshikawa, T.; Yajima, M.; Seki, T. Novel Thermally Conductive Sheet Applying Orientation Control of Graphite Particles. J. Jpn. Inst. Electron. Packag. 2010, 13, 462−468. (5) Hore, M. J. A.; Composto, R. J. Functional Polymer Nanocomposites Enhanced by Nanorods. Macromolecules 2014, 47, 875−887. (6) Liu, Y.; Mills, E. N.; Composto, R. J. Tuning Optical properties of Gold Nanorods in Polymer Films through Thermal Reshaping. J. Mater. Chem. 2009, 19, 2704. (7) Thompson, R. B.; Ginzburg, V. V.; Matsen, M. W.; Balazs, A. C. Predicting the Mesophases of Copolymer-Nanoparticle Composites. Science 2001, 292, 2469−72. (8) Balazs, A. C.; Emrick, T.; Russell, T. P. Nanoparticle Polymer Composites: Where Two Small Worlds Meet. Science 2006, 314, 1107− 10. (9) Bockstaller, M. R.; Mickiewicz, R. A.; Thomas, E. L. Block Copolymer Nanocomposites: Perspectives for Tailored Functional Materials. Adv. Mater. 2005, 17, 1331−1349. (10) Listak, J.; Bockstaller, M. R. Stabilization of Grain Boundary Morphologies in Lamellar Block Copolymer/Nanoparticle Blends. Macromolecules 2006, 39, 5820−5825. (11) Koski, J.; Hagberg, B.; Riggleman, R. A. Attraction of Nanoparticles to Tilt Grain Boundaries in Block Copolymers. Macromol. Chem. Phys. 2016, 217, 509−518. (12) Sides, S. W.; Kim, B.; Kramer, E.; Fredrickson, G. Hybrid ParticleField Simulations of Polymer Nanocomposites. Phys. Rev. Lett. 2006, 96, No. 250601. K


Article

The Journal of Physical Chemistry B (38) Lennon, E. M.; Mohler, G. O.; Ceniceros, H. D.; García-Cervera, C. J.; Fredrickson, G. H. Numerical Solutions of the Complex Langevin Equations in Polymer Field Theory. Multiscale Model. Simul. 2008, 6, 1347−1370. (39) Vorselaars, B.; Stasiak, P.; Matsen, M. W. Field-Theoretic Simulation of Block Copolymers at Experimentally Relevant Molecular Weights. Macromolecules 2015, 48, 9071−9080. (40) Koski, J. P.; Riggleman, R. A. Field-Theoretic Simulations of Block Copolymer Nanocomposites in a Constant Interfacial Tension Ensemble. J. Chem. Phys. 2017, 146, No. 164903. (41) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813−826. (42) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Phase Equilibria by Simulation in the Gibbs Ensemble. Mol. Phys. 1988, 63, 527−545. (43) Riggleman, R. A.; Fredrickson, G. H. Field-Theoretic Simulations in the Gibbs Ensemble. J. Chem. Phys. 2010, 132, No. 024104. (44) Eitouni, H. B.; Balsara, N. P. Thermodynamics of Polymer Blends. In Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Springer: New York, NY, 2006; Chapter 19, pp 339−356. (45) Huh, J.; Ginzburg, V. V.; Balazs, A. C. Thermodynamic Behavior of Particle/Diblock Copolymer Mixtures: Simulation and Theory. Macromolecules 2000, 33, 8085−8096. (46) Nikoubashman, A.; Register, R. A.; Panagiotopoulos, A. Z. SelfAssembly of Cylinder-Forming Diblock Copolymer Thin Films. Macromolecules 2013, 46, 6651−6658. (47) Wang, Q.; Nealey, P.; de Pablo, J. Monte Carlo Simulations of Asymmetric Diblock Copolymer Thin Films Confined Between Two Homogeneous Surfaces. Macromolecules 2001, 34, 3458−3470. (48) Yang, Y.; Qiu, F.; Zhang, H.; Yang, Y. Cylindrical Phase of Diblock Copolymers Confined in Thin Films. A Real-Space Self-Consistent Field Theory Study. Polymer 2006, 47, 2205−2216.

L


Field-Theoretic Simulations of the Distribution of Nanorods in Diblock

Recommend Documents