Non-spherical nanoparticles in block copolymer composites: nanosquares, nanorods and diamonds

A hybrid block copolymer(BCP) nanocomposite computational model is proposed to study nanoparticles(NPs) with a generalised shape including squares, rectangles and rhombus. Simulations are used to study the role of anisotropy in the assembly of colloids within BCPs, ranging from NPs that are compatible with one phase, to neutral NPs. The ordering of square-like NPs into grid configurations within a minority BCP domain was investigated, as well as the alignment of nanorods in a lamellar-forming BCP, comparing the simulation results with experiments of mixtures of nanoplates and PS-\textit{b}-PMMA BCP. The assembly of rectangular NPs at the interface between domains resulted in alignment along the interface. The aspect ratio is found to play a key role on the aggregation of colloids at the interface, which leads to a distinct co-assembly behaviour for low and high aspect ratio NPs.


Introduction
Polymer nanocomposite materials containing anisotropic nanoparticles have attracted a lot of attention due to their ability to create functional materials with enhanced properties 1 . Block copolymers can serve as perfect matrices to produce highly ordered polymer nanocomposite materials due to the well-defined, periodic structures that block copolymer melts have. This periodicity makes block copolymers excellent matrices to host nanoparticles (NPs), which can be localised in specific regions of the phase-separated block copolymer 2, 3 .
The co-assembly of block copolymer melts and colloidal nanoparticles has long been studied with theoretical 4,5 and computational [6][7][8][9][10][11] methods, as well as experiments [12][13][14][15] . The resulting nanocomposite material is not merely the addition of a passive filler. Instead, the properties of the polymeric matrix can be modified by the presence of colloids and the assembly of colloids is deeply intertwine with the diblock copolymer matrix 16 . For instance, colloids which are wetted to be compatible with one of the phases will drive the BCP into a morphological transition 6,17,18 .
While spherically shaped colloids have been widely studied -both experimentally and theoretically 10the additional difficulty of modelling anisotropy have resulted in comparably less work devoted to it. Nevertheless, the orientational degree of freedom of non-spherical colloids introduces new possibilities of BCP/NP co-assembly, given the intrinsic ordered structures of the neat block copolymer (lamellar, cylindrical, etc). Complex non-spherical nanoparticles such as nanorods(NR) have attracted a lot of attention as constituents of functional polymer nanocomposite materials 1 . For instance, gold nanorods orient along the lamellar domain when confined in one of the symmetrical phases 19,20 . Similarly, gold nanorods template the direction of the cylindrical domains in an asymmetrical diblock copolymer mixture 21 . Ordered arrays of aligned nanorods were achieved experimentally 22,23 in the co-assembly of BCP and highly anisotropic NPs, where NRs were organised in an end-to-end configuration.
Theoretical and computational works have studied the self-assembly of BCP and anisotropic NP. Dissipative Particle Dynamics (DPD) has been largely used, thanks to the ability to combine several beads into rod-like sequences. The phase behavior of such systems has been studied along with the orientation of nanoparticles [24][25][26] and the effect of shear in the global orientation 27 . Osipov et al [28][29][30] used Strong and Weak Segregation Limit Theory to determine the distribution of anisotropic particles in a diblock copolymer, with a low fraction of nanoparticles present in the system. A mesoscopic Ginzburg-Landau approach has been used to study Janus nanorods mixed with BCP 31 .
In this work our previous hybrid BCP/NP model 9 originally developed for spherical nanoparticles is extended for a vast family of NP shapes. The shape of the NPs is motivated by experiments. For example, considerable attention has been given to nanocubes 32,33 as an example of faceted particles that can lead to the formation of crystals 34 . In this work we make use a considerably fast, coarse grained Cell Dynamic Simulation (CDS) method to describe the block copolymer dynamics, while colloidal NPs are described using Brownian Dynamics.
The Cell Dynamic Simulation method has been used extensively both in pure block copolymer systems [35][36][37][38] and nanocomposite systems 9 , reproducing experiments such as aggregation of incompatible colloids 16,39 and NP-induced phase transitions 40 . Its relative computational speed makes it suitable to study properties which involve large systems over extended times, while the phenomenological approach in its model limits its validity in the microscopic realm. This hybrid method permits to explore the high NP concentration regime, in which the presence of NPs introduces huge differences over the neat BCP matrix, such as, morphological phase transitions.
In the following section a generalised model for anisotropic NPs will be presented. To our knowledge, this work represents the first mesoscopic simulations of BCP/NP assembly for cuboid, rectangular or diamond shaped nanoparticles. Contrary to previous works 31 , we aim to capture the finite width of rectangular NPs, thus allowing to study the role of the aspect ratio in the high NP concentration regime.

