A density-functional theory investigation of cluster formation in an effective-potential model of dendrimers

# A density-functional theory investigation of cluster formation in an effective-potential model of dendrimers

## Abstract

We consider a system of particles interacting via a purely repulsive, soft-core potential recently introduced to model effective pair interactions between dendrimers, which is expected to lead to the formation of crystals with multiple occupancy of the lattice sites. The phase diagram is investigated by density-functional theory (DFT) without making any a priori assumption on the functional form of the density profile or on the type of crystal lattice. As the average density is increased, the system displays first a transition from a fluid to a bcc phase, and subsequently to hcp and fcc phases. In the inhomogeneous region, the behavior is that found in previous investigations of this class of cluster-forming potentials. Specifically, the particles arrange into clusters strongly localized at the lattice sites, and the lattice constant depends very weakly on , leading to an occupancy number of the sites which is a nearly linear function of . These results are compared to those predicted by the more widespread approach, in which the DFT minimization is carried out by representing the density profile by a given functional form depending on few variational parameters. We find that for the model potential studied here, the latter approach recovers most of the predictions of the unconstrained minimization.

cluster formation, dendrimers, soft-core interactions, density-functional theory

## I Introduction

Bounded repulsive potentials arise in the description of effective, coarse-grained interactions between polymers. Because of their noncompact nature, mutual overlap between polymers does not require the actual superposition between the constituent monomers and is thus allowed, while still implying a steric repulsion due to the reduction of their available configurations. The prototypical example of these effective soft-core interactions is the repulsive Gaussian two-body potential, whose thorough study (1) has shown how its equation of state, correlations, and phase diagram differ markedly from those of hard-core fluids.

In recent years, much attention has been given to another class of soft-core repulsive interactions which occur, for instance, in the modelization of ring polymers (2) or dendrimers (3); (4). These interactions differ from the simple Gaussian potential mainly for being softer at small interparticle separation , possibly displaying even a minimum at . As a consequence, their Fourier transform is not positive everywhere, and has a negative minimum at non-vanishing wavevector  — a property sometimes referred to as characteristic of the class of repulsive interactions (5), which causes the tendency for the system to self-assemble into regular structures consisting of clusters of particles with a period , basically independent of density (6). This scenario is similar to that giving rise to microphases in hard-core fluids with competing attractive and repulsive interactions (7); (8), although the physical origin of clustering is different in the two cases. In fluids with competing interactions, clusters occur because on the one hand the attraction favors particle aggregation, while on the other hand the repulsion prevents the aggregates from growing beyond a certain size. In repulsive potentials, instead, clustering is due to the fact that particles are allowed to superimpose to each other at the price of a very small mutual repulsion.

The formation of cluster phases in potentials has been investigated in detail for the generalized exponential model (GEM) with interaction of the form (6); (9):

 vGEM(r)=ϵexp[−(r/R)4]. (1)

In this case one finds that, when the density is increased at sufficiently high temperature, the homogeneous fluid is replaced first by a cluster bcc phase, then by a cluster fcc phase at higher density. Below a triple-point temperature , the bcc phase disappears, and the only transition left is that between the fluid and the fcc cluster phase (9); (10) until, at even lower temperatures, the fcc phase splits into a sequence of fcc phases connected to each other by isostructural transitions, in which the lattice constant jumps discontinuously (11); (12).

Another family of potentials of this class has been proposed as a modelization of the effective interactions between amphiphilic dendrimers with a solvophobic core and a solvophilic shell (3); (4). In particular, for a specific choice of the parameters which determine the potentials between adjacent and non-adjacent monomers of the dendrimer, the interaction between two isolated dendrimers can be accurately represented by a fit to a sum of two Gaussians of opposite sign (4):

 βvD(r)=ϵ1exp[−(r/R1)2]−ϵ2exp[−(r/R2)2], (2)

where , being the Boltzmann constant, and , , , , being the diameter of the core monomers. A plot of this potential is displayed in Figure 1. Its Fourier transform in the neighborhood of the minimum is shown in the inset.

