Nematic Ordering of Anisotropic Nanoparticles in Block Copolymers

Block copolymer melts have been previously used to control the position and alignment of anisotropic nanoparticles. In this work, 2D and 3D mesoscopic simulations are used to explore the phase behavior of block copolymer/nanoparticle systems. The method combines a time‐dependent Ginzburg‐Landau for the polymer and Brownian dynamics for the anisotropic nanoparticles. Rhomboidal and spheroidal shaped particles are simulated in two and three dimensions, respectively. It is found that the nanoparticle nematic order aligned by the block copolymer domains enhances the lamellar phase of the block copolymer, due to an anisotropy‐driven phase transition. Additionally, anisotropic nanoparticles within circular‐forming block copolymer leads to a competition between the nematic colloidal ordering and the hexagonally ordered mesophase. At large concentrations, the nematic order dominates, deforming the block copolymer mesophase.


Introduction
Block copolymer(BCP) melts have been largely used as templates to control the position of colloidal nanoparticles (NPs) due to the ability of BCP to self-assemble into periodic ordered structures such as lamellae, hexagonally ordered cylinders or www.advancedsciencenews.com www.advtheorysimul.com for both isotropic nanospheres [35,36] and anisotropic NPs. [25,37] This implies that the BCP does not act as a mere passive template and the anisotropic NPs co-assembles with the BCP template. In order to achieve optimal NP alignment, we will explore the interplay between the BCP morphology and the NP assembly. Furthermore, experiments have shown additional changes in the BCP morphology due to the NP shape [25] which motivates a detailed study of the NP anisotropy and its relationship with the BCP mesophase. Moreover, previous experiments and simulations have focused on dilute or moderate concentrations of NPs. In this work, we use computer simulations to study the alignment of anisotropic NPs within BCP melts. We will explore considerably large NP concentrations, in which competition between the BCP morphology and NP nematic order emerges. For these objectives, it is suitable to make use of a relatively fast mesoscopic computational method that allows to reach large system sizes.

Model
In this section we present the model used to study block copolymer melts with anisotropic NPs. This is a mesoscopic hybrid model which permits to achieve relatively large box sizes compared to previous works. In a previous work, [34] we have presented the full 2D model, which is here summarized and briefly extended to 3D. Furthermore, this model has been exploited and compared with experiments [37] obtaining a considerable match.
The evolution of the BCP/colloids system is determined by the excess free energy which can be separated as  tot =  pol +  cc +  cpl (1) with  pol being the free energy functional of the BCP melt,  cc the colloid-colloid interaction and the last contribution F cpl 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 (r, t) = A (r, t) − B (r, t) + (1 − 2f 0 ) ( 2 ) with respect to the relative volume fraction of A monomers in the diblock, f 0 = N A ∕(N A + N B ). The concentration fluctuations follow diffusive dynamics, described by the Cahn-Hilliard-Cook equation (CHC), [38][39][40] (r, t) where M is a phenomenological mobility constant and is zero mean and unit variance white Gaussian random noise which satisfies the fluctuation-dissipation theorem. [41] The amplitude of the random noise is given by noise . The copolymer free energy is a functional of the local order parameter which can be expressed 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, [42,43] H( ) = 1 2 where 0 , A, v, u are phenomenological parameters [43] which can be related to the block-copolymer molecular specificity. Previous works [43][44][45] describe the connection of these effective parameters to the BCP molecular composition. ′ = − 0 + A(1 − 2f 0 ) 2 , D and B can be expressed [45] in terms of degree of polymerization N, the segment length b and the Flory-Huggins parameter.

Coupling between the Block Copolymer and Nanoparticles
The presence of N p particles is introduced by the coupling free energy term, [44] where 0 is the affinity controlling the preference of the NP with the medium which is experimentally related to the coating of the NP surface. [1][2][3] Meanwhile, the strength of the interaction is controlled by . The shape of the particle is specified by the tagged function c which describes the range of the NP-polymer interaction, the interior/exterior [46] and the orientation of the particle. [34] It is a monotonically decreasing function in s, with a cut-off in the surface of the NP, www.advancedsciencenews.com www.advtheorysimul.com to model particles with three main semi-axis a, b and c with b = c < a. For the sake of clarity, the x, y, z cartesian axis are chosen to align with the body-fixed coordinate system.

