Free energy calculations for rings and chains formed by dipolar hard spheres

We employ a method based on Monte Carlo grand-canonical simulations to precisely calculate partition functions of non-interacting chains and rings formed by dipolar hard spheres (DHS) at low temperature. The extended low temperature region oﬀered by such cluster calculations, compared to what had been previously achieved with standard simulations, opens up the possibility of exploring a part of the DHS phase diagram which was inaccessible before. The reported results oﬀer the unique opportunity of verifying well-established theoretical models based on the ideal gas of cluster approximation in order to clarify their range of validity. They also provide the basis for future studies in which cluster–cluster interactions will be included.


Introduction
][24] Single-domain magnetic nanoparticles, having an intrinsic magnetic moment, effectively behave like nanomagnets.With sizes ranging from 10 to 50 nm, these magnetic nanocolloids, once in solution, experience brownian motion.In addition, provided the magnetic crystallographic anisotropy of a nanoparticle is energetically strong enough, the rotation of its magnetic moment is coupled to the rotation of the particle as a rigid body. 25he main interaction between the dipoles is a long-range anisotropic dipolar interaction.This interaction determines a preferred head-to-tail orientation of dipole moments.[28][29][30][31] However, the self-assembly scenario is not limited to linear chains: magnetic particles can also form ring structures 32,33 and networks. 34,35Experimentally, all these clusters have been observed in solution using cryogenic electron microscopy 36 and atomic force microscopy of assemblies at cross-linkable oil-water interfaces. 37arious simulation studies described self-assembly in magnetic nanocolloids by using simple models of dipolar hardor soft-spheres. 34,38,396][47][48][49] However, recent simulation studies 32 suggest that dipolar interactions alone are not able to sustain a gas-liquid like phase separation, at odds with the pioneering predictions of de Gennes and Pincus. 50ecently, we showed that there is a structural chain-to-ring transition in a highly diluted gas of DHS as temperature decreases 51 and extended this study to higher concentrations where various branched structures were also considered. 52,53hese investigations provided a possible solution to a longlasting debate about the non-monotonic temperature dependence of the initial magnetic susceptibility in suspensions of magnetic nanoparticles. 54,55These conclusions were based on results of Monte Carlo simulations and density functional theory.However, the main ingredients of the system free energy -the partition functions of various clusters -have never been extracted from the simulations and were only calculated approximately, albeit rather carefully, in analytical models.
In the present study we focus on providing a precise evaluation of the partition functions of DHS chains and rings as a function of their sizes.In order to do so, we first build on the simulation method proposed in ref. 56-58.This method is based on constrained grand canonical Monte Carlo single-cluster simulations.Secondly, we compare the available analytical expressions for the same quantities to the simulation data.Quite remarkably, the theoretical parameter-free expressions quantitatively reproduce the numerical data at low density.
The manuscript is organised as follows.In the next section we thoroughly describe our simulation approach.Next, in Section 3, we verify previously used analytical expressions for partition functions, and the density functional approach describing DHS self-assembly.The last section contains a summary of the work and some perspectives.

