Equilibrating high-molecular-weight symmetric and miscible polymer blends with hierarchical back-mapping

Understanding properties of polymer alloys with computer simulations frequently requires equilibration of samples comprised of microscopically described long molecules. We present the extension of an efficient hierarchical backmapping strategy, initially developed for homopolymer melts, to equilibrate high-molecular-weight binary blends. These mixtures present significant interest for practical applications and fundamental polymer physics. In our approach, the blend is coarse-grained into models representing polymers as chains of soft blobs. Each blob stands for a subchain with Nb microscopic monomers. A hierarchy of blob-based models with different resolution is obtained by varying Nb. First the model with the largest Nb is used to obtain an equilibrated blend. This configuration is sequentially fine-grained, reinserting at each step the degrees of freedom of the next in the hierarchy blob-based model. Once the blob-based description is sufficiently detailed, the microscopic monomers are reinserted. The hard excluded volume is recovered through a push-off procedure and the sample is re-equilibrated with molecular dynamics (MD), requiring relaxation on the order of the entanglement time. For the initial method development we focus on miscible blends described on microscopic level through a generic bead-spring model, which reproduces hard excluded volume, strong covalent bonds, and realistic liquid density. The blended homopolymers are symmetric with respect to molecular architecture and liquid structure. To parameterize the blob-based models and validate equilibration of backmapped samples, we obtain reference data from independent hybrid simulations combining MD and identity exchange Monte Carlo moves, taking advantage of the symmetry of the blends. The potential of the backmapping strategy is demonstrated by equilibrating blend samples with different degree of miscibility, containing 500 chains with 1000 monomers each. Equilibration is verified by comparing chain conformations and liquid structure in backmapped blends with the reference data. Possible directions for further methodological developments are discussed.