Interparticle Potential
The calculation of forces and torques requires a choice of a suitable intercolloidal pairwise additive potential. In order to assure non-overlapping of 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 . The specific form of the colloid-colloid potential depends on the NP shape, in particular, the subfamily of super-ellipses determined by the value of n. In 2D, in this work we are restricted to rhomboidal particles with n < 1. A completely repulsive potential which is proportional to the overlapping area between two arbitrary-shaped rhomboids has been presented in reference [34] and is used in this work. This potential captures the NP shape and is limited to instances of overlapping, which are determined exactly via the separation axis theorem, [48] which states that the overlapping of two convex objects in 2D is prevented if a line can be drawn separating both objects. The overall scale of the potential is tuned to prevent overlapping.
In the particular case of 3D ellipsoids (n = 1) the Gay-Berne(GB) potential, a Lennard-Jones potential adapted to ellipsoidal geometry, is widely used. [49] While this potential possess appropriate features-capturing the spheroid shape and possessing attractive and repulsive terms-we select a totally repulsive, Yukawa-like potential, as we are mainly interested in the NPpolymer interaction. In order to capture the anisotropy of the particles, we substitute the Yukawa diameter for the orientational dependent GB (r ij ,û i ,û j ) parameter present in GB potentials, leading to a potential of the form where r ij is distance between two particles' center of mass. This potential only acts upon particles when two colloids overlap. The modified diameter GB captures the NP shape [49] and, in the case of an ellipsoidal NP with n = 1, it describes the same shape as the polymer-colloid coupling shape function, c (s), as expressed in Equations (7) and (8).

Brownian Dynamics
Colloids undergo diffusive dynamics described by the Langevin equation in the overdamped regime. The center of mass of each colloid follows Brownian dynamics, that is, where f cpl i = −  cpl ∕ r i is the coupling force derived from Equation (6) and f cc i = −  cc ∕ r i is the colloid-colloid force derived from Equation (10). 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 = −  cpl ∕ i and M cc i = −  cc ∕ i , respectively for the coupling and colloidcolloid torques. In 3D the orientation of particle i can be generalized to the Euler angles i , i , i . Generally speaking, the translational diffusivity is a tensor that accounts for the anisotropic diffusivity of arbitrary shaped particles [50] along the parallel and perpendicular main axis of the colloid. Nonetheless, we are interested on the equilibrium properties of the assembly of complex shaped colloids. Thus, we can assume t and R to be a scalar. [51]

Computational Method
The dynamic CHC Equation (3) is solved with the Cell Dynamic Simulation(CDS) method which has been widely used to model the time evolution of BCP melts. [43,[52][53][54] While the detailed features of this approach have been widely shown before, [44] it is important to notice that the high isotropy of the Laplacian in the CDS scheme is due to the numerical expression  (14) in 3D, with NN, NNN, NNNN meaning nearest neighbors, nextnearest neighbors, and next-next-nearest neighbors, respectively. In 2D instead, In this work, we use standard parameters used in hybrid CDS/Brownian dynamics methods: the BCP parameters are 0 = 0.35, B = 0.002, D = 1.0, u = 0.5, v = 1.5, A = 1.5 and M = 1.0. A moderate noise term noise = 0.05 is introduced following the algorithm by Ball et al. [41] The CDS time and space discretization are set to t = 0.1 and x = 1.0.
The NP parameters are k B T = 0.1 and = 1.0. The affinity of the NP is set to 0 = −1, indicating that the NPs are totally miscible within one of the BCP phases.
The main length scales in this work are BCP interface thickness, L interface , the BCP periodicity H 0 and the various NP semiaxis. Lengths are displayed in units of the BCP periodicity H 0 . To facilitate reproducibility, the interface thickness throughout this work is L interface ∕H 0 = 0.09 which indicates a considerable degree of segregation. A particle-polymer interface width can be approximately defined based on the decay of the shape function c (s) is the distance such that c (a 0 ) = 1∕2 in along main axis. The length a 0 (b 0 ) approximately indicates the hard-core interior of the NP while the interface width a w indicates the width of the polymergrafted corona of the NP. Throughout this work, this length scale www.advancedsciencenews.com www.advtheorysimul.com is a w ∕L int ≈ 2.0(2a∕H 0 ) (b w ∕L int ≈ 2.0(2b∕H 0 )). Additionally, time scales are shown in units of the BCP microphase separation time scale t BCP = L 2 interface ∕(M ).

Observables
In this work, we make use of several quantities to quantify the coassembly of the hybrid system. For the sake of clarity, we provide a categorized list of the observables.

Block Copolymer Quantities
A domain analysis method is used to quantify the BCP morphology in 2D. The BCP interfaces can be numerically identified as the grid points such that the absolute value of the BCP field is smaller than a threshold value | (r)| < → r interface . Then, cluster analysis can be used to identify the different clusters, followed by a determination of the centroids of cluster (domain) and the shape of the domains. The number N domains and shape of the BCP domains can be calculated for each simulations to determine the fraction of total domains that are circular below a certain threshold, Φ circle . Thanks to the domain analysis method, we can obtain the center of mass of each BCP domain.
The hexatic bond order parameter is usually employed in colloidal physics to characterize the hexatic-that is, sixfold-ordering. It can be expressed as where ⟨⋅⟩ is an average over all particles i = 1…N. For each particle i, we determine the number of first neighbors N i through a cutoff, and calculate the bond angle ij . It takes values Ψ 6 = 1 and Ψ 6 = 0 for perfectly hexatic and random configurations, respectively. This parameter can be used to characterize the BCP hexagonal mesophase by applying Ψ 6 to the coordinates of the center of mass of each BCP domain, as determined by the domain analysis method. Furthermore, the bond angle parameter can be generalized to n-fold symmetry Ψ n = ⟨1∕N i ∑ j exp(ni ij )⟩. In 3D, domain analysis is computationally expensive and the phase behavior is richer than in 2D (non-planar morphologies include hexagonally ordered cylinders, gyroid and body-centered cubic spheres). We make use of the standard Minkowski functionals approach to study the topology of the BCP morphology. The Minkowski functionals are the volume, surface, curvature and Euler characteristic, and are calculated by converting the 3D field (r) into open and closed voxels. [55] The Euler characteristic provides topological information on the connectivity of BCP domains: It takes values = 1 for an isolated sphere, = 0 for a torus, and negative values for highly connected structures. [56] The 2D mean curvature of the BCP interfaces can be calculated by defining a local normal unit vector to the interfacen = ∇

||∇ ||
where || ⋅ || is the modulus of a vector. The tangential vector iŝ t and its change can be calculated with the tensorC = ∇t with components C ij = t i x j . Its projection into the tangential direction is related to the local curvature k(r) =C ⋅t. Finally, the total cur-vature is averaged over space with a weight ||∇ || to minimize the bulk contributions.

NP Quantities
The NP properties can be separated in two groups: particleparticle and particle-polymer order.
The colloid-colloid nematic order is extensively used throughout this work. The presence of the BCP field limits the emergence of global nematic order. Instead, we calculate the nematic order of a system as the average over all particles S cc = ⟨S i cc ⟩ where where N i is the number of first neighbors to particle i, which is specified by the particle-particle distance r ij < R * ij . The cut-off is chosen as R * ij = 3a. The nematic order indicates the orientation of particles relative to each other: a value S cc = 1 indicates parallel alignment, S cc = 0 random alignment and S cc = −1 indicates perpendicular alignment.
Similarly, we use a fourfold tetratic order parameter that introduces additional symmetry, S 4,i cc = (S i cc ) 2 . This parameter takes value 1 for both perpendicular and parallel alignment-that is, tetratic ordering-which is useful for particles with square shape which have rotational ∕4 symmetry.
Finally, S cpl (r) the NP orientation with respect to the center of the BCP domains is calculated as follows: first, the domain analysis method is used to determine the center of mass of each BCP domain, R . For each NP, we calculate the nematic contribution of particle i with respect to the BCP domain , and average over all particles and BCP domains where the unit vectorr i = (r i − R )∕||r i − R || is the relative vector between the NP the BCP domain center of mass. A radial histogram is produced for the distance r = ||r i || of the NP to the BCP domain center of mass.

Results
In this work we study the interplay between the BCP morphology and NP anisotropy. First, we will explore the alignment of rhomboidal NPs within BCP melts, focusing on the lamellar morphology and explore the emergence of nematic order templated by the BCP domains. Second, we will consider the case of anisotropic NPs within circle-forming BCP phases and the competition between the BCP mesophase and the NP nematic order. In both cases, we additionally consider 3D simulations to provide support for the conclusions drawn from the 2D simulations.

Alignment and Morphology of BCP/Anisotropic NPs
In order to study the cooperative co-assembly of anisotropic NPs and BCP morphology, we explore the overall phase behavior of Figure 1. Phase behavior of a BCP/rhomboid system with aspect ratio e = 1∕2. In (a) the phase behavior of the BCP is explored for a BCP composition f 0 and an NP volume fraction p . The phase of the BCP is marked as circles for the hexagonal circular phase, with the color indicating the majority phase, with blue for A as majority and red for B as majority phase. Black squares indicate the lamellar phase. Additionally, regions of high nematic order are marked as green asterisk. A snapshot of a system with high nematic order and NP-induced lamellar morphology is shown in (b) with a detail view in (c). The BCP is shown in gray and white, respectively for majority and minority monomer densities, while the NPs are colored according to the local nematic order S cc of each particle (see colorbar in the right). Snapshots correspond to the final step, with a dimensionless time t∕t BCP = 1.2 × 10 4 .
the system in terms of the NP volume fraction p and the BCP composition f 0 , respectively controlling the effect of NPs and the BCP morphology. Isotropic NPs have been shown to induce phase transitions in the BCP morphology due to the changes in the effective concentration of the hosting domains. [25,36,57] Nonetheless, in this work we focus on the collective alignment of NPs with the BCP matrix and with themselves. In 2D, rhomboidal particles with n = 0.6 are used. The NP shape is set to 2a∕H 0 = 0.211 and 2b∕H 0 = 0.105 (e = 1∕2), while the role of anisotropy will be explored in the following sections. Figure 1 shows the change in the BCP morphology induced by NPs. In Figure 1a, the phase of the overall system is determined via the number of BCP independent domains, displaying hexagonal circles (A-majority), lamellar and hexagonal circles (B-majority) for blue circles, black squares and red circles, respectively. The presence of a volume fraction p of NPs is shown to modify the BCP morphology. Furthermore, the colloidcolloid ordering is characterized by the nematic order parameter S cc , with green asterisks marking regions of high nematic ordering.
The BCP phase is shown to transition from hexagonally ordered A-majority circular to lamellar (blue circles to black squares) and from lamellar to B-majority circles (black squares to red circles) in the presence of NPs. While this behavior is consistent with simpler isotropic NPs [25,36,57] and homopolymer additives, [58,59] two additional distinct features can be observed. First, the lamellar phase space is enhanced in the presence of anisotropic NPs, which is due to the NP anisotropy that effective swells the BCP domains. Second, regions of high nematic ordering are achieved at high volume fractions, but only within the lamellar-phase. Rhomboidal NPs do not assemble nematically within the B-majority circular phase (red circles), despite the larger volume fraction of NPs in the system than in the lamellar phase at, for instance, f 0 ≈ 0.4. Rhomboidal NPs are better accommodated nematically in the lamellar phase, as opposed to the two circular phases. Within the lamellar phase, the high nematic regions can grow in time in the direction of the lamellar axis, which is not possible in a more isotropic circular phase. This indicates the cooperation between the BCP symmetry-elongated lamellar domains with a preferential axis-and the NP shapeelongated rhomboidal.
Furthermore, Figure 1b-with f 0 = 0.36 and p = 0.2-shows the co-assembly as the BCP serves as a template for the NP orientation and positioning: NPs are aligned along the direction of the lamellar domain and segregated within the minority domains. NPs are colored according to the local nematic ordering S i cc , with red, yellow and blue particles indicating parallelnematic, random and perpendicular-nematic ordering, respectively. While considerable nematic ordering is present (long redcolored areas), deviations from this behavior can be correlated with defects in the BCP morphology: In Figure 1b randomly oriented yellow-colored particles can be found in the high-curvature regions of the BCP domains, and at the intersection between two or more domains. In these instances, the nematic order cannot be preserved as NPs with different director vector merge. On the other hand, in Figure 1c NPs with perpendicular-nematic ordering (blue-colored) are found at domain ends. We note that the pure BCP morphology for f 0 = 0.36 is circular, and it only acquires a lamellar morphology after a NP-induced phase transition following the addition of a concentration p of NPs.
The time evolution of the co-assembly shown in Figure 2a additionally demonstrates the correlation between the nematic order and the lamellar mesophase: As the nematic order S cc grows, the curvature of the BCP interfaces decreases, which can can be seen in (b), (c), and (d), respectively for times t∕t BCP = 12.3, 122.5, and 1225.0 as intermediate stages preceding the final co-assembly shown in Figure 1b,c. Figures 1b and 2 have shown that despite the defects in the nematic ordering, uniaxial particles acquire a good alignment with lamellar-forming BCP. Contrary to that, Figure 1a top-right region shows lack of nematic ordering at same volume fraction as (b) and (c). Furthermore, in Figure 1 the lamellar phase space is enhanced compared to isotropic NPs, [57] which motivates the study of the role of NP anisotropy in the co-assembly.
The NP shape has been experimentally shown to modify the BCP morphology, [25] where the higher effective concentration of anisotropic NPs can drive the phase transition, as NRs effectively occupy a larger area. Figure 3 shows the phase behavior of the BCP melt in the presence of a volume fraction p of rhomboidal particles with aspect ratio e and fixed area per particle, consistent with Figure 1. The pure BCP morphology is hexagonally circular (f 0 = 0.35) with NPs segregating within the minority domains, marked as blue circles and determined by the domain analysis method. At moderate concentrations the BCP morphology changes drastically depending on the NP shape: highly anisotropic rhomboids promote a transition toward elongated lamellar domains, marked with black squares. On the other hand, isotropic square-shaped NPs (e ≈ 1) swell the circular domains maintaining a circular morphology while the elongated rhomboidal particles promote the connection between domains, merging into a lamellar-like morphology. In addition, anisotropic rhomboidal NPs can easily match the lamellar phase, as they tend to form high nematic order regions, as will be shown in the following section. Right-top and right-bottom snapshots show the colloid-colloid alignment for two extreme cases: more isotropic NPs display random nematic orientation (see NP coloring, which is generally random); on the other hand, anisotropic NPs display high nematic order regions (red-colored particles), following the same principles as in Figure 1b and c, that is, correlation between defects in the BCP mesophase and the nematic NP ordering. Regions of high nematic ordering S cc are marked as yellow asterisk, while regions of high fourfold tetratic ordering S 4 cc lacking twofold nematic ordering S cc are marked as white crosses. The change in BCP morphology shown in Figure 3 is consistent with the experimental results in ref. [25] figure 2 i and  figure 3 h, where an equal NP loading the BCP morphology is respectively cylindrical and lamellar, which can be compared with Figure 3 right-top and right-bottom, respectively. In conclusion, the BCP morphology and NP alignment are clearly intertwined, with nematic ordering promoting lamellar phases.
Similarly, prolate spheroids (n = 1, a = 0.158, b = c = 0.053 and e = 0.333) in 3D promote elongated domains as the NP aspect ratio e decreases. In Figure 4a the topological quantity Euler characteristic indicates a decrease in the number of spherical domains, as the concentration grows, signaling the transition from spherical ( >> 1) to a highly connected morphology with < 0. Nonetheless, for the same volume fraction, this decrease is more pronounced in anisotropic NPs (e = 0.3) than in Figure 3. BCP phase behavior depending on the NP anisotropy, from squares (e = 1) to rhomboids (e = 0.2), and a volume fraction p of particles. The BCP phase is marked as blue squares for circular morphology and black squares for lamellar phase. The colloid-colloid alignment is marked as yellow asterisk for simulations with high twofold nematic order S cc and white crosses for low twofold nematic order but high fourfold tetratic order S 4 cc . NPs are colored according to their local nematic ordering. spherical particles (e = 1). Visually, in (b) and (c) the contrast in the BCP morphology can be observed, respectively for e = 1 and e = 0.3: anisotropic NPs promote more connected and elongated BCP domains, which translates to a smaller value of . These results indicate that both in 2D and 3D, anisotropic NPs promote highly connected BCP morphologies and that this tendency is relatively universal, as the specific NP shape plays a secondary role.
These 2D and 3D simulations show that the NP anisotropy plays a crucial role in the NP assembly and furthermore, can modify the BCP morphology. The alignment of elongated NPs within lamellar-forming BCP has been shown experimentally [15,16,25] and in simulations. [26,34,60,61] Nonetheless, the NP nematic ordering in Figures 1b and Figure 3 is confined to the BCP lamellar domains, and the NP alignment is dictated by the BCP orientation, as the BCP morphology and the NP shape are both uniaxial.

Global Nematic Ordering of NPs
In this section we study the emergence of orientational/translational ordering in non-lamellar morphologies, where there is significant contrast between the NP shape and the BCP morphology. Majority-compatible NPs dispersed in a circle-forming BCP (f 0 = 0.6) are disordered within the continuous phase of the BCP at low concentrations, as shown in Figure 5 for p = 0.1 (b). The circular domains are organized in an hexagonal lattice, while the presence of NPs induces a large density of defects.
As the volume fraction of NPs is increased, the NP local concentration becomes considerable, as the interior of the BCP domains is excluded for the NPs. In order to maximize packing, NPs are driven toward BCP domains, acquiring a tangential ori-entation with respect to the BCP interface, as can be seen in the scheme in Figure 5 for a moderate concentration p = 0.25 (c). In this regime the BCP templates the orientation and location of NPs in a different way than shown in Figure 1b.
We quantify the orientational order with respect to the BCP by calculating the coupling nematic order with respect to the center of the BCP domains S cpl (r), that is, with respect to the direction shown in the scheme in Figure 5c. In function of the distance to the center of domains r, Figure 5a indicates the collapse of the curve S cpl (r) into a single curve with open symbols, given by the BCP periodicity, for several NP concentrations. This can be understood as the consequence of the BCP templating the NP orientational and positional order. Furthermore, a clear negative peak is found in the immediate vicinity of the BCP interface, where NPs acquire tangential orientation. This peak is followed by a positive one, corresponding to the presence of a first domain neighbor, thus it corresponds to the BCP periodicity.
For higher concentrations, the tendency of uniaxial particles to organise nematically dominates over the BCP mesophase. This allows for the formation of large, elongated areas of high nematic order, which can be seen as red stripes in Figure 5d. The global director of one of these stripes is shown as a white arrow. These areas of global nematic ordering are, nonetheless, incompatible with the hexagonal circular mesophase of the BCP. The stretching of the BCP periodicity is clear in (d), as the distance between BCP domains is enhanced in the direction perpendicular to the nematic stripes (yellow dashed line) and shrinks in the same direction as the global nematic director (blue dotted line). Quantitatively, Figure 5a shows the change in regime as the p = 0.45 regime deviates from the general behavior: as the BCP changes its periodicity, the secondary peaks become smaller in absolute value and in its position.
In Figure 5d we can observe the NP phase separation into denser and dilute regions, corresponding to nematically and randomly oriented regions, respectively. This is a consequence of anisotropic NPs effective packing fraction being considerably different in the random (high effective area) and nematic order (high packing).
The progressive competition between the BCP hexagonal order and the NP nematic order is best quantified in Figure 6 where the BCP domain bond order parameters Ψ 4 (fourfold) and Ψ 6 (sixfold) are shown for different volume fractions, as well as the colloid-colloid nematic order. As the concentration grows above p > 0.25, the nematic order increases, which is followed by an increase in the fourfold order. The correlation of these two quantities suggests that the system tends to self-assemble into a rectangular lattice in order to better accommodate the nematic ordering. Meanwhile, the sixfold ordering decreases continuously, which can be attributed to both the kinetic trapping induced by small concentrations of particles and the formation of large nematic areas at large concentrations, which distorts the hexagonal mesophase.
Furthermore, the time scales of the co-assembly are largely affected by the NP concentration. At moderate concentrations p = 0.25, in Figure 7b we can observe that the colloidal ordering is rapidly reached, with S cc reaching a plateau. On the other hand, the BCP mesophase occurs in a slower time scale, corresponding to the time scale in which defects in the hexagonal ordering slowly disappear as Ψ 6 grows, shown in (a). Contrary to that, in the nematic-dominated regime at higher concentrations p = 0.45 NPs acquire nematic order over a much slower time scale (b), while the BCP mesophase ordering tracked by Ψ 6 can be seen to reach a plateau in (a). This suggests that NPs at sufficiently enough concentrations can prevent the BCP mesophase from forming with a slow nematic order being developed. Additional simulations have been performed where NPs are initialized with a given global direction, which has no effect in the final global orientation of NPs. This can be attributed to the role of thermal fluctuations in the NPs and specially the effect of the BCP microphase separation.
In Figure 8 we can observe that the NP size plays an important role in the decrease of the sixfold symmetry of the BCP mesophase and the emergence of fourfold order. As the concentration grows, we can see that relatively large rhomboidal particles (e = 0.5, 2a∕H 0 = 0.84) promote a rapid increase in the fourfold order almost immediately, in contrast to smaller NPs. Since these larger NPs are comparable in length with the BCP periodicity, the presence of even a small concentration can produce significant changes in the BCP morphology, as they cannot be easily accommodated within the BCP phase. Contrary to that, Figure 5a shows that small anisotropic NPs are disordered within the continuous matrix and it is only at moderate concentrations that nematic order emerges. This can be contrasted visually, comparing   the snapshot in Figure 8b for larger NPs with Figure 5c for a similar concentration.
In 3D, the templating effect of the BCP can be observed for prolate ellipsoidal particles-2a∕H 0 = 0.84, 2b∕H 0 = 2c∕H 0 = 0.21 (e = 0.25)-within a BCC-forming BCP-f 0 = 0.7-as in Figure 9 where, in addition to the clear peak in the colloidal alignment, a density peak can be observed in the immediate vicinity of the BCP domain, which indicates the presence of a shell of particles with tangential orientation. The generality of the alignment regardless of the specific shape of the particles is a common feature in liquid crystal physics and suggests that the present tangential alignment is considerably generic, both in 2D and 3D.

Conclusions
In this work anisotropic NPs are shown to co-assemble within the BCP mesophase to produce highly ordered systems. Anisotropic NPs are shown to align along the direction of BCP lamellar domains, reproducing several experiments. Moreover, NPs can acquire local nematic order with each other, with the BCP domain axis acting as seed for the orientation. The phase diagram of such systems shows the cooperation of the lamellar BCP phase and the nematic ordering of NPs. This leads to an enhancement of the lamellar phase in the presence of asisotropic NPs due to an anisotropy-induced phase transition in circular-forming BCP melts. For the same NP concentration, the BCP does not undergo a phase transition for isotropic NPs.
Anisotropic NPs within circle-forming BCP, on the other hand, leads to the competition between the nematic ordering and the circular BCP phase. At moderate concentration the BCP templates the position and orientation of the NPs, with NPs aligning tangentially to the BCP domains. At a high volume fraction of particles, the NPs can additionally disrupt the BCP morphology and produce large NP regions of coherent orientation, which cannot be easily accommodated by the BCP morphology.
In conclusion, BCP melts can be used as templates to control the position and orientation of NPs, but the overall structure of the composite system is the result of the interplay between the colloidal and polymeric properties, which may be cooperative (nematic ordering within lamellar phases) or competitive (nematic ordering within circular phases). These results highlight the co-assembly possibilities and richness of BCP/anisotropic NPs systems.