Hierarchical self-assembly of telechelic star polymers: from soft patchy particles to gels and diamond crystals

The design of self-assembling materials in the nanometer scale focuses on the fabrication of a class of organic and inorganic subcomponents that can be reliably produced on a large scale and tailored according to their vast applications for, e.g. electronics, therapeutic vectors and diagnostic imaging agent carriers, or photonics. In a recent publication (Capone et al 2012 Phys. Rev. Lett. 109 238301), diblock copolymer stars have been shown to be a novel system, which is able to hierarchically self-assemble first into soft patchy particles and thereafter into more complex structures, such as the diamond and cubic crystal. The self-aggregating single star patchy behavior is preserved from extremely low up to high densities. Its main control parameters are related to the architecture of the building blocks, which are the number of arms (functionality) and the fraction of attractive end-monomers. By employing a variety of computational and theoretical tools, ranging from the microscopic to the mesoscopic, coarse-grained level in a systematic fashion, we investigate the crossover between the formation of microstructure versus macroscopic phase separation, as well as the formation of gels and networks in these systems. We finally show that telechelic star polymers can be used as building blocks for the fabrication of open crystal structures, such as the diamond or the simple-cubic lattice, taking advantage of the strong correlation between single-particle patchiness and lattice coordination at finite densities.


Introduction
Technological advances strongly rely on the development of novel materials, characterized by specific physical-chemical properties. Recent discoveries in the field of tunable nano materials have paved the way for the creation of compounds with unique properties, that can be controlled and manipulated at different length scales to satisfy specific application-determined requirements. In order to be able to do so, a deep understanding of, and full control over, how structure and composition dictate properties and performance is indispensable.
In past years, much effort has been dedicated to creating tunable building blocks that can self-assemble and stabilize complex structures with particular acoustic or photonic properties, for example the diamond lattice [1]. This has motivated materials scientists to launch a largescale analysis of building blocks of various shapes [2] and functionalizations [3][4][5][6][7]. A wellknown illustration is found in colloids functionalized with attractive or repulsive regions (patches) [8,9] which show an interesting self-assembly behavior that can be tuned either by changing the shape of the patches or their relative orientation [6]. State of the art methods for the synthesis of such particles include lithography [10,11], microfluidics [12] or glancing angle deposition [13,14].
We propose and extensively study a novel class of highly versatile tunable building blocks, the so-called telechelic star polymers (TSPs), as an ideal type of candidate molecules to be used for the self-assembly of novel and designable soft-materials. The TSPs are star polymers whose arms are formed by diblock copolymers. Every arm is tethered to a central anchoring point, whereby the internal monomers (heads) are solvophilic and the end monomers of each arm (tails) are solvophobic. To grasp the general properties of TSPs and not to be impeded by finite size effects, we will focus here on star polymers whose arms are in the so-called scaling regime. Hereto, each arm necessarily consists of a large number of monomers M 1. Consequently, the single star properties, such as the intra star aggregation behavior, only depend on the percentage α of attractive monomers, the temperature of the system, the quality of the solvent and the number of arms f . The rest of this paper is organized as follows. In section 2 we present the details of the model and explain the procedure which led from its microscopic basis to its coarsegrained rendering. In section 3 we discuss our findings concerning the self-aggregation properties of the system, focusing on the self-organization of the single molecules into patchy particles. Thereafter, in section 4 we describe the quantitative characteristics of intermolecular aggregation into a percolating gel, whereas in section 5 the morphology and connectivity of the molecules in the ordered diamond phase are analyzed. Finally, in section 6, we summarize and draw our conclusions, also briefly discussing the outlook for future work.

The model and the coarse-graining strategy
The system we are considering, as well as the questions we aim to address, require a detailed knowledge of the phenomena that take place at both the microscopic and mesoscopic levels, and therefore the best strategy is formed by computer simulations, employing various special techniques. On the atomistic level, the systems we address are composed of a huge number of monomeric units; it is therefore essential to use a coarse-grained representation to describe the system. Such a methodology needs to be precise, to be back-tractable onto the microscopic system and at the same time it has to provide access to large scale analysis both in terms of the number of molecules in solution as well as in terms of the number of microscopic constituents.
The coarse-graining methodology used to address the problem at hand is a first-principles approach, which generalizes the considerations previously introduced successfully for the case of linear diblock copolymers [15,16]. In the same spirit, we adopt a multi-blob representation to describe each diblock copolymer arm of the TSP.
In the microscopic representation, each AB-diblock arm, of total degree of polymerization M, contains M A monomers of kind A and M B monomers of kind B. A microscopic implicit solvent model is defined by specifying the solvent-averaged interactions between A−A, A−B and B−B pairs of monomers, as well as the bonding potentials between adjacent ones. Solvent selectivity is reflected in the details of these various interactions. The microscopic description commences with a monomeric-resolved model, in which the A block behaves as a random walk on a simple-cubic, microscopic lattice of spacing a, corresponding to good solvent conditions. The B block, on the other hand, is a random walk on the same lattice with nearest neighbors interactions (attractions), corresponding to -like or poor solvent conditions. In what follows, the results to be presented pertain to a temperature that is slightly below the -point of the nearest-neighbor square-well-model that describes the terminal blocks.
The key parameters to describe each of the f -arms in a TSP are the total number of monomers per arm, M = M A + M B and the asymmetry ratio: Within the multiblob representation, each of the two blocks, A and B, is mapped onto n A and n B blobs respectively, where every blob contains a number m A = M A /n A or m B = M B /n B of monomers of the underlying microscopic model. Each blob will have a radius of gyration r gγ = ab γ m ν γ γ with the subscript γ = A, B denoting the monomer species, ν γ being the critical exponent of the γ species, a the monomer size and b γ being model-specific numerical coefficients of order unity. Rewriting equation (1) in terms of blobs and their radii of gyration, we obtain a relation between the asymmetry ratio, the number of A and B blobs, and the radii of gyration r gA , r gB . In fact, without loss of generality, we can require the A-and B-blobs to have a common gyration radius, r gA = r gB = r g to obtain Equation (2) is then complemented with the conservation law for the total number of monomers M, i.e.
so that equations (2) and (3) can be solved for any given combination of α and M and for fixed values of r g , b A and b B to determine the numbers of blobs n A and n B . Thus, the information on the asymmetry only depends on the number of blobs n A and n B that is needed to represent the A and the B blocks. Blobs interact with each other via effective potentials, which are computed by means of interacting blob-dimers at infinite dilution, taking into account interactions that involve up to four bodies [15] at the blob level. Such a methodology has been extensively tested in the last few years, and it was proven to reproduce microscopical predictions, both when tested on full monomer simulations in complicated geometrical constraints (e.g. homopolymer brushes [17]) or when compared to experimental results for diblock copolymers in the semidilute regime [15]. In the case of linear diblock copolymers, the minimum number of blobs n = n A + n B to be used within the semi-dilute regime was fixed by the ratio between the density of polymers in solution and the overlap density of the blobs, ρ * b ≡ 3/(4πr 3 g ). In the case of TSPs, there exists a wide range of densities within the semi-dilute regime for which the monomer concentration in solution is lower than the intramolecular monomer concentration, due to the fact that the molecular architecture forces a high local monomer concentration close to the star cores. Under such conditions, the minimum number of blobs, n min , to be used to represent a star is dictated by the intramolecular monomer density. In order to assess the value of n min , several simulations were performed and key intramolecular properties, i.e. the radius of gyration of the star, of the arms, of the attractive and repulsive parts of the arms, were tested upon increasing the number of blobs. The quantity n min was determined as the minimum number of blobs beyond which the aforementioned structural quantities became independent of the number of blobs employed in the coarse-graining procedure [18].
In order to arrive at the coarse-grained model description, the various blob-blob interactions have to be introduced. At this level of detail, interactions are always between pairs of blobs, i.e. A-A, A-B and B-B. Hence for each pair α and β, there is a direct interaction by means of the pair potential V αβ (r ). This does not yet take into account that neighboring blobs in a polymer chain are connected as well, therefore the pair potential must be supplemented by three tethering potentials between bonded blobs, namely ϕ AA (r ), ϕ AB (r ), ϕ BB (r ). Here, the quantity r denotes the separation between the centers of mass (CM) of the corresponding blobs, which are employed as the dynamical degrees of freedom in the coarse-grained simulation, replacing the underlying coordinates of the microscopic monomers.
These effective interactions between blobs are determined by a first principles inversion procedure, that was derived for diblock copolymers [15,16], generalizing the method used earlier for the simple dumbbell representation of the same [19][20][21]. The three intramolecular ϕ AA (r/r g ) ϕ AB (r/r g ) ϕ BB (r/r g ) Figure 1. The non-bonded potentials V AA in black, V AB in red and V BB in blue and the tethering potentials ϕ AA in black, ϕ AB in red and ϕ BB in blue between the CM of the various blobs of type A and B. The radius of gyration of the A and B blobs is the same and labeled as r g .
tethering potentials, consisting of a superposition of repulsive and attractive terms, are determined from microscopic, full-monomer Monte Carlo (MC) simulations of an isolated α-β diblock copolymer, α, β = A, B. The distribution function s αβ (r ) of the separations between the CM of the two blocks is estimated by averaging over a large number of monomer-configurations, and the corresponding tethering potential follows from ϕ αβ (r ) = −k B T ln s αβ (r ) .
In order to determine the intermolecular pair potentials, we consider the six possible combinations of αβ and γ δ dimers and calculate the corresponding blob-blob pair correlation functions h αγ (r ), as functions of the distance r between the CMs of the α-block of dimer 1 and the γ -block of dimer 2. This is achieved from MC-generated histograms, by averaging over allowed monomer configurations for fixed values of r , according to the usual Metropolis algorithm. The functions h αγ (r ) are mapped out by varying the distance r , i.e. by gradually moving the CMs toward each other. The effective pair potentials are extracted by inverting the pair correlation functions according to an exact cluster expansion, valid for an isolated pair of dimers, i.e. in the low density limit [15,16,[19][20][21][22]. The effective intermolecular potentials V αβ (r ) and the intramolecular tethering potentials ϕ αβ (r ) obtained by the inversion procedure are shown in figure 1, as functions of the CM-CM distance r , reduced by the common blob radius of gyration r g . The molecular weights of the polymer segments comprising each blob, used in the simulations, are sufficiently large (m A , m B 100) to justify the statement that the data shown in figure 1 indeed correspond to the scaling limit.
We performed MC simulations in which each elementary move consists of a single blob displacement. The maximum move extension was chosen to be r ≈ 1.5r g , where r g is the radius of gyration of each blob. This particular choice has been made in order to mimic as closely as possible a molecular dynamics simulation featuring the soft blobs as fundamental entities.

Dilute regime: building the elementary blocks
The prominent characteristic of TSPs rests on the possibility to control the inter-and intramolecular association processes at different length scales. With a bottom-top strategy it is possible to start from the synthesis (computationally and experimentally) of the single telechelic stars, increase the concentration in solution and obtain specific solid phases or gel/network structures. We define as 'gels' all those systems in which stars assemble in amorphous percolating structures and where the motion of the anchoring point of the stars is extremely slow. The system, however, is not fully arrested during the simulation run. This is evidenced by the fact that, for instance, the voids of the gel do change their size, and that we see holes in the gel open up and their sizes fluctuate. In this sense, the gel is an active, physical gel with bonds that rearrange within it. Nevertheless, our terminology here is not meant in the dynamical sense, since we cannot access the true dynamics, but rather in the sense of an amorphous, multiply connected patchy polymer network that spans the whole simulation box. The analysis performed in this paper is done by means of MC simulations; nevertheless the small MC moves used to equilibrate the system allow to see pseudo dynamical properties of the stars. We found that there is a clear time separation between the intra-and inter-star 'dynamics'. In particular, for the 'gel systems', individual star arms frequently disconnect from one cluster and join another while the 'dynamical' motion of the anchoring point appears to be extremely slow.
In [18], we have shown how to build the elementary blocks: the self-aggregating patchy behavior can be fully controlled by selectively tuning either or both the number of arms per star f and the fraction of attractive monomeric units α at the free ends of the arms. By choosing, in a selective solvent, the right combination of functionality f and percentage of attractive monomers α, a single TSP self-assembles into a soft patchy particle that preserves its character (number of patches, average angle between the patches, average extension of the patches from the anchor point) upon increasing density, from the extremely dilute case up to the semi-dilute regime. Figure 2 sketches the single-molecule conformational state diagram. As already put forward in [18], in the case in which the thermal, solvophobic tails lie slightly below their -point, one finds, for a fraction of attractive monomers per arm α smaller than 0.3, that the macromolecules remain in an open star configuration even when the number of arms is augmented. In contrast, for α 0.3 the molecules self-assemble into soft patchy particles. Moreover, for a fixed functionality f , the number of patches p can be tuned by augmenting the percentage of attractive monomers α. Alternatively if α is kept constant, the number of patches can be selected by changing the functionality of the star. The aggregation into these soft patchy structures is fully determined by the balance between the entropic contributions arising from the intra-star repulsions in the core of the stars, and the competing, enthalpic terms arising from the solvophilic nature of the terminal monomeric-or blob-units.
In the second stage of the bottom-up process, we increase the concentration of the stars. As we will demonstrate in the next sections, the self-assembled TSP soft-patchy particles, while aggregating, maintain their individual properties such as size, number and arrangement of From an experimental point of view, this behavior has an important consequence for applications, as it is not necessary to use laborious preparation techniques to bring about the desired patchiness, because they simply are the results of the self-organization of the particles at the molecular level. TSPs therefore appear as a novel class of extremely promising selfaggregating adjustable soft functionalized particles, experimentally accessible on a large scale, and, as we will discuss below in detail, are able to stabilize exotic crystal structures using a novel mechanism.

Gel phases: topological and structural analysis
The behavior of TSPs at low densities shows a well-defined self-patchiness as a function of the functionality f and the asymmetry ratio α. In fact, the TSPs have a strong preference for a discrete number of patches ordered around the center of the star for quite a wide range of the star parameters, suggesting that each self-organization of the star might be robust against Table 1. A summary of the intra-and interparticle aggregation properties for star polymers with f ∈ [3,15] and fixed asymmetry α = 0.4, over a wide range of densities. The infinite dilution case, ρ/ρ * = 0, pertains to intramolecular aggregation and the number of patches formed is given by the entries in the columns labeled with p. The same quantity at finite concentrations, ρ/ρ * > 0 as indicated, expresses the number of patches per particle in the gel phase, where strong intermolecular aggregation takes place. The first column denotes the correspondence between the values of the densities ρ 1 , ρ 2 , ρ 3 and ρ 4 for which various structural characteristics are plotted in figures 5, 7, 9 and 11, and the ratios ρ/ρ * for the quoted functionalities f .  Table 2. Same as table 1 but for the asymmetry ratio α = 0.6. In the first column, the values ρ 1 , ρ 2 , ρ 3 and ρ 4 are now those for which structural quantities in the gel are shown in figures 6, 8, 10 and 12. external changes, such as an increase in the density of the stars in solution. In other words, if a TSP would suddenly feel an additional arm coming from a second star, there is a good chance that this would not disrupt the local organization of the patches. In order to test for such robustness we performed a screening of TSPs solutions at various densities and for a wide range of ( f, α)-combinations. Hereto, simulations have been performed on systems with four different star-number densities labeled by ρ 1 , ρ 2 , ρ 3 and ρ 4 , that correspond to star-polymer numbers N = 50, 75, 100 and 150, respectively, enclosed in a cubic simulation box of volume V = × × , with box size = 1000a. Since the radius of gyration R g of the stars depends explicitly on the functionality and the asymmetry ratio, the star overlap density ρ * ≡ 3/(4π R 3 g ) will be different for the various systems. Accordingly, the corresponding densities in terms of these units are summarized in tables 1 and 2 and they are also denoted as ρ 1 , ρ 2 , ρ 3 and ρ 4 for future reference. The first feature that is observed on increasing the density is the capability of the TSP particles to form aggregates that, as the concentration is augmented further to moderately high values, cluster into percolating, gel-like networks, see e.g. figure 3. It is interesting to notice that systems with f = 3 present an average patchiness p 1. As a patch is defined as at least two arms merging together, stars with functionality f = 3 arms can form at most one patch assembled by two or three of their arms. As there also exist configurations for which each of the three arms is separated from the other two, and these correspond to p = 0, the expectation value of the patchiness is less than unity. Still, these stars can form a gel by merging the solvophobic ends of the isolated arms, therefore forming bridges between their anchoring points.
There are various ways to analyze the structure of ordered and disordered systems. One of the more obvious approaches is that of considering the radial distribution function in real space, which is mostly employed in computer simulations or in experimental observations with, for instance, confocal microscopy. Its counterpart in Fourier space, the structure factor, can also be considered and can be obtained from scattering experiments. The radial distribution function by default is, however, a two particle correlation function only, and hence every multi-body correlation is absent from it. As a consequence, it sometimes lacks the possibility to distinguish between similar structures and is not always sensitive to the precise form of structure.
In a recent study on polar dumbbells, it was found that such particles have a transition from a fluid to a low density, string-like network on lowering the temperature with some characteristics of a gel [23,24]. In order to analyze the open network structure, for which the radial distribution function is dominated by the presence of linear particle chains, the authors proposed to make use of the Euler characteristic as a topological fingerprint. The Euler characteristic χ is one of the four scalar Minkowski functionals that characterize a given surface embedded in three dimensions, the others being the enclosed volume, total surface area and the integral mean curvature. The Euler characteristic itself is proportional to the integral Gaussian curvature and its value is not affected by any continuous, topology-preserving deformations of the surface. What is special about this quantity is that, contrary to the other invariants, it can only assume an integer value χ (I) for any instantaneous configuration I of the system. This is related to an alternative interpretation of the Euler characteristic by means of a combination of number of disjoint surfaces, handles and cavities of the total surface that is considered.
A suitable surface that could be considered for computing the Euler characteristic is the interface between the clusters and/or dense gel with the empty spaces in the simulation box. However, this would hardly be useful in describing the morphology of the physical-association sites (i.e. of the patches) of the molecules in an adequate fashion. Instead, we introduce a family of surfaces parameterized by an appropriate length scale R that can be constructed from a given configuration of particles. To this end, we consider the surfaces A(R) formed by spheres with a radius R located at the CM of the attractive, B-blobs in the system, where only those parts are considered that do not lie within any other sphere. Hence, the surface is the full set of points whose shortest distance to any center of mass is given by R and uniquely assigns each point in space to a single surface only. For small values of R this will be just a collection of disjoint spheres, but on increasing the value of R, some of the spheres will touch, merge and later form rings and cavities. Eventually, the collection of spheres will fully occupy space and the surface will vanish. By computing the Euler characteristic for every surface A(R), we obtain the Euler characteristic χ (R) normalized by a factor 2N , with N the number of spheres that is considered. It should be noted that the choice made for the family of surfaces is not unique; in fact different choices of families of surfaces will result in different values of the Euler characteristic. It is, however, a choice that reflects the spherical nature of the particles in the simulation as well as one that allows for an efficient computational implementation. The complexity of the analyzed surfaces based on the particle positions automatically includes a sensitivity of this measurement with respect to many-particle correlations and provides insights in both the local and global structure. For a more detailed description of the measurement and its implementation we refer the reader to [24,25].
Examples of appropriately normalized χ(R) at finite densities are given in figure 4 for systems of TSPs with f = 3, 5, 7 and 10 arms and α = 0.4, 0.6 and 0.8. For these particular graphs, only the attractive blobs have been considered in the calculations. Each curve is obtained from a single configuration, and since the computation is exact, no errors are indicated. The noise and minor fluctuations visible are caused by inhomogeneities within such a single configuration and would disappear on averaging over several ones. With the normalization that is imposed, we find that by construction χ(R = 0) = 1, because in the limit of small R, the surface A(R) consists of N disjoint spherical contributions. On increasing the value R, the Euler characteristic of a particular configuration is only affected if the topology of the surface changes. This occurs if and only if one of three possible events takes place. The first of them is the external touching of two spheres, in which case either the number of disjoint surfaces is diminished or a chain of joined spheres forms a loop. In the second type of event, three spherical surfaces, that were joined to form a loop or torus-like shape, become so large that the hole or handle in the surface between them disappears. The last type of event occurs when on increasing the value of R a cavity formed by at least four spheres vanishes. The occurrence of the first or third type of event decreases the Euler characteristic by 2, whereas the second type of event will increase it by 2, with our normalization this will be ±1/2N .
Variations in the slope of the Euler characteristic can therefore be interpreted in terms of these different types of events and their relative importance with respect to each other. In the  range R/r g 10, the Euler characteristic drops from unity to a small, possibly negative value. On these length scales the attractive blobs find one or more neighboring particles. It is, however, the absence of variations in the Euler characteristic that is the most informative in this particular system. This behavior can be observed in the plateaus found in range 10 R/r g 50, and is most pronounced for the f = 3 functional TSPs. This behavior is indicative of the formation of well-separated clusters of attractive blobs, and is caused by the self-associating, patch structures formed by individual TSPs that now possibly have merged with other star patches as well. The total number of thus formed clusters in the system can be expressed in terms of the height of the plateau as N χ (R), and therefore the average number of blobs per cluster as 1/χ(R). The location of the beginning of the plateau is related to the density of particles within a cluster, whereas the location of the end of a plateau is a measure for the nearest neighbor distance between different clusters. In order to make these relations more quantitative, one would need to make some assumptions on the shape and regularity of clusters, which can not be inferred from the Euler characteristic itself. Beyond the plateaus, the Euler characteristic can be described as a system where the clusters act as new, individual particles. Pair distribution functions g(r ) (not shown here) computed for gel-forming systems assembled by stars with different functionalities f and patchiness p did not show any significant features. Apart from trivial differences linked to the star-sizes, they appear to be insensitive to patchiness and therefore fail to convey any crucial information. It should be stressed that the Euler characteristic, as presented above, cannot reveal anything about the properties of individual TSPs at higher concentrations. Therefore, in order to confirm the preservation of the single star properties a more precise and detailed cluster analysis on the patches, their size, composition and relative position is required. In order to do so, we grouped the solvophobic particles of each arm into clusters defined by introducing a cut-off radius R cut : particles are classified as neighbors if their distance is below R cut and clusters are then built by grouping particles that share at least one neighbor. R cut is a free parameter, that is chosen in a range of values around R cut r g , where r g is the radius of gyration of the blob. The clusters are well separated in solution, consequently the optimal R cut can be chosen as the value such that the number and size of the clusters is not affected by a small variation R cut of the cut-off radius. The clusters are identified for each system configuration generated by the finite density simulations. Additional, detailed information on the distribution of key structural quantities in the gel phase for different combinations of functionalities and asymmetry, as well as for various concentrations ρ, is presented in figures 5-12. We measure the number of patches that belong to each star, their geometrical arrangement around the center of the star, and the degree of connectivity present in the system. These quantities are then averaged over all configurations. Considering the usefulness of TSPs as building blocks that retain their single-particle properties on augmenting the density, it is desirable that also the single-molecule properties obtained in the gel are similar to those obtained in the crystals. In fact, a comparison of these quantities for both phases reveals that the average properties are identical, and the main difference that can be observed is that in the gel the distributions are slightly more smeared around the average values than in the crystal phase. The minor difference in this behavior has its origin in both the softness of the particles and the inhomogeneous distribution of particles in space that is expected of a disordered gel phase. The cluster analysis is performed both at zero density and at finite density to fully characterize the properties both of the single star and of the solution.
As we have demonstrated in section 3, stars with given f , α parameters collapse either onto single or multiple patch configurations provided α 0.3. The actual number of patches formed by each particle is fully determined by the choice of the f , α combination. In order to demonstrate that the average single molecule properties are preserved upon augmenting density, the starting point of our density dependent analysis of the gels is the calculation of the average patchiness of the TSPs in solution.
The determination of the number of patches formed by any star was performed as follows. Firstly, a star-patch requires the presence within a cluster (defined in what follows) of at least two arms of a given star; isolated arms can participate in a patch formed by other stars but do not constitute a patch of the given star by definition. Secondly, by making use of our cluster counting algorithm, we analyzed the stored configurations obtained during the simulations at various densities, and we identified the groups of solvophobic blobs lying with respect to one another closer than a certain cut-off. The cut-off value was empirically identified as the distance above which the measured number of clusters did not change. This feature is supported by the Euler characteristic analysis and by visual inspection, both ascertaining that clusters are well-separated in space. Finally, the average number of patches is then defined as the number of clusters in which more than one arm per star participates, averaged over the sampled configurations. Evidently, the average number of patches formed by any star is not the same as the total number of clusters divided by the number of stars, since more than one stars can participate to a given cluster. In figures 5 and 6 we plot the distribution P( p) of the number of patches p that each individual TSP in solution displays upon self-assembling for four different star-number densities. In both figures, p is the number of patches per star averaged over configurations. Systems of stars for five different functionalities f and two different asymmetries α are analyzed, namely f = 3, 5, 7, 10 and 15 and α = 0.4 and 0.6. In this paper we will mainly present results for the two aforementioned representative asymmetry ratios, but all other systems with α ∈ [0.3, 0.8] analyzed do exhibit similar behavior. The data shown in figures 5 and 6 convincingly demonstrate that stars maintain their self-assembled average number of patches for the whole range of finite densities, since the results from different values of the latter essentially collapse on top of one another.
To characterize the average patchiness, each star is labeled and its properties are averaged over the complete simulation run, typically 10 20 MC steps per particle, where a MC step is defined as the attempt to move one soft blob. The averaging leads to the possible appearance of non-integer numbers of patches per star. The width of the distribution of the average patchiness, as shown in figures 5 and 6, increases for higher functionality. However, the fluctuations about the mean number of patches are not influenced by the number of star in solution, demonstrating that the low density self-aggregated structures are indeed robust upon increasing density. For completeness, the average number of patches per star and their standard deviations for the different functionalities and asymmetries are listed in tables 1 and 2.
To further characterize the single star properties, we compute the average distance between the CM of each patch and the central, anchoring point of the chains of the corresponding star. We refer to this quantity as the average patch elongation L. The results obtained for L at five different functionalities and two asymmetry ratios are displayed in figures 7 and 8 as a function of the density. Not surprisingly, the stars experience a moderate compression upon increasing the density. Nevertheless, the average properties of the self-assembled patchy structure remain for practical purposes unaffected by the increase in concentration. It is, however, interesting to notice that for the higher asymmetry ratio α = 0.6 the stars tend to be considerably more rigid than the one with α = 0.4, which is demonstrated by the better preservation of the average elongation and width of the distribution upon increasing the density. This can be understood in terms of the fact that stars with bigger attractive patches are also more compact and thus less prone to deformations. As stars with α = 0.4 and 0.6 can be selected to have the same number of patches by an appropriate choice of the functionality, the change in rigidity with α is a property that could be exploited to produce self-assembling soft patchy particles of controllable flexibility. It would be interesting to see to what extent the increased rigidity of the stars result in a sturdier gel phase. On similar grounds, the rigidity of the stars grows with the functionality f ; yet, they remain soft, fluctuating and deformable, and thus they differ drastically from their hard-patchy colloids counterparts.
In figures 9 and 10 we analyze the orientational distribution of the patches around the center of the stars, by means of the relative angle ω formed between the directions in which patches are located with respect to the center of the corresponding star. Although the histograms show a wide breathing motion of the self-aggregated patches around the equilibrium angle, which tends to go from stretched for low functionalities to tetrahedral at the highest functionalities considered, it is quite striking to observe that the distributions are not significantly affected by the density. This is rather counterintuitive for a soft-system, as this implies that for TSPs the internal degrees of freedom of each star remain unaffected by large changes in density. The peak in the distributions P(ω) visible at small angles is caused by configurations where individual arms are fluctuating in and out of a particular patch and are therefore temporarily identified as separate nearby clusters.
As was mentioned before, upon increasing the density the stars form a gel, the aggregation of this phase takes place via the fusion of the patches of stars into larger communal patches shared by two or more stars. Accordingly, to obtain insight into the intermolecular orientation of the aggregating patches, we define another angle ψ. The latter now has its apex at the center of the common, fused patches shared by two stars and its sides extend from this point toward the centers of the two stars that participate into this physical association site. The analysis of the distribution of bond angles ψ between pairs of stars reveals once more an invariance of the angular distribution for all densities and star parameters. This can be seen in figures 11 and 12 where we plot the bond angle distribution P(ψ) for all the systems that were simulated. In all these cases, the distribution are very similar and are centered around ψ = 100 • . In summary, the combination of the topological and structural analysis of the gels formed by TSPs offers important information on the self-organization of the same. The gels are held together by the irreversible, chemical association sites of the individual stars, on the one hand, and reversible, physical association sites formed by fused patches between neighboring stars on the other. The latter are well-separated by one another by characteristic distances that can be discerned through the marked plateaus in the Euler characteristic of the amorphous gel. This morphological signature of the gel allows then for a clear, detailed structural analysis of the patchiness of each particle, as well as of the size-and orientational-properties of its constituent molecules, which all show a remarkable robustness to the concentration.

Structural analysis of the diamond crystal
In addition to the amorphous gel, TSPs at sufficiently high concentrations can give rise to ordered crystals, whose coordination is compatible with the patchiness p of the individual building blocks. This important property of the system has been demonstrated in [18]. We have already shown in section 4 that upon augmenting the density, stars maintain their average patchiness, the average patch elongation L and the average angle ω between patches. For such reasons TSPs can be thought of as soft patchy particles and their properties can be used to assemble crystalline phases. The mechanism that stabilizes the crystalline phase is the hierarchical self-assembly of the TSPs first into individual, soft-patchy particles, which thereafter behave as functionalized building blocks that allow for a mechanical stabilization of crystals with coordination numbers that are compatible with the patchiness of the TSPs. Molecules that have shown to spontaneously self-assemble at zero density in a tetrahedral phase (namely f = 15 and α = 0.4 and 0.6), have been pre-arranged in a simple cubic, bcc, fcc and diamond phases. All phases not compatible with the tetrahedral arrangement melt within the simulation time, while the diamond phase, whose coordination is compatible with the tetragonal arrangement of the four fold soft patchy particles, remains mechanically stable, resulting, e.g. in mean-square displacements of the star centers that saturate to a plateau. The same procedure performed with the other lattices, which do not possess the fourfold-coordination of the diamond Table 3. Properties of the TSPs in the diamond lattice for two different tetrahedrally coordinated molecules, namely f = 15, α = 0.4 and f = 15, α = 0.6 over a wide range of densities ρ/ρ * . The quantity p represents the average number of patches, ω the angle between the directions to CM of patches seen from anchoring point of the corresponding star and L the patch elongation (the distance from the anchoring point to the center of mass of the patch) measured in units a denoting the bond length in the underlying monomer-resolved system. crystal, leads to a rapid melting of the system, a property that highlights the selectivity of these soft, patchy particles in stabilizing distinct ordered structures. In a similar fashion, p = 6-TSPs stabilize a simple cubic lattice and all other attempted crystals quickly melt.
With the goal of getting insight into the (unusual, as it turns out) structure of the ordered phases and to the mechanisms leading to their stabilization, it is thus pertinent to perform a similar structural analysis of the periodic crystal. In this section, we fully characterize the single star properties within the diamond phase; we focus on the very same two asymmetry values that we analyzed in the gel phase, namely α = 0.4 and 0.6. We report in table 3 the average single star properties of stars with f = 15 and α = 0.4 for a wide range of densities, ρ/ρ * ∈ [0, 7.4], and for α = 0.6 in a range of densities ρ/ρ * ∈ [0, 10]. The table entries clearly demonstrate that the average patchiness p is a quantity that is extremely stable; in fact, within the crystal its value is almost identical to the one at infinite dilution, therefore spatial order reduces the fluctuations of the patchiness compared to those encountered in the gel phase. Stars that at zero density self-assemble in a fourfold structure maintain such a patchiness even when density is dramatically augmented. That this structural stability is not restricted to the patchiness only, is evidenced by average angle ω between the patches, the average patch elongation L, as well as the width of either distribution, which can all be read-off from table 3. It is important to note that this behavior extends over a wide range of densities, all the way from the infinite dilution limit of a single, isolated TSP to densities ten times the overlap concentration of the stars.
The full distributions of structural quantities of the stars within the diamond lattice are shown in figures 13 and 14 for the parameter combinations ( f = 15, α = 0.4) and ( f = 15, α = 0.6), respectively. In figures 13(a) and 14(a), the probability distribution of the previously defined interpatch angle ω is drawn. For both α = 0.4 and 0.6 the tetrahedral conformation of the soft patchy particle is confirmed by the average value of the angle, ω ∼ = 104 • ± 36 • between the patches, which remains quite unaffected by an increase of the density of stars in solution. The strong fluctuations around the mean are caused by connections that any given star has to not only its nearest neighbor but also to more distant ones, as will be shown below. Nevertheless, the distribution of the angle ω is narrower in the crystal than in the gel, since in the former the underlying tetrahedral orientation imposes stricter constraints on the orientational fluctuations. This robustness of the star conformations is further confirmed by the insensitivity of the distribution of the patch elongation, L, to the crystal density, shown in figures 13(b) and 14(b). In figures 13(c) and 14(c) we plot the distribution of the aforementioned bond angle ψ between star centers connected via a common patch, which is rather broad since the patches can bind with different neighboring stars at different densities. This feature is shown explicitly in figures 13(d) and 14(d), which demonstrate this very important result on the peculiar self-assembling mechanism of TSPs stabilizing the diamond crystal structure. There, the distributions P(s) of distance s between two star centers that share a patch is presented. The distributions show characteristic scales and are separated in bands that are peaked at locations compatible with successive coordination shells in the diamond lattice. Indeed, in the perfect diamond the distances s i between a given point and its first few, ith coordination shells are s 1 /d = 0.43, s 2 /d = 0.7 and s 3 /d = 0.82, where d denotes the size of the conventional, cubic unit cell of the diamond crystal. As the density grows, the distribution of the peaks shifts toward more distant neighboring shells, indicating that the crystal will respond to compression by simply reconnecting the network of patch bonds, compatible with the single star patch  figure 13 but for functionality f = 15 and asymmetry ratio α = 0.6. arrangement, to more distant stars in the crystal. For the particular case of tetrahedrally coordinated patches examined here, the diamond has the proper symmetries to allow for such a process to occur without frustrating the natural structure of the stars. This is the underlying reason for the fact the diamond was found to be mechanically stable for such a wide range of densities [18]. This particular mechanism can be at work only for soft and fluctuating patchy particles and it brings about the most distinct property of these building blocks that sets them apart from hard patchy colloids.
Although we restricted ourselves here to considering patches with tetrahedral symmetry that are able to stabilize the diamond structure, similar observations have also been made for the case of TSPs of functionality f = 20 and asymmetry α = 0.5 [18]. These particular stars exhibit a patchiness p = 6, which enables them to stabilize a simple cubic crystal in much the same fashion as the diamond structure described above. Also in this case, the single star properties are maintained on augmenting the density and the cubic crystal is found to be stable in a wide range of concentrations.

Summary and concluding remarks
We have presented extensive numerical evidence to set forth the proposal that TSPs constitute novel building blocks for the formation of well-controlled phases in soft condensed matter systems. Pretty much in the same fashion in which athermal star polymers bridge between polymer chains and hard colloids [26], TSPs play an analogous role in bridging the gap between block copolymer chains and hard patchy colloids [27]. At the same time, they seem to be having a decisive advantage with respect to the latter in terms of synthesis, since no cumbersome preparation techniques are necessary to form them [28,29]. Moreover, they are able to mechanically stabilize pre-arranged crystals for a very wide range of densities, by taking advantage of their soft character and by therefore gaining entropy through the possibilities to merge their patches with those of a large number of nearby neighbors [18]. Further, the stability of the structures with respect to polydispersity of the building blocks is guaranteed, since for broad distributions of f and α, TSPs self-organize to the same pattern of patches. Finally, our approach is robust and reliable since it has been based on a systematic coarse-graining that is back-tractable, and which does not contain any ad hoc parameters or functional forms for the effective interaction potentials [15,16,18].
There exist, at the same time, formidable challenges for the future. One category of open questions pertains to the possibility of further steering the behavior of the system by temperature changes. Indeed, in this case the interactions that involve the terminal units will change, and so will the patchiness of the particles as well. A systematic investigation is necessary to analyze how this will affect the properties of phases that have already been assembled at a given temperature [30]. A second, very important issue is that of spontaneous self-assembly into targeted ordered phases, without the artificial pre-arrangement on lattice sites invoked here for the purposes to establish the stability of the crystals. The spontaneous emergence of spatial order in the system may very well require the influence of suitably organized external agents, such as the presence of a patterned surface. Work along these lines is currently in progress.