The main qualitative difference between and is that in the latter, corresponds to a minimum rather than a maximum, and the maximum is assumed for . This is due to the fact that at small separation between the dendrimers, the overlap between the solvophobic cores contrasts that between the solvophilic shells, thereby leading to a decrease in the repulsion with respect to that felt at intermediate separation when the cores do not overlap, but the shells still do.

Even though the formation of cluster crystals in dendrimers has already been addressed by computer simulation (13), to our knowledge the phase diagram of a system of effective particles interacting by the pair potential (2) has not been considered yet. In the present work, we address this problem by density-functional theory (DFT). We have used a simple, mean-field like free-energy functional, which has already been adopted many times to describe cluster formation in both the soft-core repulsion (6); (9); (14) and the competing interaction scenario (8). Here, we have given special attention to the procedure followed to minimize the functional. In fact, instead of describing the density profile by a given analytical expression depending on a small number of free parameters, we have not made any a priori hypothesis as to its form save for assuming that it is periodic, and the Bravais lattice vectors which determine the periodicity have themselves been obtained as a result of the minimization. Such a method is computationally more demanding than that more commonly adopted, but its unbiased nature is an asset, since it does not require any guess on the structures into which the system is most likely to assemble.

The plan of the paper is as follows: in Sec. II we briefly introduce the free-energy functional which we have used, and describe the method by which it has been minimized; in Sec. III we apply this method to the study of the phase diagram of a system of dendrimers, described as effective particles interacting via the potential (2); in Sec. IV, we try to get further insight into the results discussed in Sec. III, by assessing up to which extent they can be recovered by adopting an analytical representation of the density profile; finally, in Sec. V we summarize our results and draw our conclusions.

## Ii Theory

In order to study the occurrence and structure of the inhomogeneous phases of a system of particles interacting via the pair potential (2), we have used DFT. According to DFT, at a given temperature , chemical potential , and volume , the density profile is obtained by minimizing the grand potential functional , whose value then yields the grand potential , where is the pressure. Henceforth, the density profile will be denoted by , while with no point-dependence indicated will refer to the average density . For a system of particles interacting via a pair potential , the exact grand potential functional is of course unknown, and many approximate formulations of DFT have been suggested. Here we have used the following simple form for :

 βΩ=∫d3rρ(r){ln[ρ(r)Λ3]−1−βμ}+12∫d3r∫d3r′ρ(r)ρ(r′)βv(r−r′), (3)

where is the thermal wavelength of a particle of mass . In the rest of the paper will be dropped, since its actual value does not affect any of the results presented here. The first term in the rhs of Eq. (3) is the exact ideal-gas contribution, while the interaction is approximately taken into account in the second, “excess” term. This expression can be considered as a generalization to inhomogeneous systems of the van der Waals mean-field approach, which is in fact obtained from Eq. (3) in the special case of a uniform density profile . A similar approximation for could be adopted also in a system of non-overlapping particles whose pair interaction contains a singular repulsive part, but this would require as a preliminary step to treat separately the singular contribution. Here, instead, coincides with the total soft-core interaction of Eq. (2).

For a homogeneous state, the stability condition with respect to a small perturbation of the density

 ∫d3r∫d3r′δ2(βΩ)δρ(r)δρ(r′)∣∣∣ρδρ(r)δρ(r′)>0 (4)

is equivalent to the requirement for every wave vector , where is the structure factor. Equation (3) gives the random-phase approximation (RPA) for :

 S(k)=[1+ρβ~v(k)]−1, (5)

where is the Fourier transform of . If belongs to the class, the above condition on is violated for , where is the wave vector at which assumes its negative minimum. The boundary between the region where the homogeneous state is stable and that where it is not corresponds to the condition of marginal stability , which identifies a locus in the - plane generally referred to as the -line. Along the -line, the RPA diverges at , thus signalling the tendency of the system to form density modulations with a characteristic length . This can be compared with the behavior on the spinodal line predicted by the RPA for systems which undergo a liquid-vapor transition: in the latter case, the divergence of occurs for , and the homogeneous phase is unstable with respect to separation into two bulk phases.

In the present case in which the potential is athermal, the -line is actually a line at constant . For the potential (2) one has , , so that . However, it should be recalled that, just like the spinodal line for the liquid-vapor transition, the -line does not generally coincide with the boundary between the homogeneous and the inhomogeneous fluid. The latter must be obtained by comparing the free energies of the two phases, and may (and usually will) lie in a region where the homogeneous fluid is still mechanically stable.