Model
Here we present a coarse grained model of the dynamics of BCP and NPs with a range of shapes and chemical properties. This model is intrinsically mesoscopic, disregarding microscopic details, i.e. the individual monomers of the polymer chains are not resolved. This coarse grained approach permits to achieve relatively large box sizes compared to previous works.
The evolution of the BCP/colloids system is determined by the excess free energy which can be separated as with F pol being the free energy functional of the BCP melt, F cc the colloid-colloid interaction and the last contribution being the coupling term between the block copolymer and the colloids.
The diblock copolymer is characterized by the order parameter ψ(r, t) which represents the differences in the local volume fraction for the copolymer A and B with respect to the relative volume fraction of A monomers in the diblock, The order parameter must follow the continuity equation in order to satisfy the mass conservation of the polymer: If the polymer relaxes diffusely towards equilibrium, the order parameter flux can be expressed in the form j(r, t) = −M ∇µ(r, t) as a linear function of the order parameter chemical potential Introducing these equations into the continuity equation and taking into account the thermal fluctuations we obtain the Cahn-Hilliard-Cook equation (CHC) where M is a phenomenological mobility constant and ξ is a white Gaussian random noise which satisfies the fluctuation-dissipation theorem 41 . The copolymer free energy is a functional of the local order parameter which can be expressed in terms of the thermal energy k B T as where the first and second terms are the short and the long-range interaction terms respectively. The coefficient D is a positive constant that accounts for the cost of local polymer concentration inhomogeneities, the Green function G(r − r ) for the laplace Equation satisfies ∇ 2 G(r − r ) = −δ(r − r ), B is a parameter that introduces a chain-length dependence to the free energy 42 and H(ψ) is the local free energy 35,42 , where τ 0 , A, v, u are phenomenological parameters 35 which can be related to the block-copolymer molecular specificity. Previous works 9,35,43 describe the connection of these effective parameters to the BCP molecular composition. τ = −τ 0 + A(1 − 2f 0 ) 2 , D and B can be expressed 43 in terms of degree of polymerization N , the segment length b and the Flory-Huggins parameter χ(inversely proportional to temperature) . Subsequently, we will consider u and v constants 44 , which define all the parameters identifying the BCP local free energy H(ψ) . As previously shown 45,46 , CDS can be used along with more detailed approaches like dynamics self-consistent field theory (DSCFT), using CDS as a precursor in exploring parameter space due to the computationally inexpensiveness nature of CDS. We can express the time evolution of ψ , Equation 6, using CDS as r i being the position of the node i at a time tδt, and the isotropic discrete laplacian for a quantity X is given by 47 where δx is the lattice spacing. Specifically, we will use NN, NNN meaning nearest neighborsand next-nearest neighbors, respectively, for the two dimensional case.
In Equation 9 we have introduced the auxiliary function and also, the map function 35,48

Coupling between the block copolymer and nanoparticles
In previous works 9 , the presence of N p nanoparticles has been introduced via an extra term in the free energy, given by In the case of a simple circular particle, s = |s|/R = (x/R) 2 + (y/R) 2 and the tagged function is related to the core/shell nature of the particle. It is a smoothly monotonically decreasing function that takes values ψ c (0) = 1 and ψ c (1) = 0. In our choice, for s < 1 and ψ c (s > 1) = 0. This expression allows us to design a nanoparticle that interacts softly with the surrounding block copolymer, and it avoids the need to explicitly consider the NP-BCP boundary.
To extend the colloidal shape from a circular particle to an elongated ellipse, one can trivially rescale the axis, such that s = (x/a) 2 + (y/b) 2 for an ellipse that is resting horizontally (vertically) if a > b (a < b). Furthermore, we can consider a generalisation of this procedure, to represent a larger family of curves, namely, superellipses. Such closed curves are given by where the exponent 1/n rescales the decay 49 of ψ c . More generally, each rotated anisotropic colloid is characterised by an unit vectorn i = (cos φ i , sin φ i ) depending on an angle φ i . Again, x, y should be related to the rotated x (φ), y (φ) versions of the variables. In Fig. 1 representative examples of super-ellipses are given for several values of n, for two values of the aspect ratio. The simplicity of the coupling free energy allows us to simulate a vast range of NP shapes. The family of superellipses in two dimensions include ellipses, rectangles, rhombus and star-like particles. Many of these shapes can serve as a faithful representation of experimentally-motivated non-spherical colloids such as nanocubes, nanospheres or nanorods.

