Eﬀects of compositional asymmetry in phase behavior of ABA triblock copolymer melts from Monte Carlo simulation

. We simulate ABA triblock copolymer melts using a lattice Monte Carlo method, known as cooperative motion algorithm, probing various degrees of compositional asymmetry. Selected order-disorder transition lines are determined in terms of the segment incompatibility, quantiﬁed by product χN , and the triblock asymmetry parameters, α and β . We correlate the results of the simulation with the self-consistent ﬁeld theory and an experimental study of polyisoprene-polystyrene-polyisoprene triblock melt by Hamersky and coworkers. In particular, we conﬁrm the mean-ﬁeld prediction that for highly asymmetric triblocks the short A -block is localized in the middle of the B -domain due to an entropic advantage. This results in the middle block relaxation and is consistent with the experimental data indicating that as the relatively short A -blocks are grown into AB diblock, from the B -block side, the order-disorder transition temperature is considerably depressed.


Introduction
Triblock copolymer melts are both interesting and useful because of their potential to self-assemble into a plethora of nanostructures [1]. This ability originates from the immiscibility of covalently bonded blocks which cannot segregate on a macroscopic scale, and therefore form ordered nanophases. The nanophase self-assembly is mostly governed by a competition of chain stretching energy (which is of entropic origin) and the enthalpic interfacial energy. A linear triblock ABC copolymer chain consists, in general, of the three distinct blocks, A, B, and C, connected sequentially. Terminal blocks A and C are often built from the same type of segments, resulting in a triblock copolymer ABA which has only 2 types of segments, A and B, as the diblock. In this paper we focus on such ABA triblocks which have been subject of numerous studies. Additional incentive for studying those triblocks is related to their ability to form loops and bridges [2,3] which significantly modifies their mechanical properties. While AB-diblock copolymer melts are known to form only a few stable nanophases (such as layers, L, hexagonally packed cylinders, C, gyroid nanostructures, G, with the Ia3d symmetry, cubically packed spherical cells S, and a recently reported O 70 -phase [4,5]), the triblock melts form tens of different phases [1]. In case of a diblock AB copolymer melt, the phase behavior is controlled by the chain composition, f (volume fraction of segments of type A), a e-mail: mbanasz@amu.edu.pl (corresponding author) degree of polymerization (total number of segments), N , and the temperature-related χ parameter [6,7]. The ordered nanophase can be dissolved into a disordered phase, for example, upon heating. Phase diagrams of such melts exhibiting order-disorder transition (ODT) lines, also referred to as binodals, and order-order transition (OOT) lines are known from experiment [8] and are successfully predicted by mean-field (MF) theories [9,10], such as selfconsistent field theory (SCFT) which is based on the standard Gaussian chain model [11], or theories including fluctuations [12,13]. Because, in the MF theories, it is sufficient to know the composition, f , and the product χN in order to determine the nanophase [9,14,15], the diblock phase diagram are often mapped in (χN , f )-plane. Similarly, the MF phase behavior of the ABA triblock melt is governed by χN and the triblock composition, which can be parametrized by 2 convenient numbers, α and β, as follows: where f i (i = 1, 2 or 3) is the volume fraction of the block of type i (labels 1 and 3 correspond to the terminal A-blocks, and label 2 corresponds to the middle B-block); f 1 +f 2 +f 3 = 1. In most previous studies [16][17][18][19] a different parametrization was used, as shown below: The first parameter, α, provides a measure of asymmetry between terminal A-blocks (f 1 and f 3 ); α = 0 corresponds to a symmetric triblock with equal terminal blocks. The second parameter, β, quantifies the volume fraction of the middle block; β = 0 corresponds to volume fraction equal to 1/2.
The ABA triblock phase diagram can be mapped in the MF approximation, as shown for symmetric (α = 0) [18] and asymmetric (α = 0) [19] cases by Matsen who calculated ODT and OOT lines as a function of τ and f A for 3 selected degrees of incompatibility (χN = 20, 30 and 40). This phase diagram exhibited a continuous change from the AB diblocks to the ABA triblocks of varying asymmetry, but the binodals were deflected for some critical degree of asymmetry. This effect was accounted for by the localization of short terminal blocks in the middle-block B-domain due to a relative entropic advantage. This was demonstrated by both numerical SCFT and a simplified analytical model based on the strong segregation theory (SST) initially developed for copolymers in refs. [20,21]. This picture is supported by the experimental work of Hamersky et al. [22] in which the phase behavior of PS-PI-PS triblocks (2 series: 9-46-A 2 , and 9-17-A 2 with the A-block grown into AB diblock from the B-block side) was reported. While T ODT was expected to increase with increasing copolymer molecular weight due to greater incompatibility (expressed by χN ), it was, somewhat counter-intuitively, depressed (up to 43 • C) as the length of the grown A-block became comparable with the other A-block, and then the trend reversed giving rise to a minimum in T ODT . Hamersky et al. explained this behavior by localization of the short A-block within the B-blocks, which resulted in relaxation of the stretched B-midblock, and thus lowering T ODT . This explanation was mostly inspired by Matsen's work [19], but the calculations of theoretical MF binodals were based on work of Mayes and Olvera de la Cruz [16,17]. It should be noted that the theories developed by Mayes and Olvera de la Cruz are conceptually different from the SCFT approach (which is, in principle, exact within the MF framework), because they are based on 4th-order expansion of the free energy as introduced by Leibler [9], and also (in ref. [17]) on the oneloop approximation (OLA), which goes beyond the MF by including some fluctuations, as proposed by Fredrickson and Helfand [12]. Since the explanation of the localization of the short A-blocks within the B-domain is based on the MF theory in strong segregation limit (χN → ∞) we also intend to use Matsen's SCFT [19] to illustrate the effect of T ODT depression. Thus we consistently use theories which are on the same level of approximations.
While the SCFT is both efficient and successful in elucidating the phase behavior of block copolymer melts, it is still a mean-field approach, and therefore it offers only an approximate solutions to the Gaussian models of block copolymer melts. In particular, the OLA theory gives the following corrections to the binodals of symmetric diblocks (f = 1/2) [12] and triblocks (α = β = 0) [17]: where (χN ) MF ODT = 10.5 and κ = 41 for the diblock, and (χN ) MF ODT = 18 and κ = 90 for the triblock. Note that MF should be exact in the strong segregation limit (χN → ∞). However, there is evidence that this theory strongly underestimates the ODT fluctuation shift in χN as shown in Monte Carlo simulations [23][24][25]. Recent advances in field-theoretic simulations [7,13] also demonstrate, that inclusion of field fluctuations (beyond the self-consistent mean-field approximation) is crucial for quantification of microphase separation in these systems. Particle simulations of block copolymers pose a considerable challenge because of the finite-size effects [26,27]. In this work we use a lattice model, because lattice simulations tend to be more efficient (compared to continuum ones [28]), while capturing the essential physics of uncharged block copolymer assembly. However, there many practical problems with obtaining the exact results in Monte Carlo simulation. Finite-size effects, competition of many relevant length scales and prohibitively long relaxation times make the simulation quite challenging. In addition, there are also effects of the underlying lattice. Finally the relatively short chain in lattice MC simulation are not Gaussian, as noted for example by Groot and Madden [29] in their interpretation of the dissipative particle dynamics (DPD) simulations of diblock melts, also in simulation of Abu-Sharkh and AlSunaidi [30] of triblock melts.
In this work we wish to confirm the SCFT predictions for the ABA triblock melt by MC simulations. In particular, we intend to verify the binodals for selected paths in the phase diagram, estimate the fluctuation effects, confirm the localization of the short A-block in the middle of the B-domain, and relate this effect to experimental study of Hamersky et al. [22].