In order to obtain further insight in the phase diagram and the structure of the inhomogeneous phases, it is necessary to turn to the minimization of the functional (3). As stated in Sec. I, we have assumed from the outset that the density profile is periodic, i.e.

 ρ(r+ai)=ρ(r), (6)

where , , , , are a set of vectors which define a Bravais lattice. Therefore, can be expanded in a Fourier series:

 ρ(r)=1v∑me−ikm⋅x^ρm, (7)

where is the volume of the unit cell, is a vector of the reciprocal lattice, and denotes a set of three indexes . The expansion coefficients are given by:

 ^ρm=∫vd3reikm⋅rρ(r). (8)

By use of Eqs. (6) and (7), the functional (3) can be rewritten in the following form:

 Missing or unrecognized delimiter for \left (9)

We observe that, both in Eq. (8) and (9), the integration in is restricted to the unit cell. Hence, we may set , where is the matrix obtained by arranging the lattice vectors into columns, and is a vector whose components vary in the interval . By doing so, it is readily seen that in Eq. (9) neither the ideal-gas term nor the Fourier components of the density profile that appear in the excess term depend on the specific kind of lattice: that is, these quantities are determined solely by the values assumed by the density profile in the unit cell, irrespective of the cell geometry. The information on the lattice enters in Eq. (9) only via the reciprocal lattice vectors at which is evaluated. This feature makes it easy to implement a numerical procedure, in which the optimization of the grand potential functional (9) is performed with respect to both and the cell geometry, i.e., the elements of . In this study, we have assumed that the vectors of the unit cell are mutually orthogonal, so that is diagonal with eigenvalues , and the reciprocal lattice vectors have the form . This assumption simplifies the calculation, but it could be released without introducing any conceptually new element. To perform the minimization, one has to solve the Euler-Lagrange equations as well as , where the functional derivative with respect to and the partial derivative with respect to are given by:

 δδρ(r)(βΩV)=1v{ln[ρ(r)]−βμ}+1v2∑me−ikm⋅r^ρmβ~v(km), (10) ∂∂hi(βΩV)=1v2∑m|^ρm|2dβ~vd(k2)(km)him2i. (11)

In the numerical solution, the functional was first discretized by sampling on a finite set of points , so as the replace the functional derivatives with the partial derivatives with respect to . The minimization was then carried out by an iterative algorithm based on the steepest descent. In the basic version of the steepest descent, and are updated recursively by moving “downhill” in the direction opposite to that of the gradient of the discretized functional :

 ρk+1n=ρkn−η∂∂ρn(βΩDV)∣∣∣k, (12) hk+1i=hki−θ∂∂hi(βΩDV)∣∣∣k, (13)

where is the iteration index, and , are the parameters which determine the size of the downhill step. In order to increase its efficiency, the above algorithm was improved by introducing preconditioning and conjugate gradients in Eq. (12), and by determining the step-size parameters , adaptively at each iteration. A detailed description of these technical features will be given elsewhere (15).

The discretization of the density profile inside the unit cell was performed on points. The cell in real space was initially chosen as a cube with edge lengh , and was then evolved according to Eq. (13). The trial density profile used to start the minimization at a given thermodynamic state was set either to a random noise superimposed to a uniform density, or to a sinusoidal modulation, or to the equilibrium of a nearby state. The need of performing different minimization runs for the same state is due to the fact that soft-core potentials such as that considered here may present many different inhomogeneous structures, each of which corresponds to a local minimum of the grand potential. Hence, the convergence of the algorithm does not guarantee per se that the global minimum has been identified: the search for the latter necessarily involves a comparison between different local minima of the grand potential. Of course, this problem is intrinsic to systems with a complex energy landscape, and many sophisticated strategies have been developed to tackle it (16); (17). However, at the present stage we have not tried to implement any of them, and have just been content with our rather naive approach of performing different runs by changing .

## Iii Results

By minimizing the grand potential (3) via the algorithm described above it is found that, as expected, the system becomes inhomogeneous at sufficiently high density . The phase diagram which we have obtained is shown in Figure 2. Since the potential is athermal, the temperature does not play any role in the phase behavior, which is then determined by the density alone.

