Michela
Ronti
*a,
Lorenzo
Rovigatti
bc,
José M.
Tavares
de,
Alexey O.
Ivanov
fg,
Sofia S.
Kantorovich
af and
Francesco
Sciortino
c
aUniversity of Vienna, Sensengasse 8, 1090 Vienna, Austria. E-mail: michela.ronti@univie.ac.at
bCNR-ISC, Uos Sapienza, Piazzale A. Moro 2, 00185 Roma, Italy
cUniversity of Rome La Sapienza, Piazzale Aldo Moro 2, I-00185, Rome, Italy
dCentro de Física Teórica e Computacional da Universidade de Lisboa, Faculdade de Ciências, Campo Grande, 1749-016 Lisbon, Portugal
eInstituto Superior de Engenharia de Lisboa-ISEL, Rua Conselheiro Emídio Navarro 1, 1950-062 Lisbon, Portugal
fUral Federal University, Lenin av. 51, 620000 Ekaterinburg, Russian Federation
gM.N. Mikheev Institute of Metal Physics UB RAS, S. Kovalevskaya Str. 18, Ekaterinburg, Russian Federation
First published on 4th October 2017
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 offered 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 offer 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.
Nowadays it is possible to functionalise the nanoparticle surfaces,13,14 tune their shape15–21 and control their internal crystalline structure.22–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.25
The main interaction between the dipoles is a long-range anisotropic dipolar interaction. This interaction determines a preferred head-to-tail orientation of dipole moments. Thus, when thermal fluctuations are weaker than the magnetic attraction, magnetic nanocolloids exhibit directional self-assembly in linear flexible chains.26–31
However, the self-assembly scenario is not limited to linear chains: magnetic particles can also form ring structures32,33 and networks.34,35 Experimentally, all these clusters have been observed in solution using cryogenic electron microscopy36 and atomic force microscopy of assemblies at cross-linkable oil–water interfaces.37
Various simulation studies described self-assembly in magnetic nanocolloids by using simple models of dipolar hard- or soft-spheres.34,38,39 Despite the apparent simplicity of these models, the low-temperature portion of the phase diagram of dipolar hard and soft spheres remains obscure, with many different hypothesis having been put forward.40–44 Concerning the low-density portion of the phase diagram, a gas–liquid critical point has been found numerically for sequences of models that asymptotically tend towards the dipolar hard sphere (DHS).41–44 For instance, if the dipolar interaction is complemented by an attractive dispersion interaction a phase separation may occur.45–49 However, recent simulation studies32 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.50
Recently, we showed that there is a structural chain-to-ring transition in a highly diluted gas of DHS as temperature decreases51 and extended this study to higher concentrations where various branched structures were also considered.52,53 These investigations provided a possible solution to a long-lasting debate about the non-monotonic temperature dependence of the initial magnetic susceptibility in suspensions of magnetic nanoparticles.54,55 These 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.
(1) |
(2) |
We define the cluster partition function ZN of a cluster made of N indistinguishable DHS particles as
(3) |
Analytically, one can easily calculate the expressions for
(4) |
(5) |
The ratio between the two has the form:
(6) |
The evaluation of the cluster partition function (or cluster free energy) requires an average of the system Hamiltonian over all particle positions ri and dipolar orientations Ω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βμ, with μ 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 of observing a cluster of size N is proportional to
(7) |
(8) |
(9) |
(10) |
Our operative definition of a bond, initially proposed in ref. 32, employs a simple distance–energy criterion: two particles i and j are bonded if the interparticle distance rij < 1.3σ and if the interparticle potential energy Vdd(rij) < 0. The 1.3σ 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. It is worth mentioning that there exist other criteria to define a bond.59–61 However, we do not expect the results to be sensitive to the precise definition of a bond at the values of T and ρ 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 ±0.05σ and simultaneously randomly rotated by at most ±0.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 σ and width 0.3σ (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-ρ region in equilibrium and to sample the conformations of the self-assembled structures. This move – originally introduced by Lal62 and independently suggested by MacDonald et al.63 and by Mandras and Sokal64 – 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
(11) |
(12) |
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 non-interacting clusters with respect to the cluster size distribution n(N), one obtains66,67
(13) |
(14) |
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 ZN of a cluster of size N can be written as sum of the partition functions of rings ZringN, chains, ZchainN and branched structures of the same size ZbranchedN, it is possible to evaluate the ring and chain partition functions by multiplying the total partition function ZN by the probability of observing a ring Pr or a chain Pc.
Fig. 1 Probability of observing a chain Pc (a), a ring Pr (b) and a branched structure Pb (c) as a function of the cluster size N. |
(15) |
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)T3/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 next-nearest neighbour interactions only.69,70 Unfortunately, such an approximation turns out to be insufficient to describe the DHS self-assembly in the low-T regime. Consequently, the partition function of an N-particle chain that was previously developed in ref. 51 and 53 to analyse the structural transitions in low-density DHS systems takes the form
(16) |
Similarly, the expression employed in ref. 51 for a N-particle ring partition function has the form:
(17) |
Even 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.
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 Plot 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). |
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.
Indeed, by minimising the free energy of a gas of non-interacting 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 obtains51,53,66,67
(18) |
If ZchainN and ZringN are known and the number density of branched structures is negligible, then the mass-balance equation provides the unknown value of g(1) as
(19) |
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.72
In 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 > 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 ∈ (10…100). Up to ρ = 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 ρ 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 ρ 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 ρ = 0.028, and the cluster GC method fails to describe accurately the fraction of rings, especially at low T.
Fig. 6 Ring size distributions at density (a) ρ = 0.00001, (b) ρ = 0.0005, (c) ρ = 0.007, (d) ρ = 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 |
From the cluster-size distributions we can conclude that the assumption of non-interacting chains and rings starts failing at ρ ∼ 0.028. This tendency can be seen even more clearly in Fig. 7, where we plot the average cluster size as a function of ρ 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 system – turns 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 ρ < 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 ρ. 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.53 However, 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 ρ = 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 ρ 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.
Fig. 8 Fraction of particles in chains (in black) and rings (in red) as functions of temperature. (a) ρ = 0.00001, (b) ρ = 0.0001, (c) ρ = 0.0005, (d) ρ = 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 |
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.
Fig. 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 ρ = 0.007.51 |
Interestingly, in the region where simulations done with the single-cluster (which are performed in the void and hence without periodic boundary conditions) and canonical MC (which take into account the long-range nature of the dipolar interactions through Ewald summations) methods can be compared, we observe a nearly-perfect agreement between the two. This reinforces the theoretical approaches that attempt to model the thermodynamic of DHS in term of topological transition73,74 of weakly interacting (through the formation of branching points) elementary clusters (e.g. chains and rings).
We have also evaluated the bulk behaviour of the DHS system under the assumption that branched structures can be neglected, such that the system can be modelled as an ideal gas of non-interacting chains and rings. The result of this evaluation has been compared with bulk simulations of DHS, confirming that at very low densities and temperatures, the system can indeed be modelled as an ideal gas of rings and chains. For larger densities, the effect of branching becomes relevant.
The availability of precise evaluations of the chain, ring and branched partition functions calls for a theoretical effort in the direction of accurately quantifying their N and T dependence. By leveraging the same approach, we plan to examine the GCMC data to extract the partition function of the various branched structures observed previously in systems of DHS.52,53 Developing such understanding will make it possible to map the (T-dependent) DHS interactions into simple models of associating fluids, which are in principle amenable of analytic solution in mean field. This hopefully will provide an accurate modelling of DHS in an extended T and ρ window.
This journal is © The Royal Society of Chemistry 2017 |