The SCFT binodals
In order to get a starting picture for the ABA triblock copolymer phase diagram, we calculate the spinodal surface where the disordered phase becomes unstable with respect to small composition fluctuations, using Random Phase Approximation (RPA) [9] presented in appendix A. While the RPA spinodals and the SCFT binodals are different, they usually (but not always) follow the same trends. The spinodals, expressed in terms of χN , are shown in fig. 1, as a function α and β. We can identify a saddle point for the symmetric triblock (α = β = 0). In this figure we show 3 lines indicating paths showing different behavior: This can be compared to the SCFT binodals as reported by Matsen [19]. Since Matsen reported binodals only for 3 values of χN (20, 30 and 40), we fit his results to a polynomial as presented in appendix B. (data for symmetric diblock and triblock are taken into account; (χN ) c = 10.5 for α = ±1/4 and β = 0, (χN ) c = 18 for α = 0 and β = 0) In fig. 2 we show the ODT surface which resembles spinodal surface in its general shape, but is shifted towards lower χN 's, as expected.

Monte Carlo simulations
The simulations are performed using cooperative motion algorithm (CMA) [31] for a face-centered cubic (FCC) lattice with the coordination number z = 12 and the bond length √ 2a, where a is the FCC lattice constant. Chain bonds are not allowed to be broken or stretched. Standard periodic boundary conditions are applied. The lattice box size is chosen to fit the chain, and the lattice sites are completely filled with chain segments -there are no vacancies. Since all lattice sites are occupied, a chain segment can move if other segments move simultaneously. An attempt to move a single segment defines a single Monte Carlo step. The interaction energy between segments of types i and j is given by ij , where AA = BB = 0, and AB = . The interaction is limited to the nearest neighbors (z = 12), and the interaction parameter, , is related to the Flory χ parameter by the following equation: Parameter serves also as an energy unit and we can define the reduced dimensionless interaction parameter * ij = ij / , the reduced energy per lattice site and the reduced temperature as where n a is the number of lattice sites. We start the simulation by equilibrating the system in the athermal limit, that is, where /(kT ) is zero. When the system reaches its thermal equilibrium, the polymer chains assume statistical conformations, random orientations, and become uniformly distributed within the simulation box. We record the translational diffusion of the copolymer chains in order to estimate the simulation time scale. In the athermal melt, it takes about 1.4 × 10 4 MC timesteps to diffuse at a distance of the order of the radius of gyration of the 8-16-8 copolymer chain. We equilibrate the athermal melt for 1 × 10 7 MC timesteps, and from the equilibrated melt state the system is quenched to a required temperature. We also verify the quality of thermal equilibration by heating the system up and cooling it down again. This procedure is fully described in [32]. For each temperature we perform 5 × 10 6 Monte Carlo (MC) timesteps. First 2 × 10 6 are to equilibrate the system, and latter 3×10 6 to sample the data. For a given T * we repeat the simulation experiment six times starting with different initial states. For a given state point, all runs yield the same type of nanostructure, and the results are averaged over all such runs. While this method works quite relatively well for high temperatures, it tends to generate long relaxation times for lower temperatures. This results in unreliable estimates of the sampled properties. To solve this problem many modifications to the standard methods were proposed, such as parallel tempering method, in which by parallel simulation of many replicas in the relevant temperature range, the energy barriers of the local free energy minima can be overcome [33,34]. In this work we do not use the parallel method.
To calculate binodals in this simulation we monitor such quantities as the energy per lattice site, E * n = E * /n a , specific heat calculated as where n c = n a /N is the number of chains, the variations of the end-to-end distance of triblock chain, R 2 , as a function of the reduced temperature. Moreover we also calculate the structure factor, S(k), by averaging over statistically independent configurations using the following equation: where n A denotes the number of segments of type A and r m denotes the position of the m-th segment of type A. The magnitude of wave vector, k, varies from k min = 2π/L to k max = 2π/b, where L is a size of the cubic lattice whereas b is the distance between nearest monomers, i.e. √ 2a. Moreover, k vectors are commensurate with the simulation box size, and this constraint limits their number and allowed lengths. Because the system may not be isotropic, S(k) is calculated by averaging over all S( k ) such that | k | is equal to k. We take the emergence of multiple peaks in S(k) as a signature of the order-disorder transition.