As one moves away from the homogeneous region by increasing , the fluid freezes by forming a cluster bcc phase. As anticipated in Sec. II, the -line actually lies inside the inhomogeneous region. The freezing and melting lines have been determined by performing a Maxwell construction, and the corresponding coexistence region is a rather narrow stripe centered at . A calculation based on the same free-energy functional used here gives for the density of the fluid-solid transition of potentials the expression (6):

 ρ=[1.393β|~v(km)|]−1. (14)

By substituting the above-mentioned value , this gives , in close agreement with the present result.

If is further increased, two more transitions are found, first from the bcc to a hcp phase, then from the hcp to a fcc phase. In both cases the densities of the coexisting phases are very close to each other, so that the width of the coexistence region cannot be resolved on the scale of Figure 2. The transition from a bcc lattice to a close-packed lattice on increasing is in agreement with the general argument put forth in (6), as well as with the phase diagram of the GEM potential (1) obtained in (9), although in that case the hcp phase was not observed.

The density profile at a density in the bcc region near the fluid-solid boundary is shown in the upper panel of Figure 3 along the direction connecting nearest neighbors as a function of the distance from a given lattice site.

In order to assess the spatial extension of the density peaks against that of the interaction potential, the latter has also been plotted in the figure. We observe that, despite the proximity to the fluid region, is already strongly localized at the lattice sites. It is interesting to compare this density profile with that of the GEM potential at a similar density. To this purpose, we have considered the interaction (1) at a reduced temperature such that , somewhere between the well and the hump of , see Figure 1. The width of has then been fixed at the value , for which the integrated intensity of is the same as that of . The lower panel of Figure 3 shows obtained by the numerical minimization of the grand potential functional (3) for the GEM potential at , as well as for the parameters given above. The most stable lattice is again the bcc, but the peaks are significantly lower. Since the density is about the same in both panels, this must be compensated by a larger width of the peaks. That this is indeed the case is shown more clearly in Figure 4, where the logarithm of the two density profiles of Figure 3 has been plotted as a function of .

Hence, the dendrimer potential leads to a higher localization with respect to the GEM potential, and it is natural to trace this back to the minimum of at . We also note that on the scale of Figure 4, a Gaussian corresponds to a straight line. It appears that, for both and , can be accurately represented by a Gaussian as long as one has , as one would expect from Figure 3. Deviations from the Gaussian profile become evident at smaller .

As one moves more deeply into the inhomogeneous region, the height of the peaks of the density profile increases very rapidly, as shown in Figure 5 which displays in the fcc domain for and . These have also been plotted in Figure 4 in order to show that the peaks become narrower and narrower on increasing .

On the other hand, by comparing Figures 3 and 5 it appears that the nearest-neighbor distance is hardly affected by density. In fact, this behavior is peculiar of crystals formed by potentials for which, even though the specific kind of periodicity of depends of course on the lattice formed by the system, its spatial extent is basically determined only by the density-independent wave vector of the minimum of the potential in Fourier space (6). At variance with atomic crystals, this leads to a nearest-neighbor distance and a volume per lattice site which are nearly unaffected by . Since the number of particles in a cluster located at a lattice site is related to the average density by the relation (6) , a nearly constant implies a nearly linear dependence of on .

The nearest-neighbor distance and the number of particles per cluster are shown in Figure 6 for the fcc phase. Although not strictly constant, is indeed found to be a very slowly varying function of density, giving the anticipated linear dependence of on . It is interesting to observe that, as is increased, is actually found to increase, albeit slightly so. At first this may appear at odds with the fact that, if the particle number is kept constant and is increased, the volume of the system must obviously decrease, but in fact no contradiction arises, since the slight increase in is more than compensated by the decrease in the number of lattice sites . In other words, a cluster crystal manages to change its volume not by changing the lattice spacing, but rather the number of lattice sites, and consequently their population . In this respect, whether increases or decreases on increasing is of little relevance, provided it remains a slowly varying function of .

## Iv Analytical representation of ρ(r)

