Cluster formation in binary fluids with competing short-range and long-range interactions

Abstract The equilibrium phase behaviour of a model binary fluid is investigated through Monte Carlo simulations and by developing a molecular thermodynamic model. Both fluid components interact through a hard core with short-range attractions (SA), but one of the components exhibits an additional long-range repulsion (SA+LR). We find that phase behaviour for this system is controlled by the cross-interaction between the two types of particles as well as their chemical potentials. For a weak cross-interaction, the system displays behaviour that is a composite of the behaviour of the individual components, i.e. the SA component can display bulk vapour/liquid phase separation, while the SALR component can display giant micelle-like clusters for a suitable combination of SA and LR interactions. For a strong cross-interaction, qualitatively different behaviour is observed, with the resulting clusters typically composed of a more equal mixture of SA and SALR particles. Moreover, these mixed clusters can exist even when the SA component by itself would be undersaturated or supercritical, and/or when the SALR component by itself would not form giant clusters. These insights should help to identify the mechanisms for clustering in experimental systems where giant equilibrium clusters are observed. GRAPHICAL ABSTRACT


Introduction
This work investigates the equilibrium phase behaviour of binary fluids where one component exhibits shortrange attractive and long-range repulsive (SALR) interactions. Such fluids are suitable for modelling the formation of giant clusters in solutions of biomolecules [1,2], and other soft matter systems where solutes can become charged in solution (i.e. by proton exchange with the solvent). In these systems, attractive interactions between solutes can arise through a variety of mechanisms, such as hydrophobicity, hydrogen bonding (or specific binding interactions) or a depletion interaction (induced At low concentrations, the behaviour of the system can be similar to that of amphiphilic solutions as long as the short-ranged attractions and long-ranged repulsion forces are nearly balanced. Thus, a critical cluster concentration (CCC) exists above which a cluster fluid phase is present, and at higher densities, a series of modulated phases are thought to occur. The picture is perhaps less clear for strongly attractive and very short-range interactions that lead to irreversible bond formation on applicable timescales. When this occurs, phenomena such as kinetic arrest and gelation obscure any emerging equilibrium properties. These issues can be present in colloidal and protein solutions (i.e. small clusters and kinetic arrest are commonplace).
A variety of single-component SALR systems have been studied using theoretical models as well as simulations. For example, Archer and co-workers [8,9] investigated modulated phases at intermediate densities using mean-field density functional theory and Monte Carlo simulations. More recently, Sweatman and co-workers [10,11] studied low-density SALR phase behaviour using a novel micelle-like thermodynamic model and Monte Carlo simulations. They found micelle-like behaviour as well as a first-order phase transition from a cluster vapour to a condensed cluster phase, driven by depletion of the surrounding SALR vapour, and proposed a simple equation for the average cluster size in these systems.
While the single-component SALR model remains very useful for revealing fundamental aspects of clustering behaviour, many real systems involve more than one solute. This is most obviously the case in biological systems, but it is also true in small molecule solutions in which giant clusters are observed. In each of these cases, all species of the main solute are observed (e.g. all the different ionic forms of an amino acid in aqueous solution), as well as, presumably, low levels of impurity. For this reason, it is important to investigate the clustering behaviour of solute mixtures where at least one of the solute components interacts like an SALR particle.
In this work, we aim to explore the equilibrium phase behaviour of binary mixtures where one of the components exhibits an SALR potential and the other interacts through simple attractive forces. We extend the thermodynamic model developed for the pure SALR fluid [10] to this binary case. A full derivation and explanation for the pure SALR case is provided in that earlier work, and therefore only brief details are provided here.
In the next section, we lay out the basics of the binary fluid model under study. In Section 3, Gibbsensemble Monte Carlo simulations are described that reveal the underlying phase behaviour of this system. For more detailed insights, the thermodynamic model for this binary system is developed in Section 4. The thermodynamic model results are discussed in Section 5, and conclusions are drawn in Section 6.

Binary fluid model
The system of interest in this work comprises a binary mixture of (a) spherical particles interacting through a hard core of diameter d plus a short-range attraction (component 1, by itself a simple SA fluid) and (b) particles identical to (a) except they also display long-range repulsion interactions (component 2, by itself a SALR fluid). For separations r > d, convenient expressions for the direct and cross-interactions are provided by Yukawa potentials, where, with β = 1/(k B T) and i,j = 1, 2. The parameters A ij are all positive and determine the strength of attractive interactions relative to k B T, while A r likewise determines the strength of the repulsive interaction between particles of component 2. Parameters z ij are all positive and determine the inverse decay length of the attractive interactions, while z r (with min(z ij ) > z r > 0) determines the inverse decay length for the repulsive interaction. For convenience, we fix the following parameters: all z ij = 1.00, z r = 0.25 and A r = 0.50, knowing they are adequate for modeling the SALR fluid interactions. Therefore, there are five independent parameters that remain to be defined: the attraction strengths (A ij ) and the ' configurational' contribution to the chemical potential for each component (μ i ).
In the remainder of this work, specific choices are made for A 11 and A 22 , and the resulting fluid phase behaviour is investigated for a range of values of the cross interaction strength A 12 , and the chemical potentials μ 1 , and μ 2 . All parameters are expressed in their reduced form, with an energy scale of β and length scale of d.

Monte Carlo simulations
Gibbs-ensemble Monte Carlo simulations were carried out with 4100 SA particles in a total volume yielding an average density of ρ b1 = 0.1907;, between the coexisting vapour and liquid densities in the absence of component 2, but with different number of SALR particles n 2 = 300, 600, 900 and 1200. The attractive interactions of components 1 and 2 are fixed to A 11 = 0.5 and A 22 = 1.75, with inverse decay lengths of z 11 = z 22 = 1.00, while the repulsive part of the SALR fluid is fixed to A r = 0.50 with z r = 0.25. Four different scenarios were simulated with increasing cross-interaction between components 1 and 2, so that A 12 = 0.00, 0.25, 0.50 and 0.75 while keeping z 12 = 1.00. The cutoff for the interactions is 15 times the hard sphere diameter used.
To enable efficient sampling of both particle and cluster degrees of freedom, a dual Monte Carlo step size was employed with maximum step sizes of 0.10 and 1.50 hard spheres, respectively (chosen randomly), and simple cluster moves were also used. These cluster moves consist of (a) randomly selecting a particle and (b) moving all particles within a specified range of this particle along the same vector. Accepting or rejecting the moves depend on the usual Boltzmann factor, although, to ensure microscopic reversibility, cluster moves are automatically rejected if the reverse move would involve a different group of particles. Figure 1 shows snapshots from these Gibbs-ensemble simulations at equilibrium. Only four types of phase are observed, with their existence dependent on the component densities and the strength of the crossinteraction; uniform vapour, uniform liquid, vapour with clusters and liquid with clusters. When A 12 = 0, we find the following sequence of equilibrium behaviour as the concentration of SALR particles (component 2) increases: (i) vapour/liquid coexistence, (ii) vapour with clusters coexisting with uniform liquid, (iii) vapour with clusters coexisting with liquid with clusters. This shows the SALR clusters preferentially partition into the vapour phase when A 12 = 0. When A 12 = 0.25, the sequence is similar except the SALR clusters appear to preferentially partition into the liquid phase, i.e. the liquid phase with clusters has become more stable. The same sequence is observed for A 12 = 0.50, except in this case, the CCC has increased, i.e. the uniform liquid has become more stable. When A 12 = 0.75, we see qualitatively different behaviour in that the liquid no longer exhibits giant clusters at all, and giant SALR clusters apparent in the vapour phase are mixed and swollen with component 1 (the SA component).
A tentative generic phase diagram for low SALR particle concentrations is sketched in Figure 2 based on these Gibbs-ensemble results. The usual vapour-liquid coexistence occurs when SALR particles are absent if (as is the case here) the SA system by itself is subcritical. Red lines denote the critical cluster chemical potential (CCμ) of SALR particles in the vapour and liquid phases, assuming that the SA and LR interactions are balanced such that they lead to SALR clusters in the absence of component 1. The curved arrows indicate that the slopes of these critical cluster lines will depend on A 12 . Likewise, the horizontal arrow indicates the position of the liquid CCμ will vary with A 12 . For low and high values of A 12 , the Gibbs-ensemble simulation results suggest it will initiate to the right of the end of the vapour CCμ. But for intermediate values of A 12 , it can initiate to the left of the vapour CCμ. For sufficiently high values of A 12 , it appears from the simulations that clusters are no longer stable in the liquid phase. Of course, being a sketch, we have used straight lines, but we fully expect all the lines in this diagram to be curved in practice.
While these simulations provide an indication of the equilibrium behaviour for the mixtures studied, a detailed exploration of the phase diagram quickly becomes time consuming. For every choice of the set of model parameters, many simulations are required to adequately explore the phase space in the cluster region and provide relevant averages for cluster sizes and compositions. Given that each simulation can take many hours, this situation is not ideal. Particularly as we wish to be able to make fast predictions of phase behaviour to assist in the design of nanoparticles synthesised via this giant SALR cluster route. Just as importantly, we cannot guarantee that these simulations are truly at equilibrium. This would require Monte Carlo moves that swapped clusters between simulation boxes.
We therefore developed a thermodynamic model, described next, for fast prediction of phase behaviour and to confirm the tentative phase behaviour shown in Figure 2.

Thermodynamic model
The thermodynamic model of Sweatman et al. [10], which was developed to investigate the cluster fluid phase of the one-component SALR fluid, is now extended to this binary mixture. The resulting model can be viewed as a kind of density functional micelle theory. An advantage of this particular approach is its basis in statistical thermodynamics and density functional theory, i.e. it relates the properties of constituent particles (the interaction model or Hamiltonian) to resulting phase behaviour, rather than taking emergent properties, such as the surface tension, as input to the model, as is more common in theories of clustering and micellization. The model is designed to treat each of the four phases seen in the simulations in Section 3.
As the one-component SALR thermodynamic model was described in detail in an earlier work [10], here we present only a summary of the binary model used in this work, stressing the key modifications for its extension to binary systems. The approach is based on the ' capillary model' approximation and treats the cluster fluid phase as a disordered collection of spherical liquid-like droplets of volume V c = π d 3 c /6 each, dispersed within a background vapour. Each cluster of particles (or droplet), with diameter d c = 2R c , is taken to be identical with uniform body densities ρ L1 = n L1 /V c and ρ L2 = n L2 /V c for the two components, respectively, where n L1 and n L2 denote the number of particles of type 1 and 2, respectively, within each cluster. The cluster density is ρ c = n c /V, where n c is the number of clusters within a system of total volume V, and the corresponding average background fluid densities are ρ v1 = n v1 /V and ρ v2 = n v2 /V, where n v1 and n v2 are the number of particles of each type in the background fluid.
Although, in general, clusters are polydisperse in their size and shape, they are statistically identical, and so in this model, they are chosen to be spherical and monodisperse. Therefore, the profile of a cluster centred at R R R is, where is the Heaviside step-function. Thus, the ' capillary' model defines a discontinuous interface between the clusters and the background fluid. The volume fraction ψ of the clusters is, where ρ b1 and ρ b2 are the total, or bulk, density of particles of type 1 and 2, respectively. Unlike for the original derivation, where the model was presented using the Canonical ensemble (i.e. in terms of the Helmholtz free energy), for the binary mixture case, here we rewrite it in terms of the Grand-potential density, ω c , where μ i is the chemical potential of component i. Following the original approach by Sweatman et al. [10], the free energy contributions are split into two contributions: energetic contributions coming from pair interactions, and entropic contributions obtained by treating particles and clusters as hard spheres. The Helmholtz free energy (f c ) can also be partitioned into 'self ' and 'mixture' contributions, In this equation, f self describes the free energy density of particles within clusters, while f mix contains contributions to the free energy density of the mixture from both the background fluid and the clusters. Simple expressions for f self and f mix are now developed in terms of the basic parameters describing particle interactions. At all times, we only consider the configurational contribution to the free energy, as the contribution from momentum (involving the thermal de Broglie wavelength) is unimportant for equilibrium behaviour.

Self-free energy
The self-free energy of the clusters is composed of two terms, In Equation (10), U self is the interaction energy between particles in a cluster, i.e. its self-energy, and S self is the entropy of the particles within a cluster. To estimate U self , we assume configurational correlations between particles within a cluster are similar to those within a bulk liquid. So the interaction energy of particles within a cluster is given by, where g Lij (r) is the radial distribution function (RDF) between particles i and j in a bulk liquid, and P dc (r) is related to the form factor of the clusters and represents the geometric convolution of two cluster distributions. The liquid RDFs are approximated by those of a hard sphere fluid, HS (r; ρ Ltot ; d), (12) where the Percus-Yevick (PY) theory [12] for hard spheres at the combined liquid density ρ Ltot = ρ L1 + ρ L2 is used. The self-entropy is approximated as, is the entropy density of a uniform hard sphere fluid mixture of components with densities ρ L1 and ρ L2 , respectively, and A c is the cluster interfacial area.
Unlike in the original model derivation [10], the model now also takes into account the entropic contribution of the primary particles at the interface between the background fluid and the clusters, represented by the right-most term in Equation (13). Assuming a planar step-like density profile at the interface, the interfacial entropy density is The excess free energy density for the gas-liquid interface f ex HS A is calculated using the Rosenfeld et al. [13] fundamental measure functional (FMF) for hard spheres, and the uniform background fluid, f ex HS g , and cluster liquid, f ex HS L , hard sphere excess free energy densities are calculated from the same in the limit of a bulk fluid mixture.
In Equation (15), there are only three independent parameters, η 2 (x), η 3 (x) and η v2 v2 v2 (x), which are calculated in the planar limit through, where the independent parameters are calculated using the following weight functions of the Rosenfeld functional, The remaining scalar weight functions are proportional to w (2) i , and the remaining vector weight function is collinear with w (v2 v2 v2) i as indicated by Rosenfeld [13], Finally, the local densities for each component, ρ j (x), on either side of the interface take on the values inside the cluster for x ≤ 0, and in the background fluid for x > 0.

Mixing free energy
We deal with the contribution of the energy density to the mixture free energy first, which is given exactly by the usual energy equation [14], where g vij (r) is the RDF between particles i and j in the background fluid, g ic (r) is the RDF between component i in the background fluid phase and cluster centres, and g cc (r) is the cluster center-cluster center RDF. Also, φ ic (r) is the effective pair-potential between component i in the background fluid phase and cluster centers, and φ cc (r) is the effective pair-potential between cluster centers. Convenient approximations for all these functioare made as follows [10] g vij (r) = g HS (r; ρ gtot , d), g ic (r) = (r − (r c + 2d))exp(−βφ ic (r)), g cc (r) = (r − (2r c + 2d))exp(−βφ cc (r)), where the combined vapour density is ρ gtot = (ρ v1 + ρ v2 )/(1 − ψ). Here it is assumed that no particle can approach a cluster (that it is not part of) closer than 2d, since then it would be considered part of that cluster. These cluster-cluster and cluster-vapour g(r) approximations are adequate for this work. In any case, more accurate approximations are not obviously available, since integral equation approaches to this problem of giant SALR clustering typically fail in the low-density region of the SALR phase diagram in which we are interested. The effective interactions are computed in k-space as, The entropy of the coarse-grained system, due to these effective interactions, is approximated in terms of a hard sphere mixture with significant negative non-additivity (see Sweatman et al. [10] for a discussion on this). As there is no satisfactory existing expression for the entropy of such a system, the background fluid and cluster excess entropies are decoupled, where d eff c is the effective cluster diameter calculated via the Barker-Henderson prescription [14], which also defines an effective packing fraction of clusters, Again, we use Rosenfeld's FMF [13], which reduces to the PY compressibility result [14] for a bulk fluid, to calculate the excess hard-sphere free energy density. The ideal gas contribution in Equation (29), omitting terms involving the thermal de Broglie wavelength of the primary particles, is simply, The term S com in Equation (29) is a centre-of-mass correction, which ensures the coarse-grained clusters have a precisely defined position (again, see Sweatman et al. [10] for a discussion on this point) and can be seen as contributing to the entropy of formation for the clusters.
For this contribution, we used an appropriately modified version of our previously derived expression [10] for the pure SALR system, In Equation (33), choosing α = 1 gives the same expression derived originally for the pure SALR fluid in earlier work. For the binary system in this work, we found for low cluster body densities (ρ Ltot ≈ 0.2) that the formation of unrealistically small clusters with less than four particles was sometimes favoured due to an under estimation of the cluster formation entropy (i.e. the S com contribution) which cannot adequately compensate the ideal-gas cluster contribution f id (ρ c ), which favours the formation of clusters. This problem arises because of an error in the derivation of this term in [10]. Indeed, it is in principle impossible to specify the formation free energy of a cluster precisely without first specifying how a cluster is defined, i.e. specifying the rules that govern when a particle is considered part of a cluster. As we do not specify those rules here, we cannot in principle exactly formulate an expression for S com . Viewed another way, the final term on the right hand side of (33) is arbitrary, and only determined precisely when a rule governing the definition of a cluster is imposed. As we do not do this here, we must regard this final term as an adjustable parameter. Thus, we adjust α to ensure that the smallest cluster is comprised of at least four particles. The parameter that satisfies this condition, and therefore the one used in the remainder of this work is α = 0.05198; (see Figure 3). This completely defines the thermodynamic model for the cluster fluid phases of the binary mixture. For a uniform fluid without clusters, the free energy is expressed simply as, where ρ btot = ρ b1 + ρ b2 .

Cluster solid model
For sufficiently high SALR chemical potentials, a cluster solid phase is expected [10]. A thermodynamic model of the cluster solid phase was not included in the earlier one-component model; instead a criterion based on the effective packing fraction of clusters, ψ eff , was used to locate this phase boundary [10]. In the present work, we instead prefer to maintain a thermodynamicapproach, and thus develop a thermodynamic model of the cluster solid phase. First we introduce the relative packing fraction, z, which is related to the effective packing fraction and is calculated as follows, with 0 < z < 1. When considering the crystal packing of solid spheres, there are two main staking arrangements that are favoured: face-centred-cubic (fcc), and hexagonalclosed-packed (hcp), since they both provide the most efficient packing (occupying about 74% of the space). These two structures have been the focus of many studies, with Sion-Chuon and Huse [15] reporting that the fcc stacking has the highest entropy. We therefore proceed to model the cluster solid phase of our binary system as an fcc solid.
In order to account for the cluster solid phase, the thermodynamic model requires two modifications. First, the cluster-cluster RDF in Equation (26), g cc (r), which is suitable for a cluster fluid, must be replaced with a suitable RDF for the fcc cluster solid. For computational efficiency, we use an empirical fit to the hard sphere fcc crystal RDF developed by Choi and Ree [16] valid for separations up to r/d < 3.3 provided 0.700 < z < 0.985. For larger separations, we take g cc = 1.
Second, the entropy density of clusters, which for the cluster fluid phase is given by the term f ex HS (ρ c , d eff c ) in Equation (29), must be modified to account for the fcc crystal phase. We use an empirical equation proposed by Speedy [17], with the fcc parameters being: a = 0.5921, b = 0.7072, c = 0.601 and S o = −2.531. This expression provides the excess entropy per hard sphere at a packing fraction z and is valid for 0.643 < z < 0.990. To employ this equation of state for our fcc phase of clusters, we use the relative packing fraction (Equation 35) and multiply by ρ c to obtain a reduced excess free energy density. As our main interest here is to locate the cluster fluid to cluster solid phase transition, we are not interested in cluster solid phases with high packing fractions and so impose the constraint z < 0.9.

Solving the model
Minimisation of the grand-potential density with respect to the primary variables (i.e. the vapour densities ρ vi , the cluster body densities ρ Li , the cluster size d c and the cluster volume fraction ψ at fixed chemical potentials (μ i ) with a defined set of interaction strengths A ij is performed via a local down-hill search algorithm. Phase space is searched on a square grid of chemical potentials. For each grid point, we compare the grand free energy densities calculated for uniform fluid, cluster fluid and cluster solid, taking the lowest free energy as the equilibrium phase. We impose the constraint M > 4 consistent with Figure 3.
In what follows we present our results in terms of reduced fugacities rather than chemical potentials, where the reduced fugacity f * i can be readily obtained from the configurational part of the chemical potential as f * i = [πd 3 hs /6] e βμ i . This reduced fugacity equals the packing fraction of hard spheres in the dilute limit.

Thermodynamic model validation
While the thermodynamic model was developed for binary mixtures, it should also reproduce the phase behaviour of the limiting cases: i.e. the pure SA fluid and the pure SALR fluid. The model is forced to mimic the respective pure fluid by setting the chemical potential of the other component to an arbitrarily low number.
The SA fluid phase diagram is presented in Figure 4. When the attractive parameter is very low (equivalent to high temperatures), the fluid shows supercritical behaviour. As A 11 increases (equivalent to decreasing the temperature in the system), the fluid reaches a critical point at about A 11 = 0.33 eventually showing the typical first-order gas-liquid phase transition with fugacity, as expected for this system (Figure 4(b)).
Results for the predicted phase diagram of a pure SALR fluid with A r = 0.50 and z r = 0.50, identical to the parameters used originally by Sweatman et al. [10], are presented in Figure 5, in which the cluster-fluid phase is shown in red and the cluster solid phase in blue. These results were obtained by sweeping A 22 in 0.025 intervals, and μ 22 in 4 · 10 −3 intervals, and they agree well with those presented by Sweatman et al. [10]. Figure 5(b) shows results specifically for A 22 = 1.85 studied by changing the chemical potential in small intervals of 5 · 10 −4 . As expected, increasing fugacity leads to the appearance of a cluster-fluid phase. Further   increasing fugacity results in increased volume fractions until, at high fugacities (i.e. high SALR concentrations), there is evidence of a discontinuous cluster-fluid to cluster solid phase transition.
We then proceed to map the phase diagram for an SALR fluid with parameters identical to those used in the Monte Carlo simulations shown in Section 3 and which form the focus of this present work, i.e. with a longer-ranged repulsive interaction, z r = 0.25. Figure 6(a) shows results obtained by incrementing A 22 in intervals of 0.02 and μ 22 in intervals of 5 · 10 −2 . Although the cluster-fluid phase now exists over a wider range of fugacities, these correspond to a smaller range of bulk densities. As in the previous SALR system, at high fugacities, there is a discontinuous phase change from cluster-fluid to cluster solid, but now there is also a first-order phase transition between the uniform and cluster-fluid phases.

Binary mixtures
From here on the discussion will focus on SA+SALR binary mixtures. In particular, two different cases will be discussed: supercritical SA fluid (represented in an ideal sense using hard spheres, A 11 = 0.00), and a subcritical SA fluid (A 11 = 0.50). A summary of the relevant parameters in both cases is presented in Table 1. We then investigate the effect of changing the strength of the cross-interaction strength (A 12 ) in each case. Unlike the pure fluid case, model solutions for binary systems require more computational effort, as each specific grid point on the f * 1 , f * 2 phase diagram involves a complex multi-variable global minimisation, which is susceptible to the initial conditions chosen. For this reason, our final phase diagrams for the binary mixtures are found by comparing four different sets of solutions obtained by sweeping the phase space up/down on each chemical potential. Phase diagrams are conveniently viewed in terms of the cluster volume fraction, φ.

Mixture I
We initially consider this mixture (parameters in Table 1) for the case where A 12 = 0.00, modelling an SALR fluid containing inert impurities. The phase diagram obtained (Figure 7(a)) shows that for very low SA fugacities the uniform to cluster-fluid phase transition approaches the pure SALR fluid case. However, despite the absence of a cross-interaction, as the SA concentration increases (i.e. increased fugacity) Figure 7(a) shows higher SALR fugacities are needed to achieve the cluster-fluid phase, i.e. the presence of an inert impurity delays clustering, but with no significant effect on the cluster size or the cluster volume fraction at the transition to a cluster solid (Figure 7(b)). Figure 7(c) shows the density and composition of background fluid and clusters for a fixed SA fugacity f * 1 = 0.29 (i.e. μ 11 = −0.60). The background fluid is formed mainly of SA particles at low concentration, while clusters are almost pure SALR with liquid-like density. As the SALR fugacity is increased, we see first-order phase transitions from the uniform fluid to the cluster fluid, and finally to the cluster solid.
The model shows some uncertainty in the cluster solid region, where it often alternates between two solutions, one depicting slightly smaller (by 4%), denser (by 5%) clusters with lower volume fractions (by 8%). This is likely an artefact of our model for the cluster solid phase, which assumes monodisperse clusters, and/or due to the solution method which is not guaranteed to find the global minimum of the free energy.

Mixture I: increasing A 12
Mixture I as presented so far may be relevant for systems with inert impurities. However, real mixtures will always have some degree of cross-interaction (A 12 ), the nature of which plays a key role in the phase behaviour of the system.
As the cross-interaction strength is increased, the behaviour remains similar to that presented in Figure 7, except for the composition of the background fluid, which presents a higher amount of SALR particles at lower overall density. Figure 8 shows the result of increasing the crossinteraction to A 12 = 1.00. Qualitatively different be- haviour is now observed. The transition to the cluster fluid phase is no longer delayed by the SA fluid, rather it is promoted by it: notice the uniform to clusterfluid transition line moves towards the left (i.e. towards smaller SALR fugacities) as the SA fugacity increases (Figure 8(a)). Both the background phase and clusters have a significant proportion of SA and SALR particles, due to the increased SA-SALR inter-particle affinity. In turn, this causes an increase in the cluster size relative to weaker cross-interactions, which then leads to a reduction in the density of clusters.

Mixture II
In this case, the SA fluid is subcritical, thus allowing us to explore the effect of a second component with strong attractive interactions when mixed with the SALR fluid of interest. As before, we first study the case without cross-interactions A 12 = 0.00. Comparing Figures 7  and 9, although we see the same sequence of transitions, from uniform fluid to cluster-fluid to cluster solid, there is a noticeable change in phase behaviour, with respect to Mixture I, when the background fluid condenses.
For low SA fugacities, the uniform phase being a lowdensity vapour (Figure 9(c)), the system behaves almost as a pure SALR fluid, as expected, except with some SA particles in the background fluid. However, at high SA fluid fugacities, when the background phase is liquid-like ( Figure 9(d)), the SALR fugacity required for clustering increases almost linearly with increasing SA fugacity (and hence with SA density). In this regime, clusters are almost pure SALR, while the background fluid is nearly pure SA.

Mixture II: increasing A 12
Slightly increasing the cross-interaction (A 12 < 0.5) shows no significant change in the shape of the phase diagram, but with the noticeable difference that the critical cluster fugacity (CCf ) at the liquid-like and vapour-like homogeneous phases no longer intersect. This is because the CCf in the liquid-like background phase, though still almost linearly dependent on μ 1 , has been pushed back to higher SALR fugacities, indicating in this case, the cluster fluid is less stable when the background fluid is liquid-like.
As the cross-interaction increases, i.e. to A 12 = 1.00 (Figure 10(a)), the transition from vapour to liquid  Table 1  occurs at smaller SA fluid fugacities and for sufficiently low f * 1 , the system behaves almost like a pure SALR fluid. But as the SA fugacity increases, clustering occurs at lower SALR fugacities, just as for the case with a supercritical SA component (Mixture I). These clusters have diameter similar to the earlier cases, with d c < 10.
But, in the region of the phase diagram where the background fluid condenses, where the fugacity of the SA fluid is sufficiently high, the model generates some unexpected results. Under these conditions, the model predicts the formation of giant clusters (d c > 30) with a low body density, not very different to the background phase. In other words, a very weakly modulated phase is formed. It is not clear if this is the correct behaviour or an artefact of the thermodynamic model. This can be further investigated by the use of a mean-field non-local DFT, but this is beyond the scope of the present work.

Conclusions
We developed a thermodynamic model to study the equilibrium phase behaviour of a model SA+SALR binary fluid. Our multiscale DFT approach generates results consistent with those obtained from Gibbs-ensemble MC simulations for the same model mixture.
The DFT results support the MC simulation results in qualitative terms, and therefore the phase diagram sketched in Figure 2 appears to show the correct trends, even for the supercritical case where there is no bulk vapour-liquid phase transition. In particular, we find the cross-interaction largely determines phase behaviour. Its strength has a direct effect on the direction of the uniform fluid-cluster fluid phase boundary, which for low cross-interaction strengths moves towards higher SALR fugacities as the SA fugacity increases, while for high cross-interaction strengths moves towards lower SALR fugacities as the SA fugacity is increased. In other words, for low cross-interaction strengths, the presence of the SA component tends to delay cluster formation, while the opposite is true for strong cross-interactions.
Moreover, giant clusters can form readily, even if they would not form for the SALR component alone, provided the SA component concentration is sufficiently high and the cross-interaction is sufficiently strong. The SA Figure 10. Mixture II with parameters given in Table 1  composition within the clusters is significantly increased with increasing cross-interaction strength, with clusters swollen by the absorption of SA particles.
A useful way to consider the effect of the SA component is in terms of an additional effective interaction between SALR particles. If the cross-interaction is only weakly attractive, then the SA component induces a weak effective repulsion between SALR particles. But if the cross-interaction is strongly attractive, then SA particles induce an effective attraction between SALR particles. In each case, the strength of the effect is sensitive to the concentration (or chemical potential) of SA particles. We expect it might be possible to formalise this observation by integrating-out the SA component to generate an effective one-component SALR system.
This observation is very important for understanding the appearance of these kinds of giant cluster in experiments. Even though charged solutes might, by themselves, not form clusters under experimental conditions, when in the presence of an uncharged solute, which could be a different species of the same solute, clusters can still appear if the cross-interaction between them is strong enough. Essentially, it can happen that the uncharged species tends to 'glue' the charged species together. A giant cluster will still form if the screened coulomb repulsion of the charged species is sufficiently strong, such that formation of the bulk liquid phase is frustrated.
To our knowledge, this is the first theoretical approach capable of investigating both a disordered cluster phase (the cluster fluid) and a modulated phase (the cluster solid). It can also investigate the effect of the addition of a simple fluid component (i.e. the SA particles) on the clustering behaviour of an SALR system. By enabling the consistent treatment of a cluster solid phase on the same footing as the cluster fluid, we observe a first-order cluster fluid to solid transition for sufficiently high particle concentrations (or chemical potentials). In some cases, the model predicts the emergence of a very weakly modulated phase with very long wavelength fluctuations. It is not clear at this stage whether this is an artefact, and therefore a failure, of the model, or whether this behaviour is correctly predicted. Investigation of this feature can be tackled with a standard non-local DFT treatment.
More generally, the model opens the door to investigating giant cluster formation in experiments with multiple solute components, exemplified by soft matter and biological systems.

Disclosure statement
No potential conflict of interest was reported by the authors.