Simulation results and discussion
First we present the relevant results for an asymmetric triblock, 1-16-15 (consisting of one A-segment, followed by 16 B-segments and terminated with 15 A-segments), melt in order to verify the effect of the localization of the short A-block (in this case consisting only of one segment) in the middle of the B-block domain, and next we show and discuss the Monte Carlo binodals along the MAX, MON, and MIN paths, as defined in the introduction.
We start the simulation of the 1-16-15 melt, in the 32 × 32 × 32 box, at the athermal limit as described in the previous section. Next we gradually decrease the temperature to obtain the strongly segregated layers at T * = 4.0. At this temperature we calculate the density profiles of segments from different blocks, as shown in fig. 3, in the direction normal to layers. As expected from Matsen's SST approach [19], the short A-block is localized within the B-domain with the maximum density in the center of the B-domain. This is also confirmed by a visual examination of the simulation snapshots as shown in fig. 4 which  clearly show the localization of the short A-blocks. Interestingly, this effect was used by Hamersky et al. to explain their experimental results for PS-PI-PS triblocks (2 series: 9-46-A 2 , and 9-17-A 2 with the A-block grown into AB diblock from the B-block side). While T ODT was expected to increase with increasing copolymer molecular weight due to greater incompatibility (expressed by χN ), it was actually depressed (up to 43 • C) as the length of the shorter block was increased until both A-blocks became comparable, and then the trend reversed giving rise to a minimum in T ODT . Hamersky et al. explained this behavior by localization of the short A-block within the B-blocks, which resulted in relaxation of the stretched B-midblock, and thus lowering T ODT . In fig. 5 we show both the SCFT binodals and the experimental ones [22] for the 9-46-A 2 series. We verify that the T ODT depression upon growing a relatively short A-block. It is worthwhile to stress that, unlike Hamersky et al. [22] who used the approach of Mayes and Olvera de la Cruz, we used the SCFT extrapolated binodals. Thus we apply a theory which is on the same MF level of approximation as the MF theory developed for explaining the A-block localization effect [19]. For the MON case we use 20-4-8, 18-8-6, 16-12-4, 14-16-2, and 12-20-0. For the MIN case we use the simulation data from our earlier work [3]. The ODT temperatures can be estimated by monitoring the reduced energy per lattice site, E * /n a , specific heat, C V , mean-squared end-to-end distance, R 2 , as well as structure factor, S(k). For example, in fig. 6 we show results for 8-16-8 triblock copolymer melt simulated on the 32 × 32 × 32 lattice lattice. We notice that at about T * ≈ 8.0 the reported quantities show a characteristic behavior. In particular, a decrease of energy is seen upon cooling, and also an increase of mean-squared end-to-end distance, and a maximum in the specific heat. We calculate structure factor from eq. (12), as we demonstrated in ref. [35]. At high temperatures we observe a single broad peak, characteristic of a disordered phase, but as the T * is lowered the multiple and narrow peaks develop which indicates the onset of the ODT. Those peaks can be used to identify the nanophase at lower temperatures. This identification can be confirmed by visual examination (snapshots) of the simulation configurations.
In fig. 7 we show the simulation ODT lines (binodals) for the MAX path and the corresponding SCFT results. It is interesting to notice that increasing asymmetry causes a decrease of the ODT in terms of χN . It means that, if we convert χ's to temperatures, increasing asymmetry stabilizes the ordered phase, as expected from the MF theory and observed in experiment. While a qualitative agreement between the simulation and SCFT binodals can be observed, we can see that the MC binodals have a maximum which is more flat, compared to the SCFT maximum. We observe ordered layers for all simulation points. The simulation binodals are shifted in terms of χN by a factor of 2, when compared to the SCFT binodals. For α = 0 this shift is larger than the OLA prediction but smaller than the Matsen's estimate [23] which shits the (χN ) ODT by a factor of 2.7. This shift is expected by the Fredrickson-Helfand theory, but the its magnitude is not easy to estimate for a variety of reasons: -while the SCFT approach is based on the Gaussian chain model, the Monte Carlo method is based on the lattice model, and therefore the results for those models may differ, particularly for relatively short chains and blocks; -the relatively short lattice chains are not Gaussian; -field-theoretic simulation [7] are the exact way to calculate the fluctuations for the SCFT.
In fig. 8 we present the Monte Carlo binodals (a) and the corresponding SCFT binodals (b) as a function of β parameter with a fixed asymmetry parameter, α = 0.1875. The SCFT binodals for are compared with simulation results. Simulation has been performed for five β's values: 0.125, 0, −0.125, −0.250 and −0.375 (not shown), which corresponds to 12-20-0 (diblock), 14-16-2, 16-12-4, 18-8-6 and 20-4-8 microarchitectures, respectively. In fig. 8 we can notice a qualitative agreement between the SCFT calculation and simulation. As β is increased along this path, the binodals monotonically decrease. Moreover, we observe a variety of ordered phases, including L, PL, S and a an unidentified bicontinuos (B) phase. The PL phase is probably the non-equlibrium phase corresponding to the equilibrium G phase [36].
It should be noted that eq. (8), relating χ to the reduced temperature, may not approximate for comparing  the results of the lattice model and the Gaussian chain model, as demonstrated by Morse and Chung [37], and earlier by Muller and Binder [38]. Morse and Chung have shown that the effective χ parameter should be extracted from the lattice simulation as follows: -the effective number of nearest intermolecular contacts has to be calculated in the athermal state for different chain lengths, N ; -this number has to be extrapolated to N → ∞ (or 1/N → 0), yielding z ef f ; -z ef f is to be used in eq. (8) instead of (z − 2 = 10).
For this simulation we estimate z ef f ≈ 7.5, and this value should be used to calculate the effective χ paramater, χ ef f , which is smaller by 25% compared to the original value, obtained with z − 2 = 10. This rescaling improves the overall agreement between the MF and Monte Carlo results by the same factor, that is 25%.
In summary, we show that the Monte Carlo binodals follow the same trends as the SCFT binodals for the 3 representative cases, that is for MAX, MON and MIN paths. However, the SCFT binodals are significantly shifted towards higher values of the incompatibility expressed by χN , and also by χ ef f N .

Conclusions
We simulate the ABA triblock copolymer melt using the CMA method, probing various degrees of compositional asymmetry, and determine binodals in term of the segment incompatibility, quantified by product χN , and the asymmetry parameters, α and β. We correlate the results of the simulation with the self-consistent field theory and the experiments of Hamersky et al. [22]. In particular, we verify that for a highly asymmetric triblock melt the short A-block is localized in the middle of the Bdomain as predicted by the mean-field theory [19]. This is accompanied by the relaxation of the middle block. Moreover, this agrees with the experimental report of Hamersky and coworkers showing that as the A-block is grown into AB diblock, from the B-block side, the orderdisorder transition temperature is depressed for short A-blocks.
We also show that the Monte Carlo binodals follow the same trends as the SCFT binodals for the 3 representative cases, that is for MAX, MON and MIN paths. However, the SCFT binodals are significantly shifted towards higher values of the incompatibility expressed by χN , and also by χ ef f N . This shift is expected by the Fredrickson-Helfand theory, but its magnitude is not easy to estimate. The SCFT approach is based on the Gaussian chain model, the Monte Carlo method is based on the lattice model, and therefore the results for those models may be difficult to relate, particularly for relatively short chains and blocks which are not Gaussian.
Finally it is worth to notice that the general shape of the RPA spinodal and the SCFT binodals agree. Since the SCFT binodals correctly describe the lowing of the T ODT upon successive growing the A-block from the B-block of the initial AB diblock, it may be possible to capture this effect by RPA alone, without resorting to the MF theories, such as SCFT and SST, not to mention more sophisticated theories which include fluctuations.