As evident from Figures 3 and 5, as increases the density profile becomes very strongly localized at the lattice sites. It is then natural to represent it by a sum of non-overlapping or weakly overlapping peaks centered at these sites. The crudest approximation of this kind consists in assuming the highest degree of localization, and modeling as a sum of -functions:

 ρδ(r)=nc∑mδ(r−Rm), (15)

where describe the positions of the lattice sites. Strictly speaking, Eq. (15) cannot be true, since it would imply an infinite free energy cost due to the loss of translational entropy. In a system with temperature-independent interactions, Eq. (15) would hold at , but then could not be chosen at one’s will, since it would be bounded to assume integer values (11). However, the above assumption still makes sense as an approximation, if one assumes that, for the true , the free energy is dominated by its energetic contribution and may thus be approximated by disregarding the entropic term altogether. One is then left with an internal energy per particle which, for of Eq. (15), has the form:

 (βEN)δ=12ζd3ρ[βvD(0)+∑m≠0qmβvD(cmd)]−12βvD(0). (16)

In the above expression, represents the interaction between particles belonging to the same cluster, and the interaction between different clusters is given by the lattice sum of the potential over successive neighbor shells. The constants and are specific to the lattice considered, and represent respectively the number of -th neighbors and their distance from a given lattice site in units of the nearest-neighbor distance . The number of particles per cluster has been expressed as , where the quantity is the volume per lattice site expressed as a function of , and again depends on the lattice: in particular, for the bcc lattice and for the fcc and hcp lattices. If Eq. (16) is truncated at the nearest-neighbor level, one obtains:

 (βEN)δ≃12ζd3ρ[βvD(0)+qβvD(d)], (17)

where is the number of nearest neighbors. Equation (17) can be regarded as a simple variational expression to be minimized at fixed with respect to the single variational parameter . Again, we must observe that this statement taken at face value makes little sense, since the minimum of Eq. (17) at fixed is actually obtained for : in other words, this expression does not prevent the unphysical collapse of the crystal at a single lattice point. What drives the system away from such a collapse is the fact that, as is decreased, more and more shells contribute significantly to the energy. This leads to a divergence of the sum in Eq. (16), but such an occurrence is lost in Eq. (17), which contains just one shell. Nevertheless, Figures 3 and 5 indicate that, at the actual , the interaction energy between different lattice sites is rapidly decreasing with their distance, so that the contribution due to the outer shells should indeed be very small. Hence, although Eq. (17) cannot be trusted at every , it might still be a reasonable approximation in the neighborhood of the value of chosen by the system, which would then be recovered as a local minimum of Eq. (17). In fact, besides the spurious minimum at , this expression is found to have a (spurious) local maximum and a local minimum. The latter is equal to for the bcc lattice and for the fcc and hcp lattices. In Figure 6 we have reported together with the value obtained by identifying the wave vector of the minimum of with the nearest-neighbor distance in reciprocal space, which yields . One sees that is a rather accurate estimate of the nearest-neighbor distance, actually more so than . As a consequence, also the slope of the vs. plot is accurately reproduced.

The value of obtained from the minimization of is clearly independent of density, and this would be true even if the summation over all neighbor shells of Eq. (16) were used instead of the truncated expression of Eq. (17). In order to describe the dependence of on density, the contribution of the translational entropy to the free energy must be taken into account, and this in turn requires to go beyond the representation of as a superposition of -like spikes. Following Refs. (9) and (6), we then replace the -functions in Eq. (16) with Gaussians of finite amplitude:

 ρG(r)=nc(απ)3/2∑me−α(r−Rm)2, (18)

where the normalization constant has been determined so that gives the number of particles at each site as requested. If is substituted into Eq. (3) and the overlap between the Gaussians at different lattice sites is neglected in the ideal-gas term, then the following expression for the Helmholtz free energy is obtained (9):

 (βFN)G=lnnc+32lnαπ−52+nc√α32π∫+∞0drr2e−αr2/2βv(r) +nc√α8π∑m≠01Rm∫+∞0drr[e−α(r−Rm)2/2−e−α(r+Rm)2/2]βv(r). (19)