Simulation method
DHS systems are characterised by two interactions: the anisotropic dipole-dipole attraction and the steric repulsion.The anisotropic dipole-dipole interaction has the form: where m k is the magnetic moment of particle k; |r ij | is the distance between particles i and j and r ˆij = r ij /|r ij |.In the present study we assume that all the magnetic moments have the same magnitude m.The steric repulsion has a standard hardsphere form: where s is the particle diameter.We define the cluster partition function Z N of a cluster made of N indistinguishable DHS particles as where b is the inverse thermal energy and L 3 ¼ l T 3 Ð dO 1 is the normalised thermal wavelength l T (arising from the integration over the translational and rotational momenta), convenient to define a spherically-averaged partition function.Here r i and O i denote the position and orientation of the i-th particle, respectively.The superscript 0 indicates that the sum is restricted to all phase space configurations for which the N particles form a cluster.A cluster is commonly defined as a collection of particles in which all possible pairs are connected via a sequence of bonds.Consequently, a preliminary step in every cluster calculation is the operative definition of bond, which will be given below.In the following we choose L 3 = s 3 as the thermal volume.
Analytically, one can easily calculate the expressions for and The ratio between the two has the form: The evaluation of the cluster partition function (or cluster free energy) requires an average of the system Hamiltonian over all particle positions r i and dipolar orientations O i for which the cluster constraint is satisfied.A numerical evaluation of the free-energy requires in general a reversible path from an analytically-known reference system to the required thermodynamic state.In the case of a cluster, it is possible to define such a path connecting the required cluster size to a monomer (a cluster of size one), for which the partition function is analytically known (eqn (4)), via a path of clusters of all intermediate sizes.This idea, presented in ref. 56, and then optimised in later publications, 57,58 requires a grand-canonical (GC) simulation (e.g. a simulation performed at fixed temperature T, volume V and activity z e bm , with m being the chemical potential) in which the particle additions and deletions are limited to configurations in which all particles in the system are part of a single cluster.
In the GC ensemble, the probability PðNÞ of observing a cluster of size N is proportional to where the GC partition function Z GC is defined as Therefore the ratio eliminates the unknown GC partition function and allows one to precisely evaluate the canonical cluster partition function as It is important to recall that Z N /Z 1 is independent of z.The value of z used in the GC calculation is thus arbitrary.A convenient choice which reduces the numerical error is the one for which PðNÞ Pð1Þ is of order one for all N.
Our operative definition of a bond, initially proposed in ref. 32

, employs a simple distance-energy criterion: two particles
This journal is © The Royal Society of Chemistry 2017 i and j are bonded if the interparticle distance r ij o 1.3s and if the interparticle potential energy V dd (r ij ) o 0. The 1.3s cut-off roughly coincides with the position of the first minimum in the radial distribution function, and hence selects particles belonging to the first neighbour shell.0][61] However, we do not expect the results to be sensitive to the precise definition of a bond at the values of T and r investigated here.
We run GC simulations with four types of MC moves: (i) the roto-translation move in which a random selected particle is displaced randomly by at most AE0.05s and simultaneously randomly rotated by at most AE0.1 rad; (ii) the pivot move, in which a random particle i and one of its bonded neighbours j are selected; if the deletion of the bond splits the cluster in two pieces, then one of the two pieces is rotated by a random angle around the dipole direction of particle i.If new bonds are generated after the rotation, the move is rejected; (iii) the insertion of one particle inside a shell of internal radius s and width 0.3s (corresponding to an insertion volume ) around a randomly selected particle; (iv) the deletion of one randomly selected particle.In all cases, the move is rejected if the final configuration is not composed by a single cluster.We find that the pivot move is particularly convenient to explore the low-T and low-r region in equilibrium and to sample the conformations of the self-assembled structures.This move -originally introduced by Lal 62 and independently suggested by MacDonald et al. 63 and by Mandras and Sokal 64 turns out to be highly efficient for studying configurational properties such as the physical size of the clusters.We note that no periodic boundary conditions are present (the volume is ideally infinite) and the potential energy of the system is composed by the sum of all N(N À 1)/2 interactions.The acceptance rules for inserting and deleting a particle have the form: 58 where n ins counts the number of bonds formed by the added particle, to compensate for phase-space overcounting (the same phase-space point could be obtained by inserting the particle in the bonding shell of any of the particle with which the inserted particle can form a bond) and DE is the change in potential energy between the final and the initial state.In the simulation, z is optimized to sample with the same statistics all values of N.
To speed-up the calculation and to benefit from multiprocessor clusters, the grand-canonical simulation is parallelised through a successive umbrella sampling technique, 65 where simulations exploring different but overlapping N intervals are run on different processors.In this study, the largest cluster investigated is composed by 200 monomers.The knowledge of the cluster partition function makes it possible to evaluate the cluster size distribution at any density under the hypothesis of an ideal gas of equilibrium clusters.Indeed, by minimizing the free energy of a gas of noninteracting clusters with respect to the cluster size distribution n(N), one obtains 66,67 Here we have used as normalization of the n(N) function P N NnðNÞ=V ¼ r, where r is the dimensionless particle number density (multiplied by s 3 ) and r 1 nð1Þ V .It is also possible to take into account the excluded volume interaction in an approximate fashion by selecting the Carnahan-Starling expression for the free energy of hard spheres as the reference free energy of the fully monomeric state by writing 58 where During the course of each GC simulation we keep track of the number of configurations in which the cluster is a ring (configurations in which all particles have two bonds) or a chain (configurations in which N À 2 particles have two bonds and two particles have a single bond).The remaining configurations are classified as branched structures.Since the partition function Z N of a cluster of size N can be written as sum of the partition functions of rings Z ring N , chains, Z chain

