Self-assembly in soft matter with multiple length scales

Spontaneous self-assembly in molecular systems is a fundamental route to both biological and engineered soft matter. Simple micellisation, emulsion formation, and polymer mixing principles are well understood. However, the principles behind emergence of structures with competing length scales in soft matter systems remain an open question. Examples include the droplet-inside-droplet assembly in many biomacromolecular systems undergoing liquid-liquid phase separation, analogous multiple emulsion formation in oil-surfactant-water formulations, and polymer core-shell particles with internal structure. We develop here a microscopic theoretical model based on effective interactions between the constituents of a soft matter system to explain self-organization both at single and multiple length scales. The model identifies how spatial ordering at multiple length scales emerges due to competing interactions between the system components, e.g. molecules of different sizes and different chemical properties. As an example of single and multiple-length-scale assembly, we map out a generic phase diagram for a solution with two solute species differing in their mutual and solvent interactions. By performing molecular simulations on a block-copolymer system, we further demonstrate how the phase diagram can be connected to a molecular system that has a transition from regular single-core polymer particles to multi-core aggregates that exhibit multiple structural length scales. The findings provide guidelines to understanding the length scales rising spontaneously in biological self-assembly, but also open new venues to the development and engineering of biomolecular and polymeric functional materials.

Simple molecular self-organization, such as aqueous micellization and emulsion formation, can be easily explained at the level of the interplay between water entropy, the relevant surface tensions, and the corresponding free energies [22]. Multiple structural length scales emerge in the presence of competing interactions [8,23], yet fundamental theoretical underpinning to deeply understand this is lacking. We present here a theoretical framework demonstrating how competition between molecular interactions can lead to spontaneous formation of structurally complex self-assembly matter exhibiting multiple-length-scale structural features and demonstrate its connection to a block-copolymer system. We employ energy minimization principles based on classical density functional theory (DFT) [24][25][26][27] to model the spatially varying average density of each component in the thermodynamic equilibrium in a multicomponent soft matter system. The interactions are approximated by a free-energy potential. This approach provides a generic phase diagram exhibiting a variety of possible equilibrium assembly configurations as a function of the system composition. We present a formulation that covers both the formation of single length-scale, FIG. 1. Generic assembly phase diagram of a 2D model binary system. Four distinct phases emerge: (i) a phase where both b and s particle densities are spatially uniform (squares), (ii) a phase corresponding to one species condensing into droplets surrounded by the other, i.e. a single-core phase (dots) (results for the large dot shown in Fig. 3), (iii) a multi-length-scale phase (stars) in which the species assemble into droplets inside droplets with the response corresponding to two subsets: (iiia) multi-core micelle type assembly (results for the large star shown in Fig. 2) where each s island is clearly separate from the other s islands and (iiib) where the s islands overlap at their outskirt regions, and (iv) a phase where the s particles form a hexagonal lattice covering the whole space while the b particles are homogeneously spread (rhombi). This corresponds to species b being soluble while s forms droplets. Purple stars are hybrid cases, where most of the islands of s particles behave as the green star states. However, some random clusters do not display the short spacing and the distribution of the s particles is similar to the one of the b particles (see Fig. 2(a)). The dotted line represents the linear instability of the length scale associated with the s particles, ks, and the dashed-dotted line the one associated with the b particles, k b . The boundary lines are guides to the eye.
For computational efficiency, we consider a twodimensional (2D) model system composed of two species in a solvent, however, the formalism is also valid in 3D. The system is characterized by soft interactions, where the effective pair potentials between coarse-grained complex molecules are designed using the generalized exponential model of index n (GEM-n). Such potentials have been extensively used to describe the effective interactions in a vast variety of polymeric systems [32][33][34][35][36][37][38][39]. GEM-n models exhibit interesting pattern formation if they are part of the Q±-interactions class [40], i.e. for n > 2. This theory applies generally to any molecular system that has three partially immiscible components. One of the components can be considered as the solvent, implicitly present in the interaction potentials of the two explicit species. Such systems are common in e.g. aqueous polymer mixtures, biomolecular and colloidal sys-tems, emulsions, and various liquid-crystal-forming systems.
The statistical distribution of the average densities ρ i (r) are obtained using the Ramakrishnan-Yussouff approximation [41] (details in the Supplementary Material (SM)). The fundamental condition for thermal stability of any self-assembled structure is that it corresponds to a minimum in the relevant thermodynamic potential. We choose a mixture of by big (b) and small (s) coarsegrained particles interacting via the effective potentials φ bb (r) = ε bb e −(r/R bb ) 4 ; φ ss (r) = ε ss e −(r/Rss) 8 ; where the subscripts refer to the interacting particle pair. ε bb > 0 and ε ss > 0 represent repulsion between pairs of b and s particles, respectively, corresponding to energy penalty for overlap due to entropic effects. ε bs imposes an attraction between b and s, whereas r is the distance between the centers of the particles. In water solutions, such effective attraction typically rises between hydrophobic components as a result of entropic contributions. It could also be generated by charge distributions on the molecules. Specifically, for charged species such as polyelectrolytes, this interaction can be easily tuned by salt. A long tail repulsion via ε + bs > 0 is useful in preventing singularities. R bb and R ss define the sizes of the particles (typically comparable to the radius of gyration for polymers). Similar theoretical approaches have been used to describe other forms of self-assembly, such as superlattice structures [42] and quasicrystals [43]. In order to choose a relevant set of parameters for Eq.
(1), we consider the linear dispersion relation branches ω ± (k), derived in Ref. [44], which characterise the growth or decay of density modulations, with wave number k, in the liquid state. The dependence of this quantity with the different parameters in Eq. (1) is rather complicated. Before investigating such a complex quantity, we start by assuming that the mixture is formed by two independent components. This can be ensured by carefully selecting the form of (1). In order to have two clearly distinct length scales we set the ratio k s /k b ≈ 20, where k s and k b correspond to the large and small particle characteristic wave lengths, respectively. We set R bb = 5 and R ss = 0.25, and use dimensionless units throughout. When considering two separate systems of b and s particles, respectively, the linear dispersion relation for each component reads [45][46][47][48][49] where D i is the diffusion coefficient and S i (k) the liquid structure factor of component i, which can be measured in experiments, and is accessible through S i (k) = 1/(1 − ρ iĉi (k)) [26], whereĉ i (k) is the Fourier transform of the direct pair correlation function (see SM for details) and ρ i is the bulk density of species i. With this FIG. 2. Dual-length-scale self-assembled phase (large star in Fig. 1, ρ b = 0.9, ρs = 0.55). For polymer systems, this refers to multi-core micelles, and for an emulsion to a thermodynamically stable droplet-inside-droplet state (multiple emulsion). (a) shows ρ b (r), (b) ρs(r), (c) the total density ρt(r) = ρ b (r) + ρs(r), and (d) a blow-up of a single peak. ρ b (r) exhibits hexagonal structure with spacing ≈ R bb . Cocentric with ρ b (r) peaks, also ρs(r) forms smaller hexagonal clusters with spacing ≈ Rss.
βε bb βε bs βε + bs βεss R bb R bs Rss 0.45 -0.45 0.02 70 5 3 0.25 approach, one can easily obtain the values of ε ii necessary to achieve linear instability for a given critical value of densities ρ c i , i = b, s, and vice versa. It is convenient (but not mandatory) to set ρ b ≈ ρ s . We choose βε bb = 0.45 and βε ss = 70, where β = 1/(k B T ), k B being the Boltzmann constant and T the temperature. This choice implies that under the hypothesis of two independent particle systems, the critical densities are ρ c b ≈ 0.55 and ρ c s ≈ 0.78. These values are obtained by finding the corresponding bulk densities equivalent to the onset of linear instability in Eq. (2), i.e. ω i (k) = 0 and dω i (k)/dk = 0. A common choice for the cross interaction radius is R bs = 1 2 (R bb + R ss ) ≈ 3. In order to avoid a singularity in the density distribution, one should restrict to 2π ∞ 0 drrφ bs (r) > 0. However, negative values of the integrated strength can be compensated by the strength of repulsion arising from φ bb (r) and φ ss (r). Furthermore, the attractive part of φ bs (r) should be strong enough to favor mixing, and the repulsive part must be small enough to avoid phase separation. Here we choose βε bs = −0.45 and βε + bs = 0.02. Designing Eq. (1) with the conditions provided above allows us to have a mixture in which ρ c i does not influence ρ c j , for i = j. This property is evident in the fact that the instability lines in the phase diagram of Fig. 1 (to be explained below) are almost vertical and horizontal, respectively. This condition is not required; however, it is very helpful in avoiding complications due to the interplay of φ bs in ω ± (k). In fact, under some conditions this interplay could, e.g., partially or fully suppress one or both instabilities. The model parameters used in Eq. (1) are summarized in Table 1. We emphasize that results similar to those obtained here are expected for different parameter sets, but also for models different from Eq. (1). The linear dispersion relation contains all the information required to determine whether or not the uniform liquid is stable with respect to any (small) perturbation. Thus, to resolve thermody- namically stable states where also assemblies with multiple length scales emerge, we address ω ± (k) and look for the crossing point of the two linear instability lines (the dashed and dashed-dotted almost horizontal and vertical lines in Fig. 1). We could have used also Eq. (2) to approximate the crossing point of the two linear instability lines, but ω ± (k) provides the exact information. Unfortunately, knowing this crossing point merely suggests the values of ρ b and ρ s that potentially exhibit equilibrium multiple-length-scale structure. The full assembly phase diagram, Fig. 1, based on Eq.
(1) and the set of parameters in the Table I, is obtained by varying ρ b and ρ s . An extensive description of this process is reported in the SM. To facilitate comparison with computer simulations on multicomponent polymer systems in a solvent, we adopt for description of the phases terminology that is specific to polymer selfassembly. The diagram consists of four distinct phases: (i) a phase where both densities are uniform (squares); (ii) a core-shell (or single-core) micelle phase (dots); (iii) a multi-length-scale phase (stars), composed by two subsets (iiia) and (iiib); and (iv) a phase in which the s particles form a hexagonal structure across the system while the b particles are homogeneously spread (rhombi). The boundary lines between the different states in Fig. 1 are meant as a guide for the eye, and should not be considered as exact.
We start the description of the different phases with subset (iiia), commensurate with multi-core-micelle formation in polymer assemblies. An example of this structure is shown in Fig. 2 (large star in Fig. 1). Here the b particles form a hexagonal structure with lattice spacing of ≈ R bb , as shown in Fig. 2(a). Simultaneously, the s particles form islands of small hexagonal clusters with lattice spacing of ≈ R ss , centered at the density maxima of b (panel (b)). Importantly, the orientations of the s islands are independent of each other at each b maximum. Thus, the s particles are locally ordered within each b maximum only. Such assembly response describes complex ordered phases in soft matter such as multi-core micelles or multiple-emulsions. Note that the perfect hexagonal order is due to periodic boundary conditions and lack of thermal fluctuations in the present calculations.
The structural change from (iiia) to (iiib) consists of the s particles islands increasing in size ( Fig. 2(b)). Increasing ρ s eventually bridges the islands and makes them indiscernible. As the s islands differ in orientation, grain boundaries will emerge. We mark the states with at least two merged islands with orange stars. The extreme case is a structure where the s particles fill completely the space (large values of ρ s ).
In contrast, decreasing ρ s starting from states in (iiia) causes the hexagonal arrangement of small islands s to melt into large s droplets centered around the maxima of ρ b (r). However, the attractive nature of φ bs forces the s particles to remain in the vicinity of such maxima ( Fig. 3(b) and inset of Fig. 3(a) show an example corresponding to the large dot in Fig. 1). This transition corresponds to moving from phase (iii) to phase (ii). The latter structure is commensurate with single-core (coreshell) micelles. These states survive for decreasing ρ s unless ρ b < ρ c b . Figure 1 shows that multi-core micelles are obtained mostly for ρ s < ρ c s , i.e. for densities at which the system is linearly stable with respect to k s . The deviation of the actual phase boundaries from the linear instability line can be explained as follows; consider a single-core-micelle state: if the average density of the s particles in each of these peaks is roughly of the same order as ρ c s , the s particles find energetically favourable to form clusters with typical spacing of ≈ 2π/k s ((ii) → (iii)). This can be seen as a local instability occurring at the level of a single b site. The two purple stars in (iiia) represent a state in which most of the s clusters behave as in the case represented by green stars (Fig. 2(b)). However, some random clusters of s particles do not display the short length scale order but instead their distribution is similar to the one of the b particles (Gaussian-like distributed) (Fig. 2(a)). These states might be metastable due to the vicinity with phase boundaries. In (iv), the system is filled by a hexagonal structure with spacing of ≈ R ss . Simultaneously, ρ b (r) is uniform. This phase can be found for ρ s > ρ c s and ρ b < ρ c b . The attractive nature of φ bs favours the orange-star states to exist down to values of ρ b < ρ c b , where one would expect to find rhombi states if only linear instabilities were considered. Cluster formation of b particles is enhanced by the attraction mediated by clusters of s particles. This can be explained in a manner similar to the discrepancy between the linear instability line of the s particles and the boundary between dots and stars. Finally, in phase (i) both species are uniformly distributed, i.e. a solution of miscible species. The slight discrepancy between the linear instability line of particles b and the boundary between the homogeneous state and the states in which the latter species is linearly unstable is due to finite-size effects and to the fact that the dispersion relation only considers the linear contributions of the dynamics. Furthermore, the size of the system has been chosen randomly which may lead to a wave number slightly different from k b . ω + (k b ) is concave with a maximum at k b at the onset of the instability: particles b will form clusters for a slightly different value of ρ b . The symmetries of ρ b (r) and ρ s (r) in the different phases are summarized in Table II. We note that the transition between phases (ii) and (iii) is discontinuous. This is due to spatial symmetry breaking between these phases. To compare the emergence of the above-discussed single vs. multiple-length-scale assembly in a realistic molecular level system, we perform dissipative particle dynamics (DPD) [50][51][52][53] simulations of a polymeric solution composed of a solvophobic polymer (A 19 ) and a di-block copolymer (A 1 B 6 ) in a solvent using the LAMMPS [54] simulation package. The segments A and B are solvophobic and solvophilic, respectively. The subscripts refer to the number of DPD beads (block lengths). Each bead represents a group of atoms which experience a force where F C describes conservative interactions, F D dissipative contributions, F R random contribution and where N is the total number of beads. The interaction forces are treated as pairwise additive and are truncated at a distance r c . In the polymer chains, adjacent beads also contribute to a spring-like force F S . These soft potentials facilitate acceleration of the numerical simulations so that realistic experimental time and length scales can be achieved. Further details and simulation specifics can be found in the SM. Figure 4 shows snapshots of different equilibrium configurations obtained from simulations by varying the molar fractions of the components. This system strikingly self- The subscript refers to the number of DPD beads. Left: single-core (core-shell) configuration for molar fractions 50% A19 and 50% A1B6. Right: multi-core configuration for molar fractions 10% A19 and 90% A1B6. The solvent is not shown for clarity. assembles into structures corresponding to a multi-corecorona and a single-core-corona states. Both cases correspond to a total solid content of N s /(N s + N w ) = 20.15% in aqueous solvent, where N s and N w are the number of solid and water beads, respectively. In the former case the molar fractions are 10 % A 19 and 90 % A 1 B 6 , whilst in the latter 50 % A 19 and 50 % A 1 B 6 . The simulations are performed in a cubic box of 100×100×100 r 3 c . These equilibrium configurations show dual and single-lengthscale structural assemblies analogous to those in our DFT based phase diagram (phases (ii) and (iiia) in Fig. 1). It is tempting to interpret a DPD polymer chain as a single coarse-grained DFT particle, but it should be noted that the surfactant-like nature of the copolymer makes a quantitative comparison between chemically predictive molecular model and DFT-assembly structures difficult. Nevertheless, the two distinct length scales emerging from the effective interactions in Fig. 2 also appear in the DPD model. Additionally, also in DPD simulations the key feature for multi-core assembly is a sufficient degree of immiscibility between the species. The formation of solvophobic cores needs to be energetically favorable and need to be sufficiently stabilized by the solvophilic polymer segments, and the solvophilic polymer cannot have too favorable interactions with either the solvent or the solvophobic polymer. If either the solvophobic polymers are too solvophobic, or the solvophilic one too solvophilic, or the two polymers readily mix, the assembly becomes core-shell or phase separates. This demonstrates how competition between interactions and mutual balance of immiscibility gives rise to multiple-length-scale assemblies -single-length-scale self-assembly is retained when any of the pairwise effective interactions dominates.
To summarize, we have presented a microscopic theory capable of explaining self-assembly within soft matter with multiple competing length scales. The theory relies on the interplay between the effective interaction potentials modelling the different constituents of the system. The significance of our findings resides in the identification of key molecular features involved in general multiple-length-scale self-assembly. Although multiple-length-scale assembly is ubiquitous, such guidelines have been lacking until now. By providing much needed insight to understanding, e.g., morphologies rising in intracellular biocondensates, multiple emulsions, or lipid bilayer microdomains, we obtain means to engineer and tune soft matter to desired multiple -length-scale structures. In addition to obvious applications to multiple-length-scale self-assembling compartmentalization, such as drug delivery, catalysis and selective multistep reaction platforms, our work provides tools to harnessing the full potential of revolutionary materials production via biological mechanisms, advanced engineered biomaterials, or complex polymer assemblies, providing crucial insight on how to tune the assembly response.

Density functional theory for mixtures and the Ramakrishnan and Yussouff approximation
The grand potential of a system composed of two types of particles is where F is the intrinsic Helmholtz free energy functional, ρ 1 = ρ 1 (r) and ρ 2 = ρ 2 (r) are the spatially varying densities of the two species, V ext i (r) is the one-body external potential acting on species i (for bulk systems V ext i (r) ≡ 0 for i = 1, 2) and µ i are the chemical potentials. The intrinsic Helmholtz free energy can be split into two terms where the first term is the ideal gas contribution, i.e.
and Λ i is the thermal de Broglie wavelength, d is the dimensionality of the system, T the temperature and k B the Boltzmann constant. The second term in Eq. (S2) is the excess Helmholtz free energy arising from the interactions between the particles. Following Ramakrishnan and Yussouff [41], we approximate this functional by a functional Taylor expansion around the homogeneous fluid states ρ 0,i . A truncation of the series expansion at second order gives where δρ i (r) = ρ i (r) − ρ 0,i and µ exc are the excess chemical potentials. The pair direct correlation functions c ij (r) are obtained via the random phase approximation (RPA). The RPA consists of assuming a simple mean-field form for the excess free energy in Eq. (S2). This leads to c ij (r) = −βφ ij (r), where φ ij (r) are the effective pair potentials, and β = (k B T ) −1 [33]. For additional accuracy, one could, for example, use the hypernetted chain (HNC) Ornstein-Zernike integral equation theory [24]. The equilibrium density profiles ρ i (r) are those which minimise the functional of the grand potential Ω[ρ 1 , ρ 2 ] and which therefore satisfy the following pair of coupled Euler-Lagrange equations for i = 1, 2.
Difference in the grand canonical potential energy Substituting Eqs. (S2-S4) into Eq. (S1) one obtains an expression for the grand canonical potential. For homogeneous system, i.e. ρ 1 (r) = ρ 0,1 and ρ 2 (r) = ρ 0,2 , the latter reduces to The difference in the grand canonical potential energy between the equilibrium state and the corresponding homogeneous state becomes dr dr δρ i (r)c ij (| r − r |)δρ j (r ).

(S7)
Dissipative particle dynamics simulations (DPD) The DPD simulation method, originally proposed by Hoogerbrugge and Koelman [52], is a mesoscale coarse-grained bead-based molecular simulation technique. Combining aspects of molecular dynamics and lattice-gas automata, DPD acknowledges the idea that different beads can overlap, modelling non-bonded interactions with soft repulsive potentials [52,53]. The soft potentials employed to describe the interactions between the beads allow the simulations to reach realistic experimental time and length scales for e.g. block-copolymer self-assembly systems. Each DPD bead represents a coarse-grained region in the molecular system (e.g. several monomers of a polymer, or a solvent region, or a group of atoms) which experience a force where F C describes conservative interactions, F D dissipative contributions, F R a random contribution, and N corresponds to the number of DPD beads in the system. The interaction forces are treated as pairwise additive and are truncated at a distance r c . The conservative force is a soft repulsive force acting along the centers of two DPD particles, given by where a ij is the maximum repulsion between beads i and j, r ij =| r i − r j | /r c is the cut-off normalized distance between beads i and j, andr ij = r ij /r ij gives the force direction via a unit vector. The coefficient a ij is connected to the Flory-Huggins mixing parameter χ ij via the relation χ ij = (a ij − a ii )/3.27 at a density ρ = 3. The choice a ii = 25 for the repulsion parameter between beads of the same species (i.e., χ ii = 0) is common and based on the compressibility of the dilute solution [52]. An a ij value exceeding the self-repulsion in magnitude corresponds to a stronger bead-bead repulsion between beads of different species. The dissipative force is given by where γ is a viscosity related parameter (γ = 4.5), ω D is a weight function that goes to zero at r c , and the relative velocity is v ij = v i − v j . The random force is given by where ξ ij is a zero-mean Gaussian random variable of unit variance and σ 2 = 2γk B T . The weight functions follow the relation ω D (r ij ) = ω R (r ij ) 2 = (1 − r ij ) 2 for r ij < r c . Consecutive beads in a polymer chain also perceive a spring force F S i defined by where κ is the spring constant, which is set to κ = 80 in this work, the equilibrium distance r 0 = r c , and where j * refers to the nearest neighbours in the chain.
The system modelled in this work is composed of a mixture of a solvophobic polymer (A 19 ) and a di-block copolymer (A 1 B 6 ) in a solvent (S) medium. The nomenclature A, B, S refers to the DPD beads such that A is a solvophobic bead, B a solvophilic bead and S the solvent bead. The subscripts refer to the number of beads in the chain. For the interactions between beads of same type, we use a AA = a BB = a SS = 25. This value is based on the compressibility of the dilute solution [52]. The interactions between the beads in the simulations are a AB = 72, a AS = 115 and a BS = 30. To simplify the DPD equations and simulations, the cutoff radius r c , the bead mass m, and the energy k B T are reduced to r c = m = k B T = 1 which leads to the time unit τ = (mr 2 c /k B T ) 1/2 = 1. The DPD simulations were performed using the LAMMPS [54] package. A modified version of the velocity-Verlet algorithm is used to integrate the equations of motion. A time step of ∆t = 0.05τ is used. The simulations are performed in a cubic box of 100 × 100 × 100 r 3 c in size. Periodic boundary conditions in all directions in 3D were used. The system is initialized with setting the polymers and the solvent beads with random placement in the box (overlap excluded in initialization). The simulations were carried out for 2 × 10 6 time steps, and equilibration checked for via analysis of the time evolution of the assembly structures.