For the potential , the integrals in Eq. (19) can be performed analytically, and one has:

 (βFN)G=ln(ζd3ρ)+32lnαπ−52+12ζd3ρ⎧⎨⎩ϵ1(αα+γ)3/2{1+∑m≠0qme−αγc2md2/[2(α+γ)]} −ϵ2(αα+δ)3/2{1+∑m≠0qme−αδc2md2/[2(α+δ)]}⎫⎬⎭, (20)

where , , and we have expressed as as before. In the limit , the interaction term in Eq. (20) reproduces Eq. (16), save for the thermodynamically irrelevant term . However, the minimization of the free energy must now be carried out with respect to both and . This requires the solution of the non-linear system of equations . If we make the change of variables:

 x = γd22(1+γ/α), y = δd22(1+δ/α), (21)

the above system can be cast in the rather symmetric form:

 x = HγδyS(y)HS(y)−3, y = KδγxS(x)KS(x)+3, (22)

where we have set:

 H = √2ζρϵ2δ3/2, K = √2ζρϵ1γ3/2, (23)

and the function is defined by:

 S(x)=x3/2⎡⎣3+∑m≠0qm(3−2c2mx)e−c2mx⎤⎦. (24)

We can now introduce the same approximation used in going from Eq. (16) to Eq. (17) by truncating at the nearest-neighbor shell:

 S(x)≃x3/2[3+q(3−2x)e−x]. (25)

As observed above for Eq. (17), strictly speaking such a truncation would lead to the collapse of the crystal at , but one expects that a physical solution may still be recovered as a local minimum of the free energy thus obtained. The results for and given by the numerical solution of Eq. (22) with the approximation (25) for the fcc lattice are again displayed in Figure 6. Both and run very close to the corresponding quantities obtained from the fully numerical minimization of the free energy described in Sec. II. It appears then that, provided the finite amplitude of the density peaks is taken into account, keeping the interactions only up to nearest neighbors is already sufficient to reproduce accurately the density dependence of , at least at high density where the fcc phase is favored. For the bcc phase, it turns out that also the second-nearest neighbor shell must be included in the calculation. This is due to the fact that in the bcc lattice the nearest and second-nearest neighbor shells are much closer to each other than in the fcc lattice. In either the fcc or the bcc, inclusion of the other neighbor shells from the third on leads to a result which is indistinguishable from that given by the fully numerical calculation, as again shown in the figure for the fcc phase.

In view of such an agreement, it is natural to ask how the phase behavior predicted by Eq. (20) on the basis of the Gaussian parameterization of the density profile (18) compares with that discussed at the beginning of Sec. III. This is a rather subtle issue, because the free energies of the different phases are very close to one another. Therefore, even a small change in their values may alter their relative balance. In fact, it is found that the truncation of the interactions at the nearest-neighbor level in Eq. (20) has rather important consequences on the phase diagram, because it would make the bcc phase more stable than both the fcc and hcp phases even at high density, thus leaving the bcc as the only inhomogeneous phase. In order to recover the stability of the close-packed structures at high density, in Eq. (20) the interactions have to be taken into account at least up to the second-nearest neighbor shell. As recalled above, the contribution of this shell is much more substantial in the bcc lattice than in either the fcc or hcp, and proves crucial in order to move the bcc free energy above that of the close-packed lattices. In Figure 7 we show the free energy per particle and the grand potential per unit volume of the bcc, fcc and hcp phases obtained via the minimization of Eq. (20) by including a rather large number of neighbor shells such that, for the densities considered, the addition of further shells would leave the free energy unchanged within numerical accuracy.

These results are compared to those given for the same phases by the fully numerical minimization. It is evident that the two sets of data are indistinguishable on the scale of the figure. Moreover, for both sets of data the free energies and grand potentials corresponding to different phases are also indistinguishable, although by enlarging the energy scale the small differences which determine the most stable phase at a certain density would be appreciated. In particular, Eq. (20) predicts that the bcc lattice becomes unfavorite with respect to close-packed lattices for , in very good agreement with the position of the bcc-hcp transition shown in Figure 2.