Interparticle potential
The calculation of forces and torques requires a choice of a suitable intercolloidal pairwise additive potential. In the particular case of ellipsoids (n = 1), a modified Lennard-Jones potential has been widely used 50 . Nonetheless, in this work we will study the particular case of rhomboid and rectangular shaped particles. To this end, we design a simplified repulsive potential that will assure non-overlapping between particles.
We introduce a colloid-colloid contribution to the free energy as with U (r i , r j , φ i , φ j ) a potential that depends on the distance between two colloidal particles, as well as their orientations φ i , φ j . Specifically, we choose a potential that is proportional to the overlapping area (volume in three dimensions) between two colloidal particles of arbitrary shape which prevents overlapping between colloids. To our knowledge, there is no general analytic expression of the overlapping area between two parallelogram of arbitrary orientation and distance. Nonetheless, the Separating Axis Theorem (SAT) is widely used in computer science to test collisions between convex polygons 51 .
In order to derive forces and torques between colloids, we require an expression of the overlapping area. To this end we numerically calculate the overlapping area between colloids and perform a fit of the results into simple analytical expressions. Sines and cosines are used to capture the appropriate symmetry of each object. Additionally, the overlapping area is approximated to be proportional to the center-to-center distance. The occurrence of overlapping can be determined exactly and -since forces and torques are only calculated when overlapping is true-making this method computationally fast. The details for different shapes, orientations and distances can be found in the Supplementary Information.

Brownian Dynamics
Colloids undergo diffusive dynamics described by the Langevin equation in the overdamped regime. The centre of mass of each colloid follows Brownian Dynamics, that is, where f cpl i = −∂F cpl /∂R i is the coupling force derived from Eq. 13 and f cc i = −∂F cc /∂R i is the colloid-colloid force derived from Eq. 16 through the fitted expression of the overlapping area in Eq. 17. Similarly, the orientation of particle i follows where γ t and γ R are the translational and rotational friction coefficients and torques can be calculated as M cpl i = −∂F cpl /∂φ i and M cc i = −∂F cc /∂φ i , respectively for the coupling and colloid-colloid torques. Generally speaking, the translational diffusivity is a tensor that accounts for the anisotropic diffusivity of arbitrary shaped particles 52 along the parallel and perpendicular main axis of the colloid. Nonetheless, we are interested on the equilibrium properties of the assembly of polyhedra shaped colloids. Thus, we can assume these γ t to be a scalar 20 .

Results
The standard values of CDS 9,35,36 will be used τ 0 = 0.35, u = 0.5, v = 1.5, A = 1.5, D = 2.0 while the BCP/NP interaction is set to σ = 1.0. The BCP time scale is set by the mobility M = 1. A moderate BCP noise strength 0.1 is chosen while the NP temperature scale is set to k B T = 0.1. A cell spacing δx = 1.0 and time discretisation δt = 0.1 are chosen. Lengths will be presented in units of grid points. Initially, both the position and the orientation of colloids is random. Furthermore, the block copolymer is initially disordered, corresponding to a quick quench at t = 0.