Understanding properties of polymer alloys with computer simulations frequently requires equilibration of samples comprised of microscopically described long molecules. We present the extension of an efficient hierarchical backmapping strategy, initially developed for homopolymer melts, to equilibrate high-molecular-weight binary blends. These mixtures present significant interest for practical applications and fundamental polymer physics. In our approach, the blend is coarse-grained into models representing polymers as chains of soft blobs. Each blob stands for a subchain with N b microscopic monomers. A hierarchy of blob-based models with different resolution is obtained by varying N b . First the model with the largest N b is used to obtain an equilibrated blend. This configuration is sequentially finegrained, reinserting at each step the degrees of freedom of the next in the hierarchy blob-based model. Once the blob-based description is sufficiently detailed, the microscopic monomers are reinserted. The hard excluded volume is recovered through a push-off procedure and the sample is re-equilibrated with molecular dynamics (MD), requiring relaxation on the order of the entanglement time. For the initial method development we focus on miscible blends described on microscopic level through a generic bead-spring model, which reproduces hard excluded volume, strong covalent bonds, and realistic liquid density. The blended homopolymers are symmetric with respect to molecular architecture and liquid structure. To parameterize the blob-based models and validate equilibration of backmapped samples, we obtain reference data from independent hybrid simulations combining MD and identity exchange Monte Carlo moves, taking advantage of the symmetry of the blends. The potential of the backmapping strategy is demonstrated by equilibrating blend samples with different degree of miscibility, containing 500 chains with 1000 monomers each. Equilibration is verified by comparing chain conformations and liquid structure in backmapped blends with the reference data. Possible directions for further methodological developments are discussed.
Keywords: multiscale modeling, polymer mixtures, polymer simulations, blends, soft models (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Understanding the properties of liquids comprised of long polymer molecules using computer simulations, frequently requires models reproducing key microscopic features of chain architecture and liquid structure. However, the preparation of equilibrated samples described with such models is challenging since microscopic features, such as hard excluded volume, in combination with mesoscopic chain connectivity lead to protracted relaxation times [1]. Thus developing nondynamic approaches for efficient sampling of configuration space of polymeric liquids is an area of active research. Several strategies are available, including advanced rebridging Monte Carlo (MC) algorithms [2,3], configuration assembly procedures [4][5][6][7], and hierarchical backmapping approaches [8][9][10][11][12][13][14][15]. This classification is not strict, e.g. it is possible to combine concepts from configuration assembly and hierarchical backmapping [16,17].
Both configuration assembly and hierarchical backmapping are divide-and-conquer strategies. In configuration assembly procedures the polymer liquid is first described with crude resolution, aiming to reproduce equilibrium longwavelength structural and conformational properties. This is typically achieved by assembling an ideal gas of chains, reproducing conformational distributions expected [4][5][6][7] in the polymer liquid. Correlations are then partially recovered, e.g. by reducing density fluctuations through a MC optimization [4,5,7]. Subsequently, the microscopic degrees of freedom are reinserted into the long-wavelength 'matrix' by recovering slowly the microscopic excluded volume (pushoff process). Because at this stage only short-wavelength properties must be equilibrated, the computations remain feasible. Although the computation time for optimizing the long-wavelength liquid structure increases with chain length, currently such techniques can equilibrate large samples of entangled chains [7]. Nevertheless, because the long-wavelength description results from an algorithmic construction rather than a strict statistical-mechanical treatment of a coarse-grained (CG) model, configuration assembly methods encounter some fundamental difficulties. In principle, generating conformations assuming that long chains in a polymer melt are ideal is an approximation [18]. An a priori prediction of conformational properties in inhomogeneous systems is even more challenging [5].
In contrast to configurational assembly, hierarchical backmapping strategies are schemes linking simulations performed in a standard statistical-mechanical framework. The polymeric material is described on different levels of resolution, where the microscopic representation occupies the last level. On the remaining levels, the material is described through suitable CG models. The sample is first equilibrated using the coarsest model and one proceeds by reinserting the degrees of freedom of the next level. The system is re-equilibrated and the procedure is repeated until the target microscopic description is reached.
The efficiency of backmapping schemes can be tremendously increased [11][12][13][14][15][16][17] by incorporating into the hierarchy mesoscopic models with soft potentials, e.g. comparable in strength to the thermal energy kT. For example, equilibrating intermediate and long-wavelength properties of a polyethylene melt using a model with a soft potential (inspired by dissipative particle dynamics), prior to switching to Lennard-Jones interactions, allowed [17] the preparation of samples containing 10 3 highly entangled chains (10 3 united atoms). Melts with 10 3 bead-spring (BS) chains with about 200 entanglements per chain were equilibrated using a multiscale approach including a lattice soft model [14]. Hierarchies based on models describing polymers as strings of soft blobs [19][20][21][22][23] are straightforward to construct [12,24], because the sequence of different resolutions is naturally generated by adjusting the amount of microscopic monomers lumped into each blob. A proof-of-concept study using sequential fine-graining within blob hierarchies reported [12] the generation of melts containing, e.g. 10 3 chains with 2 × 10 3 beads. The method is not limited to the sample-sizes reported in that study, since the efficiency of the scheme depends on system volume but not on chain length [12].
In contrast to homopolymer melts, hierarchical backmapping based on soft models has been rarely employed to generate microscopically described samples of multicomponent systems with long polymer chains [15]. Such samples are required, for instance, as starting configurations for studying rheological properties of highly entangled miscible blends. Understanding rheology of miscible blends presents tremendous interest for materials development, motivated by the remarkable improvement of toughness in multiple network gels and elastomers [25][26][27]. Miscible polymer blends could serve as a basis for studying these improvements, utilizing their different reactivities and solubilities against cross-linking chemical agents as in [28,29]. Comparing to single component polymer liquids, developing hierarchical backmapping schemes for multicomponent systems is significantly more challenging. Partially, this stems from the fact that the properties that must be reproduced on the mesoscopic levels of the hierarchy are more complicated and difficult to quantify. The degree of mutual compatibility of different components is one of these properties.
Motivated by the current state of art, we develop here a simple hierarchical backmapping scheme based on blob models, which enables the equilibration of highly entangled miscible blends. This scheme is evolved from a backmapping procedure developed earlier for homopolymer melts [12]. We consider blends of generic A and B homopolymers, described on the microscopic level through a standard BS model [30]. Though generic, this model captures key features of actual polymeric liquids including hard excluded volume, strong covalent bonds, and realistic monomer density. Only symmetric systems are addressed here, where the two comp onents are differentiated exclusively via interactions between A and B monomers. For the mesoscopic representation of these blends we adapt the simple blob-based model [23] which has been employed in the hierarchical backmapping method for homopolymers [12]. Even for this simple multicomponent system, the development of a hierarchical backmapping strategy requires the consideration of several important methodological questions including (a) the parameterization of the blob-based model, (b) fine-graining procedure, and (c) verification of equilibration of backmapped samples.
Our manuscript is organized as follows. Section 2, introduces the microscopic and the soft blob models used to describe the symmetric blends. In section 3, we address the question of parameterizing the soft blob models. The procedure of the backmapping is explained in section 4. An application of the backmapping procedure is presented in section 5. We provide some concluding remarks and an outlook in section 6.

Models
The workflow of the backmapping scheme used to equilibrate the blends in this work, presents a three-level hierarchy as illustrated in figure 1. First, a blend described with a crude blob-based model, figure 1(a), is equilibrated. The sample is fine-grained to obtain a configuration described by a blobbased model with doubled resolution at each step, figure 1(b). The last configuration serves as an input for the reinsertion of microscopic details as demonstrated in figure 1(c). Below we discuss the details of the microscopic and blob-based models incorporated into the backmapping scheme.

Microscopic model
We adopt as target microscopic systems for our back-mapping method, blends described by a standard BS model [30]. Specifically, we consider blends of A and B homopolymers with the same chemical structure, described by chains of beads linked by finite extensible nonlinear (FENE) springs. The potential corresponding to a FENE spring of length r is given by: There are no angular potentials between the springs, so that our microscopic model corresponds to setting the stiffness of the angular potential in [7,12] to k θ = 0. In total there are n molecules in the blend (sum of A and B homopolymers) and each chain, irrespective of its chemical species, contains N BS beads.
Non-bonded interactions are present between every pair of beads and are defined through the Weeks-Chandler-Andersen (WCA) potential: (2) Here r is the distance between interacting beads, while α and β specify their type (A or B). We set σ 0 and AA = BB = 0 as the length and energy units. K = 30 0 /σ 2 0 and R 0 = 1.5σ 0 are used for the FENE potential. The number density of beads is fixed at 0.85/σ 3 0 . The time scale of the model is given by τ = mσ 2 0 / 0 , where the mass of the bead, m, is set to unity. The strength of interactions between beads of different type is set by AB , which acts as a free parameter.
With this choice of parameters, a single-component homopolymer melt described by the microscopic model is characterized by an entanglement length of N e ≈ 87 monomers [7]. In the following, this estimation facilitates the choice of simulation protocols treating the microscopic model and setting the resolution of blob-based descriptions.

Soft blob-based model
As a coarse-grained model for back-mapping, we adopt the soft-blob model studied in [12,23,24,32]. Since the blobbased model has been described in detail elsewhere [23,32] only a brief summary is presented here. Specifically, the microscopically described polymer molecules are coarse-grained into chains containing N SB = N BS /N b blobs. N b stands for the number of microscopic beads in a subchain represented by each blob. The radius, σ, and the coordinates of the center, R, of the blob match the instantaneous radius of gyration and position of the center-of-mass (COM) of the underlying subchain. The potential associated with the degree of freedom σ is defined as: The connectivity of the chains of blobs is described by bond and angular potentials defined as: For a single component system (homopolymer melt) the energy of non-bonded interactions between blobs is given by: Figure 1. Concept of hierarchical backmapping using mesoscopic blob-based models of polymer chains. The configuration of a binary blend described with a crude blob-based model is presented in panel (a). Each blob in this sample is replaced by a pair of smaller ones, to obtain a configuration described with higher resolution panel (b). The fine-grained configuration serves as a framework for introducing microscopic details panel (c). The configurations have been created using the VMD software [31].
Here ρ expresses [23] the density distribution of microscopic monomers underlying a blob. This distribution is assumed to be Gaussian and is defined as: The second term in the right-hand side of equation (6) describes the blob 'self-energy' due to interactions between underlying monomers. The self-energy does not depend on R and can be seen as an additional potential, coupled to the size of the blob. The determination of the parameters a 1 , a 2 , a 3 , b, k, and ˜ 0 will be discussed in section 3.
To describe symmetric binary blends, we simply generalize the definition of non-bonded energy to: The quantity ˜ AB presents an additional parameter, controlling the compatibility of blend components.
To describe the blends, as illustrated in figures 1(a) and (b), we employ two soft blob-based models, SB25 and SB50, with resolutions set by N b = 25 and 50, respectively. All simulations involving these blob-based models employ an efficient particle-to-mesh MC method elaborated earlier [32]. The choice of the two resolutions is semi-empirical. Previous studies [12] demonstrated that N b = 25 is 'sufficiently large' for the simple blob-based description to be a valid coarsegrained representation of the generic microscopic model. Specifically, SB25 reproduces reasonably well [12] conformational and structural properties of microscopically-described homopolymer melts on scales comparable or larger than the size of SB25 blobs. In parallel, SB25 fulfils the condition N b < N e , which ensures [12] that the blend after the reinsertion of microscopic monomers must be re-equilibrated only on scales smaller than the tube diameter. This regime is governed by fast Rouse dynamics [1] and the computations are, therefore, feasible. One extra level is added to the hierarchythe SB50 model-to speed-up equilibration of SB25 blends (through fine-graining of SB50 configurations). For our study using a hierarchy with only two blob-based models is sufficient. If required, when modeling other systems, it is straightforward to include into the hierarchy models based on larger blobs [12].

Transferability assumption
The transferability ansatz plays a key role in our parameterization of the blob-based model. Specifically, we assume that the parameters defining the molecular architecture, a 1 , a 2 , b, and k, self-interactions, a 3 , and non-bonded interactions between blobs of the same species, ˜ 0 , are independent from the blend composition. Thus we set these parameters to values already known from an earlier study [12] of single-component melts of polymers, which are identical to our A (B) homopolymers. Table 1 recapitulates the parameters used in this work for SB25 and SB50 models, taken from the Supplementary Information of [12] The transferability ansatz can be rationalized taking into account theoretical predictions concerning the dependence of chain size on the strength of incompatibility between different components in miscible blends [33,34]. This incompatibility is quantified via χN , where χ is the Flory-Huggins (FH) parameter and N stands for the number of repeat units per chain. Simple scaling arguments and more elaborated calculations predict [33] that polymer chains should contract in the homogeneous phase, i.e. in the miscible regime, with increasing χN . Within the one-loop correction to Flory-Huggins mean-field theory, [34] the difference between the mean square radii of gyration of polymers in a homopolymer melt, R 2 g(m) , and miscible blend, stands for the radius of gyration of the chain if it were an unperturbed (ideal) coil. The quantity N = (ρ c R 3 e(m) ) 2 ∼ N is the invariant degree of polymerization (ρ c and R e(m) are, respectively, the number density of chains and their root mean-square endto-end distance in the melt). The coefficient C depends on χN . The maximum incompatibility for the symmetric blends considered in our study is χN = 1.7, so that C −0.4, as estimated with the help of figure 6 from [34]. For a blend of chains with N BS = 1000, which corresponds to N 1/2 ≈ 60.1, to be below 1%. Due to the weak dependence of chain conformations of long polymers on χN in miscible blends, it is plausible to assume that at least those parameters that are directly linked to coarse-grained bonded potentials (a 1 , a 2 , b, and k) can be estimated from the limit χN = 0.

Parameterizing interactions between blobs of different species
Within the transferability ansatz, ˜ AB is the only parameter that must be determined on the basis of reference configurations of blends described with microscopic detail. Because the coarse-graining procedure conserves the amount of polymer molecules, χN is an invariant quantity [35][36][37]. Therefore ˜ AB must be chosen such that the mesoscopic description of the blend conserves the χN of the microscopic blend to be generated via backmapping.
In blends described with the microscopic model from section 2.1, χN is proportional [38] to ( AB − 0 )N BS . To obtain the specific relationship between χN and ( AB − 0 )N BS , we initially prepare a set of symmetric mixtures of A and B homopolymers composed of chains with lengths N BS = 50, 100, 200 and 300. We only focus on cases where ( AB / 0 − 1)N BS 3.0 in order to remain in a miscible state [38]. Due to moderate chain lengths, these blends can be equilibrated through brute force molecular dynamics (MD) simulations. All MD simulations used the Langevin thermostat and were performed with the efficient ESPResSo ++ software package [39]. When designing the MD runs, we took into account that the Rouse relaxation time τ e of the entanglement strand N e ≈ 87 (see section 2.1) in a homopolymer melt described by the microscopic model is about 10 000τ . With this estimate in mind, the length of the MD simulations for the unentangled blends with N BS = 50 and 100 was set to about 5τ e , which is longer than the Rouse relaxation time of the chains. For the weakly entangled systems with N BS = 200 and N BS = 300, the run time was set (respectively) to 60τ e and 75τ e , which exceeds the reptation time of the chains. For all cases, we confirmed that the correlation u start u end of the end-to-end unit vectors of the chains between the initial configuration and the last configuration has decayed to zero. The chains in the samples amount to 400 for N BS = 50 and 100, 200 for N BS = 200, and 100 for N BS = 300.
After the initial samples of blends with symmetric composition are relaxed using brute force MD, they are subjected to hybrid MD/MC simulations in the semi-grand canonical ensemble [33,37]. The MC part of the scheme consists of standard identity-changing moves attempting to switch the chemical species of a randomly chosen chain. The proposed move is subjected to a Metropolis criterion under a given exchange chemical potential µ ex . After attempting n/2 type exchanges, the configuration is subjected to a short MD simulation corresponding to about 1600τ . This cycle of MC moves and MD runs is repeated 300 times for N BS = 50 and 100, 400 times for N BS = 200 and 300. From the last two third of these cycles, we estimate the average volume fraction of the A homopolymer, φ.
The mean-field Flory-Huggins theory predicts that µ ex , φ, and χN , are connected as: Overall, we observe that for the blends considered in our study the relationship described by equation (13) is fulfilled rather well. As an illustration, figure 2(a), presents an example for blends with N BS = 100, considering several values of AB . For different chain lengths, χN can be estimated [33,37] from the slope of the linear plots, as those presented in figure 2(a). Figure 2(b) presents χN extracted in this way, as a function of ( AB − 0 )N BS . It can be observed that the slope of the linear relationship between χN and ( AB − 0 )N BS changes slightly with chain length, e.g. the average slopes for N BS = 50 and N BS = 300 differ by about 8%. Such chain-length effects have been extensively investigated in previous theoretical studies [34,37,[40][41][42]. It is currently accepted that within a first order approximation the dependence of χN on chain length can be expressed as: χN ≈ χ e N 1 + C 0N −1/2 . As has been discussed previously, N denotes the invariant degree of polymerization in the melt, while C 0 is a constant. For structurally symmetric blends theoretical arguments based on perturbation theory predict [41] that this constant is universal and equal to C 0 2.64. For detailed discussions the interested reader is referred [34,37,[40][41][42] to several reviews and original studies. We mention here briefly, that this expression for χN includes corrections due to short-range fluctuation effects and correlations (e.g. local liquid packing), as well as long-range fluctuation effects. Short-range effects are included [42] in the 'effective' χ e N . By taking these effects into account, χ e N differs [37,42] from the 'bare' χ o N that one would obtain within a simple mean-field approximation. Effects of longrange fluctuations are captured by the C 0N −1/2 term. To distinguish conceptually [42] χN from χ e N , it is useful to name the former as 'apparent' [42].
In this work, we establish the apparent χN for the needs of our backmapping procedure and do not focus on the physics of its dependence on various system parameters. For our purposes, it is more convenient to approximate χN in the microscopic model as: where C 1 and C 2 are unknown parameters. Based on equation (14), the data from figure 2(b) are re-plotted in figure 3(a) as a χN/( AB − 0 )N BS versus 1/ √ N BS graph (solid circle symbols). It can be seen that the data points obtained for the same chain length but different AB scatter. We have observed no systematic dependence on AB concerning the location of the scattered points with respect to each other. Hence for each 1/ √ N BS , we consider a single χN/( AB − 0 )N BS (red crosses) obtained as an average over the corresponding scattered points. The standard deviation of these averages sets the width of the error-bars in the linear fit used to correlate these data (dashed line). From this fit we obtain C 1 = 0.56 and C 2 = 1.05. Interestingly, if we recast C 2 / √ N BS in terms of the invariant degree of polymerization N , we conclude that the calculated C 2 is equivalent to C 2 ≈ 2.0 which is not too far from the theoretically predicted universal value [41].
Using simple algebraic manipulations, it is instructive to transform equation (8) into: We have defined: as the blob volume fractions. Here ρ SB = nN SB /V is the number density of blobs. Taking into account equation (15), it is straightforward to write the energy per blob in a homogeneous phase with composition φ A and φ B , within a mean-field approximation where the instantaneous variables φ A,B are replaced with their mean values φ A,B : The last term in equation (17) is written as a sum of independent contributions of each blob and can be included into a Hamiltonian with the bonded potentials of the soft blob systems. Following the simple mean-field Flory-Huggins approach, we assume that chain conformations are not affected by mixing and neglect fluctuations of blob size. Within these approximations, we obtain a rudimentary mean-field expression for the free energy F in the homogeneous state: where the first term is the translational entropy of the chains. Note that we have set k B T = 1. As in other soft models of blends [37], the last term in equation (18) is the only contribution to the mean-field free energy which depends on the composition and will be, therefore, relevant for the phase behavior [37]. The simple mean-field free energy suggests that the apparent χN in our mesoscopic blends is proportional As in the case of the microscopic model, we identify χN for the soft blob models SB25 and SB50 via semi-grand canonical MC (casted in terms of the particle-to-mesh algorithm [32]). Figures 2(c) and (e) (SB25 and SB50 systems, respectively) demonstrate that the relationship given by equation (13) is fulfilled for these two mesoscopic blends as well. The dependence of χN on ρ SB (˜ AB −˜ 0 )N SB for chains with number of blobs ranging from N SB = 10 to 80 is presented in figures 2(d) and (f) (SB25 and SB50 systems, respectively). Overall, the plots in figures 2(d) and (f) are consistent with the expected proportionality of χN on ρ SB (˜ AB −˜ 0 )N SB . Comparing to the microscopic description, the dependence on chain length is somewhat weaker, especially for the SB50 system.
To obtain an analytical expression for the dependence of the apparent χN on chain length for the blob-based models we use the expression: Based on equation (19), the data from figures 2(d) and (f) are presented in figures 2(b) and (c) (SB25 and SB50 systems, respectively) as a χN/ρ SB (˜ AB −˜ 0 )N SB versus 1/ √ N SB plot. For each value of 1/ √ N SB , the scattered points are treated as in the case of the microscopic model. From linear regression of averaged points (crosses) we obtain for SB25 C 1 = 0.67 and C 2 = 0.26, while for SB50 we obtain C 1 = 0.78 and C 2 = 0.13. Once the parameters C 1 and C 2 in equations (14) and (19) are known, matching χN in the microscopic and a blob-based model is straightforward. For a chosen χN , equation (14) is used to estimate the required AB in the microscopic model. Substituting the same χN into equation (19) and solving with respect to ˜ AB delivers the required strength of interactions between A and B blobs.

Back-mapping procedure
After the soft blob models are parameterized, large samples of blends with long chains can be equilibrated using a slightly modified version of the hierarchical backmapping procedure developed in [12]. The approach can be subdivided into two major stages (see section 2 and figure 1) presenting: (a) equilibration of a blend described with the SB50 model and subsequent fine-graining into a blend described with the SB25 model and (b) reinsertion of microscopic details into the mesoscopic configuration described by the SB25 model and equilibration of reinserted microscopic degrees of freedom.
Stage (a) was performed as described in [12]. We summarize briefly that the blend described with the SB50 model is equilibrated with the particle-to-mesh MC. Subsequently, every single SB50 blob is substituted by a pair of blobs of the SB25 model. For this substitution, the COM and the gyration radius of the two SB25 blobs match the position and the size of the coarser blob, respectively. After reinsertion, the SB25 blend is equilibrated by a short particle-to-mesh MC simulation.
In stage (b) each of the n chains in the SB25 blend is replaced by a molecule described with the microscopic model. This substitution is performed [12] so that the COM and the radius of gyration (squared) of every subschain with N b microscopic monomers match the COM and the radius (squared) of the corresponding blob. In practice, this task is accomplished by manipulating [12] the conformation of each reinserted chain with special potentials. We emphasize that during this first reinsertion of microscopic details, all intermolecular and intramolecular (apart from 1-2) non-bonded interactions are deactivated. Thus the blend presents an ensemble of independent chains placed in external fields (the potentials used to manipulate the chain conformations) and the reinserted monomers can overlap with each other.
Comparing to [12], in our study the ensemble of independent chains obtained after the first reinsertion of microscopic details is subjected to an additional optimization of chain conformations. To motivate this step, let us consider for simplicity a homopolymer melt (equivalent of a blend with χN = 0) with N BS = 1000, which has been already equilibrated in [12]. For this system figure 4 presents the internal distance plot, which constitutes one of the most sensitive conformational descriptors [4]. This quantity is obtained by plotting R 2 (s) versus s, where R 2 (s) is the mean square distance of monomers in the same chain and s is the difference of their ranking numbers along chain contour. We consider first the internal distance plots at two stages (green and black line, respectively): (a) just after the first reinsertion of microscopic details and b) after recovering all non-bonded interactions and final equilibration [12]. In the latter case, the plot (black line) is indicative of a melt equilibrated at all scales [12] and is in excellent agreement with a reference curve calculated from melts equilibrated through configuration assembly [7] (red line). At the same time, it can be observed (green line) that just after reinsertion the conformations are distorted at short scales (s 100). For scales smaller than the blob size (s 25) this distortion is linked to the fact that the conformations of short subschains N b = 25 during reinsertion are constrained to reproduce only two observables (COM location and radius of gyration. For s > 25 the conformational distortions stem from imperfections [12] of the simple blob model. The additional optimization step mitigates conformational distortions, facilitating the final recovery of the liquid structure (see below). The optimization presents a very short MD run in the ensemble of independent chains, where the auxiliary potentials controlling the COM and radius of gyration of the subchains are now removed. Instead, auxiliary angular, V bend (cos θ), and non-bonded, V s (r), potentials acting between intramolecular beads other than 1-2 neighbors are introduced. These interactions are given by: and The parameters take the values k θ = 0.912 0 and A = 3.0 0 . The value of k θ is chosen such that the mean square endto-end distance, R 2 e of an isolated chain subjected to the two auxiliary potentials matches the R 2 e in a homopolymer melt. For this purpose we benefit from the approximate relationship [4]: . (23) Here the bond length is set to l = 0.97σ 0 . Figure 4 presents the internal distance plot calculated after the optimization step (blue line), demonstrating substantial reduction in distortion of chain conformations.
The missing non-bonded interactions are introduced into the optimized samples following a standard push-off procedure [4,7,12]. The WCA potential in equation (2) is 'force capped' and the configuration is subjected to MD. The maximum allowed force is gradually increased [7], as monomer-monomer overlaps are becoming less severe, to recover finally the original WCA interactions. In the microscopic model, non-bonded interactions are differentiated depending on whether the interacting monomers are of the same or different chemically species. In general, one could adjust the protocol for removing the force-capping during the push-off to the type of the potential. Here we use the simplest approach, where all WCA potentials are subjected to the same force-capping. We will return to this point in section 5, when presenting an application of backmapping. After the push-off, the blend is subjected to a standard MD simulation which relaxes locally the reinserted microscopic degrees of freedom on length-scales corresponding to a couple of blobs.

Application example
Here we apply our backmapping procedure to generate samples of blends containing n = 500 chains with N BS = 1000 beads. The composition is set to φ = 0.5. Different degrees of miscibility are considered, corresponding to χN = 0.57, 1.1, and 1.7 ( AB = 1.001, 1.002, and 1.003).
As with any new computational method, it is desirable to compare the results of the hierarchical backmapping scheme with reference data obtained from established simulation techniques. The backmapping procedure was previously verified [12,24] for homopolymer melts, using reference data generated via well-established configuration-assembly procedures [4,7]. Currently, such techniques are not available for multicomponent polymer systems. Therefore, to generate reference data for blends with N BS = 1000 we rely on a variant of the hybrid MD/MC simulations, introduced in section 3). For these simulations we use as initial configurations equilibrated samples of single-component homopoly mer melts generated in [12], containing n = 1000 chains with N BS = 1000 monomers. The homopolymer melt configurations are converted into symmetric blends, assigning the polymers randomly A or B identity. These initial blend configurations are subjected to a hybrid MD/MC procedure. Specifically, we employ MC moves attempting to exchange identities between two randomly chosen A and B polymers, setting µ ex = 0. Every n/2 attempted identity-exchange MC moves, the blend configuration is subjected to MD run lasting 500τ . In total, we perform 120 such MD/MC circles, after which the samples are subjected to a MD simulation for about 50 000 τ which is several times longer than τ e in single-component melts. All MD steps are carried out using the efficient software ESPResSo ++ [39].
Formally, the starting and the final configurations of a N BS = 1000 blend prepared using the hybrid MD/MC simulations are not fully de-correlated. For example, the orientation de-correlation function of polymers calculated from these configurations has decayed to u start u end 0.88. Nevertheless, we argue that these configurations are still useful for checking the backmapping procedure. Indeed, if miscible blends generated by the hybrid MD/MC and hierarchical backmapping (which are different techniques) have the same properties, then this strongly supports the validity of both methods. Figures 5(a)-(c) present the internal distance plots calculated, respectively, in backmapped blends (red lines) with χN = 0.57, 1.1, and 1.7. The plots are compared with their counterpart (black symbols) calculated in samples prepared using the hybrid MD/MC. In all cases, the deviations between data obtained using the two methods is below 2%. The internal distance plots for the backmapped blends from figures 5(a)-(c) are re-plotted in figure 5(d) together with the internal distance plot for χN = 0. The tails of these plots are shown magnified in the inset. The deviations of the internal distance plots from each other are essentially within the statistical noise of the data and we do not observe any systematic dependence on χN . This behavior agrees with the theoretical predictions (see section 3) that in the region of studied χN , for N BS = 1000, the chain dimensions should change by less than 1%.
To verify the equilibration of the liquid structure of the backmapped blend we compare in figures 6(a)-(c) the inter-molecular radial distribution functions calculated in the backmapped blends (lines) for the same, g AA (r) and different, g AB (r), types of beads with their counterparts (lines) in the samples from hybrid MD/MC. On small scales around r ≈ σ 0 , the deviation of the radial distribution functions is at most 2.5%. As one moves to larger distances the deviation becomes even smaller and drops below 2%.
We have observed that using longer MD runs to equilibrate the reinserted microscopic beads in backmapped blends, reduces the differences between the radial distribution functions calculated in these systems and their counterparts in hybrid MD/MC. This trend of convergence supports our expectation that the blends prepared through hybrid MD/MC simulations can serve as reliable reference systems.
Overall, we observe in blends that the equilibration of microscopic degrees of freedom reinserted through the standard push-off procedure developed for homopolymer melts [4,7] requires run-times on the order of 2.8 τ e in order to reproduce the long-wavelength structure with deviations smaller than 2%. We have verified that the moderate increase of equilibration times in blends, comparing to homopolymers, should not be attributed to imperfections of the simple blobbased model describing the mesoscopic samples (used as an input for the backmapping procedure). For this verification we considered configurations of blends described with microscopic detail, prepared using the hybrid MD/MC approach. These configurations were coarse-grained into a blob-based representation, replacing each subchain with N b monomers by a blob placed at the COM of the subchain. The radius of each blob was set equal to the radius of gyration of the underlying subchain. In this way, we obtained an 'exact' blob-based representation of the miscible blend, which e.g. incorporates all complex effects of multibody interactions. These effects, in part, are not captured by the simple blob model defined in section 2.2). The 'exact' blob-based description was used as a starting configuration for the reinsertion of microscopic details following the algorithm described in section 4. The equilibration times required for g AA (r) and g AB (r) to converge to the reference data, were similar to the case of reinsertions starting from configurations described with the simple blobbased model. As expected, the local structure depends on the reinsertion procedure and prolonged equilibration times are due to the simplicity of the push-off algorithm. In all cases we observe that g AA (r) and g AB (r) converge to the reference data from above and below, respectively. These trends demonstrate that the simple push-off generates a blend with microscopic structure representative of a χN which is marginally higher comparing to the target value. To increase the efficiency of the backmapping scheme, blends require more advanced procedures to control the degree of miscibility during push-off. Improving such details will be a subject of future work.
The backmapped blends were subjected to hybrid MD/MC simulations in the semi-grand canonical ensemble (described in section 3.2) to calculate the χN parameter and demonstrate that it indeed matches the targeted value. Figure 7 presents (crosses) the calculated χN as a function of the targeted χN demonstrating that they match closely. We add the χN estimated for the backmapped blends as a χN/( AB − 0 )N BS versus 1/ √ N BS plot to figure 3 (green crosses). The results are consistent with the linear fit (black dashed line) calculated in section 3.2 for blends with smaller N BS . It is instructive to consider the values of AB corresponding to the χN parameters realized for the backmapped blends and estimate the solubility from the one-fluid approximation [37,38,40]: where g inter (r) is the inter-molecular radial distribution function for the case of χN = 0 (homopolymer melt). The (χN) rdf estimated via equation (24) is also presented in figure 7 (open symbols) and agrees well with the χN obtained from the semi-grand canonical simulations.
At a first glance, the backmapping and the hybrid MD/ MC procedures required comparable resources to equilibrate the N BS = 1000 blends. To equilibrate mixtures with n = 500 chains, the reinsertion process (push-off and final  re-equilibration) in the backmapping procedure required 72 h on 16 processors (2.6 GHz). The time for generating the blobbased configuration into which the reinsertion was performed, required only 5 h on a single processor and is negligible. To equilibrate systems with n = 1000 chains, the hybrid MD/ MC required 120 h on 64 processors of comparable architecture. These CPU times are equivalent to 2.3 × 10 −3 and 7.6 × 10 −3 CPU h / particle, respectively. But, in fact, the hybrid MD/MC method requires as an input equilibrated configurations obtained from hierarchical backmapping of long homopolymer melts. Therefore, the hybrid MD/MC is not a stand-alone equilibration technique. Moreover, comparing to hierarchical backmapping, this procedure is significantly less versatile, e.g. the identity-exchange move cannot be applied to blends with different lengths of chains or strong asymmetries in molecular architecture [43]. So far, hierarchical backmapping is the most optimal way for preparing configurations of highly entangled polymer blends with a large number of molecules. Because the equilibration time in hierarchical backmapping does not depend on chain length (it is on the order of τ e ), the method can be advantageous even for equilibrating weakly entangled blends, where N BS equals a few N e . In contrast, in brute force MD the CPU time is proportional to the reptation time τ rep , which describes the relaxation of the entire chain. Because [1] 3 , a small increase of chain length in the weakly entangled regime (see runs with N BS = 200 and N BS = 300 in our case) has tremendous effects on the efficiency of the brute force MD. At the same time, backmapping procedures for relatively short chains might be more sensitive to approximations in simple blob-based descriptions, e.g. chain-end effects increase and differentiation of blobs according to their position in the chain might be necessary.

Concluding remarks and outlook
We extended an efficient hierarchical backmapping scheme, developed earlier for homopolymer melts [12,24], to polymer blends described with a generic microscopic BS model [30,38]. For a first method development, we focused on the simplest case of symmetric miscible blends. The idea of the method can be summarized as follows. Blends are equilibrated using a model where polymers are described as chains of soft blobs. Each blob represents a subchain with a large number of microscopic monomers, i.e. the model is very coarse. The equilibrated coarse-grained blend is sequentially fine-grained, reinserting at each step the degrees of freedom of a finer blobbased model. Once the blob-based description is sufficiently detailed, the microscopic monomers are reinserted. Namely, the hard excluded volume is recovered through a push-off procedure [4,7] and the blend is re-equilibrated with MD simulations corresponding to time-scales on the order of the entanglement time, τ e .
The potential of the backmapping strategy was demonstrated by equilibrating samples of blends containing 500 chains with 10 3 microscopic monomers. For the flexible BS model used our study these degrees of polymerization lead, on the average, to about ten entanglements per chain [7,12]. For the same length of chains, the number of entanglements can be easily increased by making the BS chain stiffer [7]. The generated blends were characterized by various degrees of miscibility, defined by χN = 0.57, 1.1 and 1.7.
The investigated systems are among the simplest examples of blends. Nevertheless, their modeling required the consideration of several methodological issues. The most important are (a) development and parameterization of blobbased models for multicomponent systems, (b) reinsertion of microscopic details into configurations described by blobbased models, and (c) validation of equilibration of samples prepared through hierarchical backmapping. Although these questions were successfully addressed for the purposes of our work, there is substantial room for further methodological developments.
Our blob-based descriptions assume that the ingredients of the model (potentials and parameters) associated with the blended A and B homopolymers are the same as in singlecomponent melts of these homopolymers. Only the interactions between A and B blobs were introduced such that the blob-based model reproduces a desired χN . The potentials in our blob-based description are inspired by generic laws of polymer physics [23,24,44] and can be seen as simple approximations of the many-body potential of mean force [45]. For instance, blob-blob mean-force potentials should exhibit long tails [22,[46][47][48][49] changing at large distances from repulsive to attractive. These features are crucial [20] for reproducing, e.g. compressibility, but are missing from the repulsive Gaussian potentials in our model. Our implementation of the simple blob-based model is justified by the observation that reinsertion of microscopic details mitigates imperfections in the description of local conformational properties and liquid structure on blob level (at least within the statistical accuracy of the simulation data). Nevertheless, we expect that studying more complex blends, e.g. with strong conformational and structural asymmetries, will require more advanced blob-based models. The properties of such asymmetric mixtures depend on composition in a non-trivial way [43] making the implementation of the transferability ansatz problematic. Earlier studies [46,50,51] in solutions and melts demonstrated that blob-based models can be developed through systematic coarse-graining techniques. Integral equation theories are a promising approach [52,53,54] for developing blob-based models, without requiring reference data from computationally demanding brute force simulations.
Reinsertion of microscopic monomers into a blend described on the blob level is an important step during backmapping. We followed the standard push-off procedure [4,7], initially developed for homopolymer melts, using the same protocol to recover the hard excluded volume interactions for all types of monomers. Due to the simplicity of the pushoff scheme, the equilibration of the reinserted microscopic degrees of freedom required somewhat longer relaxation times comparing to homopolymer melts. In the future, the efficiency of the final equilibration step can be increased by developing reinsertion protocols tailored for multicomponent systems.
Validating equilibration of backmapped blends required reference samples with long chains. Since configurationassembly procedures are not available for blends, the reference samples were prepared through a hybrid MD/semi-grand canonical MC technique. The latter is not a stand-alone technique and we used as starting configurations blends created by random labeling of chains in homopolymer melts equilibrated by hierarchical backmapping [12]. The systematic agreement between conformational and structural properties in blends prepared using backmapping and hybrid MD/MC, demonstrates that both methods can equilibrate the generic blends in our study. However, the identity-switch moves in the hybrid MD/MC scheme are impractical in cases of strong asymmetry between components [43]. In contrast, hierarchical backmapping can be applied to generic blends with strong asymmetries, such as polydisperse systems. Altogether, the preparation of reliable reference data for parameterizing blob-based descriptions and validating equilibration remains an important methodological issue for hierarchical backmapping. In certain cases the preparation of reference data could benefit from advanced MC methods, including rebridging algorithms [2,3].