There is, however, a difference between the phase behavior given by the unconstrained minimization discussed in Sec. III and that given by Eq. (20), namely, while the former predicts the occurrence of the hcp phase between the bcc and fcc phases, according to the latter the free energy of the hcp phase is always higher than that of the fcc phase. Hence, the hcp phase is not observed, and as the density is increased the sequence of phase transitions is fluid–bcc–fcc instead of fluid–bcc–hcp–fcc as in Figure 2. Nevertheless, it should be remarked that in both cases the free energies of the hcp and fcc phases are everywhere extremely close to each other: specifically, their relative difference is of the order of , while that between the bcc and fcc phases is . This is not surprising, in the light of the fact that the free energies of the fcc and hcp lattices become distinct only from the third-neighbor shell on. Therefore, even though the numerical accuracy of the calculation is high enough to allow resolution of such tiny differences, speaking of a “most stable” phase in this context is somewhat academic, and one should more realistically regard the fcc and hcp phases as nearly degenerate. One also expects that whether the contributions to the free energy from the high-order shells will stabilize one phase or the other will most likely depend on rather specific features of the decay of the potential at long distance.

## V Conclusions

We have employed density-functional theory to study the phase diagram and the density profile of a system of particles interacting via an athermal, penetrable-core potential aimed at modeling the effective pair interaction between dendrimers (4). The potential belongs to the so-called class (5), which is expected to lead to the formation of cluster crystals at high enough density. The results discussed in the paper confirm this prediction, and the main features of the behavior of the system are those already encountered in the study of the generalized exponential model potential carried out some years ago (6); (9). Specifically: i) in the inhomogeneous phases, the particles arrange into clusters localized at the lattice sites, such that each cluster contains many particles. ii) The lattice spacing is nearly independent of density, which in turn implies that the number of particles per cluster is a nearly linear function of density. In the present system the localization is even stronger than that observed in the GEM potential, because of the local minimum presented by at . iii) As the density is increased, one observes a transition from the bcc structure to close-packed structures. Here we find both hcp and fcc phases, while for the GEM potential the hcp phase was not observed.

The density functional which we have adopted has the same mean-field form already used in the study of the GEM potential (6); (9). The minimization of the functional has been performed by a parameter-free procedure, in which both the values assumed by the density profile and the lattice structure were determined as the outcome of the calculation. The only assumptions which we made are that the density profile is periodic, and that the unit cell can be represented by a set of mutually orthogonal vectors. Although not absolutely general, this assumption does not necessarily rule out the occurrence of structures whose primitive cell cannot be represented in this way, e.g. the aforementioned hcp lattice. The results thus obtained have then been compared with those given by the more common procedure, in which one considers a priori a set of possible lattice structures, represents as a superposition of Gaussians centered at the lattice sites, and identifies the most stable structure among the selected ones by inspection of their free energies once the minimization has been carried out. The comparison shows that the Gaussian parameterization reproduces most of the results obtained by the unconstrained minimization. In fact, for the densities at which the system arranges into a close-packed lattice, its properties can be already described very satisfactorily by further simplifying the treatment, and disregarding the contribution to the excess free energy beyond nearest-neighbor lattice sites, although this contribution is necessary in order to account for the transition from the bcc to the fcc structure. The only qualitative prediction which we did not recover by the Gaussian representation of is the occurrence of the hcp phase between the bcc and fcc, the fcc phase being always favored with respect to the hcp. However, in both approaches the relative differences in free energy between the two structures are extremely small, so that labelling either of them as the more stable one is of little practical relevance. In fact, these differences are due to the interactions from the third nearest-neighbor shell on, which we expect to depend on rather specific details of the decay of the interaction at long distance. In the light of these considerations, for the system considered here the unconstrained minimization may admittedly appear as an overkill. One can easily imagine, however, that this may not be the case in other situations, in which the density profile is less amenable to being described by a superposition of Gaussians. We plan to report on such an instance in the near future (15).