Results and discussions
The calculated DHS partition functions can either be compared to known theoretical models or exploited to evaluate structural characteristics of the system.The first approach is very challenging, as only approximate partition functions can be obtained analytically.As for the cluster properties, they are usually estimated for finite particle concentrations.Consequently, one needs to separate the influence of inter-cluster interactions carefully.Here, we propose a combined approach.We use our numeric results to find the range of validity of known theoretical models for the partition functions of isolated clusters and for the free-energy density functional.

Chains, rings and branched structures: relative relevance
The analysis of the conformations generated during the GCMC simulation allows us to evaluate the relative fraction of configurations which are in chains, rings and in branched structures for all investigated T. As alluded previously, these fractions are proportional to the partition functions of rings Z ring N , chains, Z chain N and branched structures of the same size, normalized by the cluster partition function Z N .The results of such analysis are shown in Fig. 1.The data shows the complex interplay that determines the thermodynamic stability of the different clusters.On cooling, initially chains become more and more abundant, favoured by the progressive increase in the number of inter-particle bonds.By contrast, upon further cooling they give way to rings, structures which are favoured by one additional bond.At the lowest simulated temperature, chains have essentially disappeared.On decreasing temperature, branched clusters enter into play and become more abundant for progressively larger aggregates.Since branched clusters can be visualised as interacting chains and rings, 53 in the following we focus predominantly on the elementary units, chains and rings.

Chain and ring partition functions
The only known method to calculate analytically the partition function of a DHS cluster with a given topology is to assume, in one way or another, a factorisation: the cluster free energy is reduced to a product of free energies of particle pairs.Thus, the basic ingredient required to analytically evaluate the partition function of any higher order cluster is the partition function of the dimer, eqn (5).There exist several ways to approximate such integral. 68Here, we use the one developed by Morozov and Shliomis in ref. 59, which is obtained through expanding the Hamiltonian in series in the vicinity of the close contact in a head-to-tail orientation and reads (coherently with the definition of Z N in eqn (3)) where v 0 = ps 3 /6 is the particle volume.From here on, T = s 3 /(bm 2 ) stands for the dimensionless temperature.To this day, the expression in eqn ( 15) is the most accurate approximation for the dimer partition function known in literature.
The first term of such expression contains the main energetic contribution and the sum in the brackets comes from the fluctuations of orientations and interparticle distance.It is worth mentioning that equivalent, precision-wise, expressions for a three-particle cluster are not available yet.In Fig. 2, we show the comparison of the simulation data and analytical predictions, eqn (15), for various orders of the expansion.It turns out that, even though we are interested in the low-T limit, fluctuations from the head-to-tail orientation still matter, i.e. even for relatively large values of magnetic interactions, the last terms in the brackets are not negligible.
Dealing with larger clusters becomes significantly more complicated due to the fact that, even though the energetic contribution can be very well approximated by exp(2/T)T 3 /3, the fluctuations are defined not only by temperature, but also by the size and the topology of the cluster.As a consequence, analytical expressions for the partition functions are usually calculated taking into account nearest neighbours or nextnearest neighbour interactions only. 69,70Unfortunately, such an approximation turns out to be insufficient to describe the where C(N) is defined as The latter expression is the result of considering all the distances between each pair of particles.Indeed, in a chain of size N there are N À 1 spheres separated by one diameter, N À 2 particles separated by two diameters, and so on, and only one pair of particles is separated by N À 1 diameters.
Similarly, the expression employed in ref. 51 for a N-particle ring partition function has the form: where The term 1/N 3n+1 captures the difference in entropy between chains and rings arising from (i) the entropic contribution of the chain ends and (ii) the N ways of opening a ring to form a chain; n = 0.588 is the self-avoiding random walk exponent, R ðNþ1Þ=2 stands for the residual of the division and [Á] has the meaning of the integer part of the expression in the brackets.This expression was obtained by calculating all pair interactions in a perfect ring with dipoles oriented tangentially. 71ven though eqn ( 16) provides a nearly quantitative description of the aggregate formation in DHS systems at low and moderate densities, a closer look reveals a sizeable discrepancy, shown in Fig. 3. Here, the results of the simulations differ rather noticeably from theoretical predictions plotted with dashed lines, albeit only quantitatively.The partition functions calculated from eqn (17)  are plotted in Fig. 4 with dashed lines.Analogously to the case of the chain partition function, we observe a small but noticeable difference with simulation data.
The task of improving the analytical expressions for chain and ring partition functions is not straightforward and will be the subject of future work.However, already at this stage we can profit from the knowledge of accurate numeric values for the partition functions and analyse the range of validity for the commonly used Density Functional Theory (see, e.g., ref. 51  and 53) to predict the DHS self-assembly and phase behaviour.

Cluster analysis
The knowledge of the cluster partition function makes it possible to evaluate the cluster size distribution at any density under the hypothesis of an ideal gas of chains and rings.The ideal gas of cluster approach includes the strong association that originates the formation of aggregates, but it neglects mutual interactions between these aggregates.In our case, this corresponds to neglecting all branched structures which can be considered as arising from the interaction between chains and rings.
Indeed, by minimising the free energy of a gas of noninteracting clusters with respect to the cluster size distributions g(N) (chains) and f (N) (rings), (where g(N) and f (N) are defined as the number density of cluster of size N and g(1) is by definition the number density of monomers, e.g.chains of length one) one obtains 51,53,66,67 If Z chain N and Z ring N are known and the number density of branched structures is negligible, then the mass-balance equation provides   (17).At the highest temperature investigated here, T = 0.155, the lack of symbols at large N is due to the chance of finding rings of large sizes, N \ 100, being very small compared to finding branched structures of the same size (see Fig. 1).
the unknown value of g(1) as Strictly speaking, our approach is only valid at infinite dilution, when it is possible to consider the system as a collection of non-interacting chains and rings.Thus, comparing the distributions calculated via Grand-Canonical Monte Carlo simulations to those obtained through full simulations, 32 we can find the DHS density at which the assumption of a system composed solely of non-interacting chains and rings breaks down.It is worth mentioning that there is one more difference between the two methods: single-cluster simulations are performed in vacuum, with no periodic boundary conditions such that the long-range nature of the interactions arising from the presence of other clusters is neglected.By contrast, the full canonical simulations employ periodic boundary conditions and hence require a special treatment for the long-range dipolar interactions, in this case performed through Ewald summations. 72n Fig. 5 we report the chain size distributions g(N)/g(1) calculated from the MC simulations performed in ref. 32 and compare them with the cluster size distributions g(N)/g(1) obtained from the GCMC simulations for several values of the DHS density.For the lowest density the chains are short and, for T 4 0.125, a clear exponential decay is present.In addition, in this temperature range, the fraction of particles in chains grows monotonically with decreasing temperature.At lower temperature the tendency changes: the fraction of particles in chains decreases, and the exponential decay is replaced by a more complex dependence, which gets more resolved upon increasing density, finally becoming a shoulder for N A (10. ..100).Up to r = 0.007, there is no influence of cluster-cluster interactions on the chain distributions, as the canonical (symbols) and single-cluster (lines) results coincide within the noise.However, the agreement worsens as r exceeds 0.028, as seen in Fig. 5(d).
The behaviour of rings is presented in Fig. 6.Here, we see that for low densities our new method provides another advantage: the fraction of rings at T = 0.155 can be accessed.In fact, the probability of observing a ring is so small that they are effectively not seen in the canonical simulations, whereas single-cluster simulations can even quantify it.At low r the fraction of rings has a clearly pronounced maximum that shift towards larger rings on cooling.The growing probability of finding a ring at low T explains the sudden decrease of the probability of finding a chain in the system.As for the chains, inter-cluster interactions become important at r = 0.028, and the cluster GC method fails to describe accurately the fraction of rings, especially at low T.
From the cluster-size distributions we can conclude that the assumption of non-interacting chains and rings starts failing at r B 0.028.This tendency can be seen even more clearly in Fig. 7, where we plot the average cluster size as a function of r for different temperatures.Here, one can see the comparison between canonical Monte Carlo results and theoretical estimates, computed through two different density functional approaches.The average cluster size -an integral characteristics of the systemturns out to be a very good indicator for the influence of the inter-cluster interaction.Indeed, we see that the limitations of the ideal gas of chains and rings manifest themselves at r o 0.001.The behaviour of the solid lines, which show the outcome of the density functional theory and cluster partition functions calculated here, demonstrate that according to the assumption of non-interacting chains and rings, the average chain length grows very rapidly with the particle concentration.At the same time, the usage of the analytically approximated partition functions (dashed lines) seems to better describe the average cluster size up to higher values of r.Such a serendipitous agreement has been previously taken as evidence to set the range of concentrations in which the cluster-cluster interactions can be neglected. 53owever, the much better agreement, compared to the old approach, of the new data with the simulation data at very low concentrations demonstrate without any doubts the superiority  of the former.As a result, the accurate calculation of the partition functions performed in this work reveals strikingly low values of particle concentration below which only chains and rings are present.At densities higher than these values branched (and eventually percolating) clusters become important and must be taken into account to correctly describe the system.This is especially important considering that real ferrofluids rarely contain less than 7-10 volume percent of magnetic nanoparticles, i.e. at least one order of magnitude more than the aforementioned limits.
The quality of the ideal-gas assumption affects more observables than the average cluster size.In Fig. 8 we present the fractions of particles aggregated in rings and chains, using the same notation convention.These plots show that the previously used expressions for partition functions (plotted with dashed lines) were responsible for the low-concentration (Fig. 8a) deviations of the DFT predictions and the results of the canonical Monte Carlos simulations (symbols).The last panel of the figure shows that already at r = 0.007 the majority of the particles in the canonical MC simulations are neither in a chain nor in a ring, demonstrating that the deviations between solid lines and points is to be ascribed to the formation of branched structures.The formation of the latter is a clear consequence of the arising of interactions between chains and rings.This effect becomes more and more important as r increases.It is worth noting that the fractions of particles in chains and rings seem to be more sensitive to the inter-cluster interactions than the cluster size distributions.
In the last part of our discussion we present the ratios of ring fractions to chain fractions as a function of cluster size.This ratio can be found in Fig. 9 for two temperatures.The comparison between the results of the method developed here and of the DFT with approximated partition functions leads to an important observation: even though the two approaches yield different absolute values of the partition functions, as shown in Fig. 3 and 4, they both well reproduce the ratio between the ring and chain partition functions.It is clear that for an infinite ring and an infinite chain the difference between the partition functions should vanish.Besides that, as the accurate calculations of partition functions performed here show, chain and ring free energies behave in a qualitatively similar fashion with growing cluster size.However, the subtle differences in interparticle correlations within these structures have a clear effect on the fractions of various self-assembled clusters.

Conclusions
In the present manuscript we reported the numerical evaluation of the partition function of clusters composed of bonded DHS, specifically focusing on chains and rings.The GC Monte    9 The ratio of ring to chain concentrations as functions of the cluster size for T = 0.125 (top, in black) and T = 0.155 (bottom, in blue).Solid lines are results of the density functional calculations with the cluster partition functions calculated with the method developed here; dashed lines are the predictions of the density functional theory with the cluster partition functions calculated analytically; symbols are obtained from canonical MC simulations for r = 0.007. 51

N
and branched structures of the same size Z branched N , it is possible to evaluate the ring and chain partition functions by multiplying the total partition function Z N by the probability of observing a ring P r or a chain P c .

Fig. 1
Fig. 1 Probability of observing a chain P c (a), a ring P r (b) and a branched structure P b (c) as a function of the cluster size N.

Fig. 2
Fig. 2 Plot of the logarithm of the dimer partition function normalised by the monomer partition function, ln(Z 2 /Z 1 ).Inset: Relative contributions of the various terms to the expansions, expressed in per cent.

Fig. 3
Fig. 3 Plot of the logarithm of the N-particle chain partition function vs. N.The log of the partition function has been divided by N À 1 to highlight the contribution of each bond in the chain.Symbols are simulation results; dashed lines are predictions of eqn (16).

Fig. 4
Fig.4Plot of the logarithm of the N-particle ring partition function normalized by the number of bonds in the ring (N).Symbols are simulation results; dashed lines are predictions of eqn(17).At the highest temperature investigated here, T = 0.155, the lack of symbols at large N is due to the chance of finding rings of large sizes, N \ 100, being very small compared to finding branched structures of the same size (see Fig.1).

Fig. 5
Fig. 5 Chain size distributions at density (a) r = 0.00001 (b) r = 0.0005, (c) r = 0.007, (d) r = 0.028 as a function of the cluster size for different values of the temperature.The lines indicate the results from the singlecluster simulations and the symbols are directly obtained from the canonical MC simulations.

Fig. 6
Fig. 6 Ring size distributions at density (a) r = 0.00001, (b) r = 0.0005, (c) r = 0.007, (d) r = 0.028, as a function of the cluster size for different values of the temperature.Lines are results from single-cluster simulations, symbols are obtained from canonical MC simulations.33

Fig. 7
Fig. 7 Average cluster size of (left) chains and (right) rings as a function of density for systems at (top, in black) T = 0.125 and (bottom, in blue) T = 0.155.Solid lines are results of the density functional calculations with cluster partition functions calculated with the method developed here; dashed lines are the predictions of the density functional theory with cluster partition functions calculated analytically; symbols are directly obtained from canonical MC simulations.

Fig. 8
Fig.8Fraction of particles in chains (in black) and rings (in red) as functions of temperature.(a) r = 0.00001, (b) r = 0.0001, (c) r = 0.0005, (d) r = 0.007.Solid lines are results of density functional calculations with the cluster partition functions calculated with the method developed here; dashed lines are the predictions of the density functional theory with the cluster partition functions calculated analytically; symbols are obtained from canonical MC simulations.51

Fig.
Fig.9The ratio of ring to chain concentrations as functions of the cluster size for T = 0.125 (top, in black) and T = 0.155 (bottom, in blue).Solid lines are results of the density functional calculations with the cluster partition functions calculated with the method developed here; dashed lines are the predictions of the density functional theory with the cluster partition functions calculated analytically; symbols are obtained from canonical MC simulations for r = 0.007.51