Selective NPs
In order to study a minimal example of anisotropic colloids in BCP, we can firstly consider squarelike NPs which are compatible with the minority phase of the BCP (ψ 0 = −1). By doing so, we are introducing a high effective volume fraction of colloids. We select a cylinder-forming BCP with f 0 = 0.35 and explore a concentration of colloids φ p with a side length 2a 0 . At low concentrations, the NPs are confined within circular BCP domains. At higher concentration, the NPs are effectively increasing the fraction minority monomers, enough to induce a cylinder-to-lamellae transition in the BCP. 6,15,17,18,39,40,53 A combination of the entropic NP-NP interaction and the enthalpic effect of the surrounding block copolymer can lead to organisation of the square-shaped colloids into a defined sideto-side assembly. In order to quantify an eventual grid-like arrangement, we can define the interparticle orientational parameter where the brackets indicate average over all particles' neighbours. This parameter is 1 for a grid-like arrangement of particles and it is clearly defined positive. It provides insight both on the orientational and translational order of the colloidal set of particles, as it captures the π/2 rotation symmetry of squares. We can explore the parameter space given by the fraction of particles in the system φ p and the size of the particle a 0 . By doing so we can explore the arrangement of particles for a given level of constrain for different sizes. In Figure 2 we plot such curves. We can observe that small-sized particles are not well oriented within the block copolymer, with T < 0.5 for any fraction of particles when a 0 = 1.5. We can conclude that small particles tend to be dominated by thermal motion and neither the block copolymer, nor the entropy associate with the anisotropy of particles can induce an ordered phase. Larger particles, instead, can achieve ordered configurations with array-like organisation of squares. For a three-dimensional system, we can expect a monolayer of well-ordered cubes in a lamellar phase.
In Fig. 3 we study the assembly of rhombus-shaped NPs by selecting an exponent n = 0.6 from Eq. 15. A value f 0 = 0.35 and φ p = 0.2 can be selected, such that the BCP morphology is lamellar in the presence of a concentration of NPs. By doing so, rhomboid-shaped NPs segregate within the minority domains and experience a considerable local effective concentration, which leads to close-packed configuration for three different rhomboid sizes: a 0 = √ 2, 2 and 2 √ 2 for (a), (b) and (c), respectively. The NP aspect ratio is kept constant b/a = 1/2. Regardless of the size of the particles, the BCP lamellar morphology and its confinement effect leads to a highly ordered NP configuration.

Neutral NPs
Homopolymer chains can be grafted to the surface of NPs to create surfactant NPs that tend to segregate at the interface between BCP domains 54 . Setting a neutral affinity ψ 0 = 0, we can study the role of anisotropy in the assembly of rectangular NPs (n = 2). Because of the colloidal shape, rectangular nanoparticles tend to orient along the interface, as can be seen in Figures 4 (a) for an asymmetric, cylinder-forming BCP. In Figure 5 (a) we can see the orientation of colloids along the lamellar domains. In both cases, a large number of particles oriented in the interface leads to a transition induced by the saturation of the interface and by the formation of bridges along domains. Figure 4 shows the transition from a cylinder-forming BCP with circular domains organised hexagonally (a) into a blue-domain lamellar with a combination of minority (red) phase with colloidal-rich region (b). This transition was already reported in simulations of circular NPs and BCP 40 , but it can be noticed that many colloids still prefer to assemble at the interface and oriented along it.
A lamellar, symmetric BCP mixed with rectangular particles can be seen in Figure 5. In (a), a small concentration of particles is oriented along the interface. At larger concentrations the colloids form bridges along domains which break them into smaller, irregular domains. The chessboard-like arrangement of BCP domains can be noticed, which was not as clear in mixtures of circles/BCP 40 . This is clearly a consequence of the anisotropy of colloids, as can be noticed in the detail image in Figure 5 (b), where the rectangular particles tend to orient along interfaces which can enhance the formation of grid-like configuration of BCP.
The anisotropy of neutral nanoparticles has been shown to result in alignment along the interface. The effect that the aspect ratio of rectangles have on the assembly of colloids and its effect on the overall co-assembly of the BCP composite system can be studied. In particular, we can consider a lamellar-forming BCP (f 0 = 0.5) mixed with a concentration φ p of neutral nanoparticles with a larger  semiaxis sized a = 3 and an aspect ratio e = b/a. At low concentration, we can expect NPs to be preferentially segregated towards the interface. Figure 6 shows that even at low concentration weakly anisotropic NPs (e = 0.67 and e = 1.0) display a radically distinct behaviour from strongly anisotropic NPs (e = 0.33). The latter are mostly attached to the interface, leaving the BCP largely unmodified. This can be tracked by calculating the number of BCP domains in the system, which remains constant up to a high concentration value in which the interface is saturated and the lamellar domains are divided into shorter domains. Contrary to that, weakly anisotropic NPs are prone to form aggregates due to their comparatively lower trapping energy at the interface. For this reason we BCP morphology suffers a steady increase in the number of domains due to the break-up of domains following the formation of NP clusters. Figure 6: Number of BCP domains in terms of the fraction of particles in the system for several particle types. b = 1, 2, 3 represent 3 types of rectangular particles with a = 3 larger semiside and b being the shorter semiside.

Comparison with experiments
Anisotropic nanoparticles possess an additional orientational degree of freedom which-when mixed with phase-separating block copolymers-can induce an orientational order. Composto et al found nanorods oriented along the lamellar domains 19 while end-to-end, in-lamella assembly has been found in thin films 23 . More recently, nanoplates have been found to align within a lamellar-forming PS-b-PMMA block copolymer 55 . While this type of assembly is intrinsically three dimensional, we can aim to study the two-dimensional slice of a nanoplate into a rectangle with sizes d 1 and h, related to the larger diagonal and the thickness of a nanoplate. Along with the lamellar periodicity, we can simulate experimentally relevant parameters of such systems, as in Figure 7 where we observe a clear alignment of rectangular nanoparticles along the lamellar domains.
The horizontal cross-section of a nanoplate can be simulated as a rhomboid particle using experimental relative sizes 55 to simulate realistic rhomboids (nanoplates in 3D) with sizes b/a = 0.629 (ratio of diagonals) and the relative size between the lamella thickness and larger diagonal tuned accordingly. In Figure 8 top-left we can see how the rhomboids are placed within the circular domains, enlarging those in which more than one nanoparticle is placed. Furthermore, the rhomboids tend to assemble in a close-packed manner. If we increase the number of particles, this behaviour is even more clear, as in top-right. For a lamellar-forming BCP (bottom row), the nanoparticles are not totally confined , only within the stripe domains. Because the lamellar thickness cannot accommodate a single NP, they tend to aggregate in order to minimise the distortion induced by the presence of NPs. Figure 7: Alignment of rectangular nanoparticles in a symmetric block copolymer. Nanoparticle size is a 0 = 12.11 and b 0 = 1.21 for the major and minor semiaxis while the BCP lamellar spacing is L ≈ 8.

Conclusions
A general model to simulate anisotropic nanoparticles immersed in block copolymer has been introduced. This flexible scheme can cover the family of curves known as super ellipses which include ellipsoids, rectangles and rhomboids. A colloid-colloid potential that is proportional to the overlapping area between two objects has been used. This analytical potential allows to calculate forces and torques and thus study the time evolution of anisotropic NPs. Due to the computational efficiency of the CDS/BD model, this hybrid method can be used to study the co-assembly of BCP/anisotropic nanoparticles, exploring a large parameter space.
The assembly of square-shaped in a minority phase of a BCP has been studied for different values of concentration and NP size, finding that the enthalpic interaction of the BCP coupling dominates over the assembly of colloids at large particle sizes, compared with the BCP period. At high concentrations colloids were able to organise into a side-to-side configuration thanks to the orientational-dependent NP-NP interaction. Similarly, minority-compatible rhomboid-shaped NPs where found to align along the domain axis when mixed with a lamellar-forming BCP.
We have been able to compare experimental conditions by considering both rectangular NPs and rhomboid-shaped nanoparticles in order to mimic three-dimensional nanoplates. Using relative sizes from recent experiments 55 , we have found alingment of anisotropic colloids within the lamella domains, reproducing experiments involving nanorods 19,23 .
Contrary to A-selective NPs in a A/B BCP, we have studied the co-assembly of neutral, interfacecompatible nanoparticles. We have found that neutral anisotropic NPs are preferentially aligned along the interface, which results in a saturation of the interface at high concentration. This behaviour is nonetheless dependent on the aspect ratio, as a consequence of the strong trapping that anisotropic particles undergo at the interface. Weakly anisotropic NPs, on the contrary, have shown a tendency to aggregate which results in a distortion of the lamellar morphology.
In conclusion, we have presented a general approach to simulate multiple nanoparticle shapes mixed with a BCP matrix. Contrary to previous models, non-spherical nanoparticles interact anisotropically both with the surrounding medium and with themselves. In combination with an efficient computational method, this allows to simulate high concentration regimes in which the NP-NP interaction is crucial at determining the overall NP configuration.