Before concluding this discussion, one more consideration is in order: since the potential has been introduced as a representation of pair interactions between dendrimers, one may ask whether real dendrimers can be expected to display the behavior described above, especially as far as cluster formation is concerned. This depends on the actual relevance of many-body interactions, which are disregarded in the present picture. The very fact that, according to the density profile predicted on the basis of , most particles concentrate at the lattice sites leading to very high values of the local density, suggests that many-body effects are indeed relevant, and may significantly influence the behavior of the system. This conclusion was reached in a computer simulation study of an assembly of ring polymers (2). For that system, the pair potential between two isolated polymers displays a profile similar to that of . In particular, it also belongs to the class, so that it is expected to promote cluster formation. However, it was found that at finite density many-body effects lead to an effective pair interaction which is substantially modified with respect to its zero-density limit, to the point that clustering does not take place. In such a situation, the pair interaction between two isolated ring polymers is clearly of little use in understanding the behavior of the system. On the other hand, a subsequent study showed that in dendrimers the formation of cluster crystals does take place (13), although many-body effects are still very important. In fact, they prevent the number of dendrimers per lattice site from growing indefinitely, as is instead predicted by the two-body potential, and inhibit the formation of the crystal at high density, where a percolated network of dendrimers is favored. Moreover, in line with this tendency of many-body interactions to curb the growth of with , the lattice spacing was found to be a slowly decreasing function of , rather than a slowly increasing one as in the present study. It should also be observed that the parameters of the interactions between monomers used in (13) are different from those used in (4), and the resulting pair potential for isolated dendrimers leads to clustering at a lower density than that obtained with the pair potential (2) considered here. Nevertheless, it is reassuring to have evidence that for this system the formation of cluster crystals does not rest solely upon the effective pair potential picture.

## Acknowledgements

This paper is dedicated to Johan Høye on his 70 birthday. The author wishes to thank Alberto Parola, Luciano Reatto, and Christos Likos for illuminating conversation.

### References

1. F. H. Stillinger and D. K. Stillinger, Physica A 244, 358 (1977); A. Lang, C. N. Likos, M. Watzlawek and H. Löwen, J. Phys.: Condens. Matter 12, 5087 (2000); C. N. Likos, Phys. Rep. 348, 267 (2001).
2. A. Narros, A. J. Moreno, and C. N. Likos, Soft Matter 6, 2435 (2010).
3. B. M. Mladek, G. Kahl, and C. N. Likos, Phys. Rev. Lett. 100, 028301 (2008).
4. D. A. Lenz, B. M. Mladek, C. N. Likos, G. Kahl, and R. Blaak, J. Phys. Chem. B 115, 7218 (2011).
5. C. N. Likos, A. Lang, M. Watzlawek, and H. Löwen H Phys. Rev. E 63, 031206 (2001).
6. C. N. Likos, B. M. Mladek, D. Gottwald, and G. Kahl, J. Chem. Phys. 126, 224502 (2007).
7. R. P. Sear, S.-W. Chung, G. Markovich, W. M. Gelabrt, and J. R. Heath, Phys. Rev. E 59, R6255 (1999); R. P. Sear and W. M. Gelbart, J. Chem. Phys. 110, 4582 (1999); A. Imperio and L. Reatto, J. Phys.: Condens. Matter 16, S3769 (2004); A. Imperio and L. Reatto, J. Chem. Phys. 124, 164712 (2006); A. Imperio and L. Reatto, Phys. Rev. E 76, 040402 (2007); A. J. Archer A J and N. B. Wilding, Phys. Rev. E 76, 031501 (2007);
8. A. J. Archer, Phys. Rev. E 78, 031402 (2008).
9. B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and C. N. Likos, Phys. Rev. Lett 96, 045701 (2006).
10. B. M. Mladek, P. Charbonneau, and D. Frenkel, Phys. Rev. Lett. 99, 235702 (2007).
11. T. Neuhaus and C. N. Likos, J. Phys.: Condens. Matter 23, 234112 (2011).
12. K. Zhang, P. Charbonneau, and B. M. Mladek, Phys. Rev. Lett. 105, 245701 (2010); N. B. Wilding and P. Sollich, EPL 101,10004 (2013).
13. D. A. Lenz, R. Blaak, C. N. Likos, and B. M. Mladek, Phys. Rev. Lett. 109, 228301 (2012).
14. B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and C. N. Likos, J. Chem. Phys. B 111, 12799 (2007).
15. D. Pini, A. Parola, and L. Reatto, in preparation.
16. D. Gottwald, G. Kahl, and C. N. Likos, J. Chem. Phys. 122, 204503 (2005).
17. C. J. Pickard and R. J. Needs, J. Phys.: Condens. Matter 23, 053201 (2011).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters