Phase behaviour of inverse patchy colloids: effect of the model parameters

The phase behaviour of inverse patchy colloid systems composed of spherical particles with two oppositely charged patches at the poles is investigated by simulation-based thermodynamic integration schemes. The interaction between the particles is derived via a coarse-grained model characterized by three system parameters: the charge imbalance between the bare colloid and the patches, the patch surface extension and the particle interaction range. Starting from a set of parameters for which a stacking of parallel layers is thermodynamically stable, the effect of each of these three parameters on the phase diagram is studied. Our results show that the region of stability of the layered solid phase can be expanded by increasing the charge imbalance and/or by reducing the interaction range. A larger patch size, on the other hand, stabilizes the layered structure with respect to the competing face centered cubic solid at high pressures but destabilizes it with respect to the fluid phase at low pressures. The location of the liquid-vapour critical point in the temperature versus density plane is also investigated: while the charge imbalance and the patch size affect mainly the critical density, a change of the interaction range has a substantial impact also on the critical temperature.


Introduction
During the last years the study of the collective behavior of colloidal particles with anisotropic shapes or interactions has attracted a considerable amount of attention. Advances in colloidal science have made the synthesis of such colloids possible and they nowadays guarantee an exquisite control over the shape of the particles and their interactions [1]. It is thus important to investigate how these particles selfassemble into high ordered structures so that they can be used as building blocks for new materials with desired properties [2][3][4]. Already numerous simulation and theoretical studies have been carried out in that direction, revealing that anisotropic particles can indeed exhibit a complex and unusual phase behaviour.
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In this contribution we will focus on the phase behaviour of the so-called inverse patchy colloids (IPCs), i.e. particles with differently charged surface regions: areas with like charges repel each other whereas regions of opposite charge are mutually attractive [5]. IPCs can be considered as simplified models representing complex systems with heterogeneously charged surfaces, such as proteins, virus capsids [6], spotted vesicles [7], etc. As a consequence of their competing attractive and repulsive interactions, these particles exhibit quite a rich behaviour [8][9][10]. Under planar confinement IPCs assemble into quasi two-dimensional structures where the particle axes are parallel to the plane and the bare equatorial regions are exposed on both sides of the monolayer, preventing thereby the growth of the assembled clusters in the direction perpendicular to the planes [8,9]. This tendency to form planar (laminar) structures has also been observed in three dimensional systems for a particular set of parameters (i.e. moderately overcharged colloids with moderate patch sizes under low screening conditions) [10]. For this particular case, the equilibrium phase diagram was calculated, revealing that the laminar structure is stable over a quite broad range of temperatures and pressures, competing at high pressures with a close-packed face centered cubic (FCC) solid that is orientationally ordered at low temperatures but becomes a plastic crystal (PC) above a given temperature. Furthermore, this model exhibits a stable liquid-vapour phase separation. Laminar phases, similar to the ones found in IPC systems, have also been observed for other simple models, such as Janus particles [11], trivalent patchy particles [12], particles with multiple patches [13] and charged colloids with an inverted dipolar field in the presence of an external field [14]. The identification of the features that stabilize the laminar structure could potentially be useful to gain insight into the behaviour of real systems that grow into similar structures, e.g. certain proteins (e.g. apoferritin) [13] and alloy nanoparticles (e.g. CdTe [15], PbS [16] and PbSe [17]).
In this work we investigate how the different IPC model parameters, namely the charge imbalance between the patches and the bare colloid, the patch size and the particle interaction range, affect the phase behaviour and the stability of the laminar phase. Only moderate changes of the parameters used in [10] will be considered here, so that it can essentially be taken for granted that no other ordered phases besides the laminar structure and the FCC (ordered and plastic) crystals are stabilized. We consider as moderate changes those large enough to lead to an appreciable modification of the phase diagram but small enough to avoid the appearance of new ordered structures. The modifications of the parameters (namely, colloid charge, patch size and interaction range) are typically between 10% and 20% of the reference value. Under these conditions, the phase diagram for the new systems (defined by the modified sets of parameters) can be efficiently computed by the so-called Hamiltonian Gibbs-Duhem (HGD) integration technique [18,19], using as starting points the coexistence lines calculated in [10]. The usefulness of the HGD method has already been proven for a large variety of systems [20][21][22].
Obviously, large changes on the model parameters of the reference IPC system [10] can lead to a complete destabilization of the laminar phase and to the appearance of new ordered structures different from those found in [10] (see, for instance, [23]). The study of those new assembly scenarios is certainly of high interest and will be the subject of future research. However in the present work the investigations focus on moderate changes in the parameters that allow us to study their impact on the stability of the laminar phase. Furthermore, also the effect of the model parameters on the critical point is of interest, as the interplay between the attractive and repulsive terms in the potential might lead to changes in the localization of the liquid-vapour phase separation with respect to the fluidsolid ones.
The paper is organized as follows. In section 2 we introduce the IPC model and specify the changes we perform on the parameters. In section 3 we present the HGD method and provide details about our simulations. Section 4 is dedicated to a thorough discussion of our results. Our final considerations are presented in section 5.

Model potential
Interactions between IPCs are modeled using the coarsegrained description introduced in [5]. An IPC is represented by a spherical, impenetrable hard core of diameter 2σ that carries a negative charge Z c . Concentric to the hard core, a larger sphere of diameter (2σ + δ) defines the interaction range δ. These particles carry two smaller spheres of radius ρ with a positive charge Z p , whose centers are located diametrically opposite to each other at a distance e (eccentricity) from the particle center. This geometry defines two regions on the particle surface: the negatively charged equatorial region (E) and two polar regions (P), each of them carrying a positive charge. A schematic representation of the model is depicted in figure 1.
In this coarse-grained description the interaction between the particles is infinite for distances less that 2σ and zero for distances larger than (2σ + δ). Between these two distances the energy can be calculated using the following expression where for simplicity we have omitted the translational and orientational dependence of U IPC , ω EE , ω EP and ω PP [5]. In this equation¯ = ( EE , EP , PP ) are the constant energy strengths of the EE, EP and PP interactions and ω EE , ω EP and ω PP are the overlap volumes between the corresponding regions of two neighbour particles normalized to the volume of the colloidal sphere. The energy strengths (that are assumed to be constant) are set by mapping the coarse-grained model to the analytical Debye-Hückel potential [5]. Throughout this article, the hard core diameter of the particle (2σ ) is used as unit of length and the minimum energy ( min ) of each model sets the energy unit. This minimum energy is attained when two neighbouring particles are at contact in the EP configuration. Pressure and temperature will be given in the usual reduced energies (p * = p(2σ ) 3 /| min | and T * = k B T /| min |).
The IPC model depends on three adjustable parameters: the interaction range δ, that can be experimentally controlled by modifying the salt concentration of the surrounding solvent, the size of the patches as defined by the half opening angle γ , that can be adjusted by modifying the location of the positive charges (i.e. the eccentricity e) and the charge imbalance between the colloid and the polar caps, Z. Note that two of these parameters are not completely independent and a modification of one of them can induce changes in the other. In particular, γ and ρ depend on e, δ and σ in the following way [5] Moreover changes in either e or δ can lead to considerable modifications of the constant energy strengths EE , EP and PP . As mentioned in the Introduction the equilibrium phase diagram for a set of parameters corresponding to moderately overcharged colloids ( Z = 30) with medium size patches (γ = 38.6 • ) and low screening conditions (δ = 0.25) has already been investigated [10]. This model, whose parameters are listed in table 1, will be referred in this work as the reference Figure 1. Schematic representation of selected IPC models. The hard repulsive core is shown in grey, while the equatorial and polar regions are shown in white and yellow, respectively. The Z60 model (first from the left) has exactly the same geometry as the reference (Ref) model, the only difference being that the colloid charge is increased. In contrast modifications of the patch size and of the interaction range induce also changes on other geometric parameters. For each model, the geometric parameters that are different from those in the reference model are specified by a prime, while dotted lines represent the reference model. The BP model (second from the left) is obtained by decreasing the value of the eccentricity e, which induces a variation of the patch radius ρ and of the half opening angle of the patches γ . The LRI model (third from the left) is obtained by increasing δ while maintaining e fixed, thus leading to a variation of both ρ and γ . Finally, the LRII model (fourth from the left) is built by changing δ and e so that γ is preserved. In the last three models the colloid and patch charges are the same as in the reference model. Opposite trends can be visualized for the SP, SRI and SRII models. (index 'Ref') model. In the low temperature-low pressure region of the phase diagram this model stabilizes a laminar structure, in which IPCs are oriented with their axes parallel to the planes and exposing the equatorial regions between the layers. Within each layer particles are tightly packed forming a grain-like pattern that maximizes the number of attractive EP interactions. On the other hand, the lateral distance between planes can be quite large depending on the thermodynamic conditions. The equilibrium inter-plane distance results from a competition between minimizing the interaction energy and maximizing both the vibrational entropy (thus favouring large distances between planes) as well as the packing fraction (that is obviously attained for short distances between planes). At high pressures, packing effects dominate over the energetics and the layered structure transforms into an FCC structure that is orientationally ordered at low temperatures and becomes a PC for sufficiently high temperatures. The energy of the ordered FCC structure is considerably higher than that of the laminar phase due to the repulsive EE interactions that can be avoided in the laminar phase, but not in the FCC. Further, as for the disordered phases, this model exhibits a stable liquidvapour phase separation with a very low critical pressure (lower than p * = 0.2 [10]). In order to study the effect of the IPC parameters on the phase behaviour, new models are designed in the following by modifying one by one each parameter while keeping the remaining ones constant; the parameters are specified in table 1.
• First, we investigate the effect of charge imbalance (i.e. the difference between colloid and patch charges). In this case the new model is obtained by simply increasing the magnitude of the colloid charge while maintaining the same charges on the patches. The new model, referred to as Z60, is designed by changing Z c of less than 10%, thus obtaining Z = Z c − 2Z p = 60 (for details see table 1). • Secondly, two new models with patch sizes that are different from those of the reference system are constructed by increasing and reducing e (i.e. the value that determines the location of the patches and consequently their size as defined by γ ). These two new models are denoted as big patch (BP) and small patch (SP) models (see table 1 for the numeric values of the model parameters and figure 1 for a schematic representation of the BP model). The half opening angle of the patches in the BP and SP models is about 10% larger and 10% smaller, respectively, than the one of the reference model. To be more specific, the respective patch sizes are defined by γ = 42.38 • for the BP and γ = 33.88 • for the SP model (to be compared to γ = 38.62 • in the reference model). • Third, the effect of the interaction range is also investigated. In the reference model, the Debye screening length is set to κσ = 2 and the interaction range is calculated assuming that δ = κ −1 , leading to δ = 0.25. Two new models with shorter and larger interaction ranges are designed by imposing that κδ = 0.8 and κδ = 1.2, respectively, while keeping the Debye screening length (given by κσ = 2) as well as the rest of the model parameters constant. The new models, denoted by SRI and LRI, have an interaction range of δ = 0.2 and δ = 0.3, respectively (i.e. about 20% lower and 20% higher than the reference model). A schematic representation of the LRI model is shown in figure 1. We also consider an even longer ranged model, referred to as LRIb, obtained with the choice κδ = 1.5. It is worth noting that modifying the interaction range while keeping the eccentricity constant leads to a considerable change of the patch radius and of the patch size (see table 1), as ρ and γ are related through equations (2) and (3). • Finally, we design two models with different interaction range but same patch size with respect to the reference model. This is achieved by solving the system of equations (2) and (3) for the desired values of δ and γ . The SRII and LRII models are designed using this strategy (see figure 1 where the LRI model is shown). With the same procedure we also consider a model characterized by an even longer interaction range, referred to as LRIIb (see table 1).

Hamiltonian Gibbs-Duhem integration
The phase diagrams for the new models were evaluated using the so-called Hamiltonian Gibbs-Duhem (HGD) integration method [20,21]. In the conventional Gibbs-Duhem (GD) method [18,19], the whole coexistence line between two phases can be traced by numerically integrating the Clapeyron equation; only a single coexistence point is required to launch the algorithm. The HGD method, on the other hand, takes advantage of the fact that in theoretical approaches and in simulations it is straightforward to modify, in addition to the thermodynamic variables, also the interaction. It is thus possible to change the interaction continuously from a reference system (with an interaction U ref ), for which the phase diagram is known, to a new system (with potential U new ), for which the phase diagram is to be evaluated. This is achieved, for example, by linearly coupling both model potentials through a parameter λ in the following way By gradually modifying λ one can pass from the reference (λ = 0) to the new system (λ = 1).
Using this coupling parameter λ as an additional intensive thermodynamic variable, a differential change in the Gibbs free energy is given by where X G is an extensive thermodynamic variable conjugate to λ Using the definition of the Gibbs free energy it is obtained or, what is the same where u is the energy per particle. This means that for each value of λ along the integration path between the reference and the target system X G can be evaluated by means of an NpT simulation in which particles interact via the potential U(λ ).
Using equation (4), X G can be rewritten as where u new and u ref are the average energy per particle of the new and of the reference models, and the angular brackets denote ensemble averages over the configurations visited during an NpT simulation for a given value of λ . Two phases I and II are in coexistence when −S I dT + V I dp + X G,I dλ = −S II dT + V II dp + X G,II dλ. (11) This generalized Clapeyron equation can be used to calculate changes in the coexistence line due to an infinitesimal variation in temperature, pressure or interparticle potential. Differential changes of either the temperature or the pressure while keeping the model potential fixed lead to the conventional Clapeyron equation where H I , H II are the enthalpies of the two coexisting phases. Additionally, equation (11) allows us to calculate the change of the coexistence pressure induced by a change in the parameter λ while keeping the temperature constant Alternatively, it is possible to calculate the change in the coexistence temperature with λ when the pressure remains Provided that the coexistence line between two phases is known for the reference potential U ref , the shift in the coexistence line induced by a change of the parameter λ (or, equivalently, of the potential) can be simply calculated by numerically integrating either equation (13) or equation (14). Depending on the slope of the coexistence curve, it can be more convenient to prefer one relation over the other. Note that, like the conventional GD integration, this method is only applicable if the integration path does not cross a phase coexistence region. If new phases emerge for the new model potential that do not occur in the reference system this method is unable to predict their stabilization.
The phase diagram for all the considered models was evaluated using the following procedure. First, we calculated one coexistence point for each of the three phase coexistences (i.e. the fluid-layer, the fluid-PC and layer-FCC ones) using the HGD procedure: to be more specific, equations (13) and (14) were integrated starting from known coexistence points of the reference model (i.e. taken from [10] and summarized in table 2). Usually about six λ-values were used (and turned out to be sufficient) to perform the integration from λ = 0 to λ = 1 (i.e. to pass from the reference to the new system). The integration was carried out using a fourth order Runge-Kutta method: to this end the derivative of either the pressure or the temperature curve with respect to λ (i.e. the right hand side of either equation (13) or equation (14)) was evaluated at each step of the integration scheme in NpT simulations. These simulations extended over 20 000 MC cycles for equilibration plus another 50 000-100 000 cycles for taking averages. Here one MC cycle is defined as N particle displacement attempts plus one volume change attempt. The number of particles in Note: For the SP and SRI models a coexistence value ( 1 ) at a lower temperature was used to avoid crossing the FCC-PC transition along the integration path. These data were taken from [10]. the simulation box was N = 500 for the fluid phase, N = 480 for the laminar phase and N = 500 for the FCC and PC lattices. Secondly, once an initial coexistence point between two phases is known the entire coexistence lines were calculated using the GD integration method, i.e. using equation (12). The FCC-PC transition was evaluated by simply heating and cooling simulation runs, as this phase transition exhibited very little hysteresis. The number of MC cycles and the system sizes used in these simulations are the same as the ones used for the HGD simulations.

Estimation of the critical point
The liquid-vapour critical point was evaluated by grandcanonical Monte Carlo (GCMC) simulations for systems in a box of size L = 11 (in units of the particle diameter). To this end we selected a subcritical thermodynamic state (specified by the chemical potential µ and the temperature T ) for which the system exhibited large fluctuations in density; then we ran ten parallel simulations (each extending over 5 000 000 MC steps) to evaluate the probability distribution of the number density ρ.
The precise location of the critical point was then obtained by searching-using the histogram reweighting technique [24]for µ and T values at which the bimodal probability distribution of ρ matched that of the 3D Ising model.

Results
In sections 4.1-4.3 we report about the effect of the charge imbalance, the patch size and the interaction range on the stability of the layered structure with respect to the fluid, the FCC and the PC phases. In section 4.4 we discuss the effects of the changes in the model parameters on the location of the liquidvapour critical point in the temperature versus density plane.

Effect of the charge imbalance
The effect of increasing the charge imbalance on the translational and orientational dependence of the IPC potential is shown in the top panels of figure 2 for three representative configurations, referred to as EE, EP and PP. For a bigger colloid charge (inducing a higher charge imbalance) the PP interaction is reduced at all distances and for all orientations with respect to the reference model, whereas the EE repulsion moderately increases. The EP interaction does not change by construction. , are used to display the potentials; the color code refers to the representative configurations as labeled. The radial dependence of the potentials is reported for two particles in a fixed mutual orientation as a function of their interparticle distance r (the particle displacement being represented by a black arrow in the corresponding particle-particle configuration); the orientational dependence of the potentials is reported for two particles at contact as a function of the rotation angle θ rot around a given axis (the rotation axis being represented by a dot in the center of the rotating particle and the rotation direction represented by the dotted circle and the arrow in the corresponding particle-particle configuration). Bottom panel: phase diagram of the reference (black lines) and of the Z60 (blue lines) model with phases as labeled. The open circles represent the coexistence pressures at zero temperature between the FCC and the layers for the reference model (black) and the Z60 model (blue); for the two models these points have been evaluated elsewhere with an evolutionary algorithm approach [10,23]. The solid dots signal the liquid-vapour critical point for the reference (black) and Z60 (blue) models.
Before calculating the full phase diagram for this model, we have checked that our implementation of the HGD method is correct and reliable by comparing the results with those obtained combining direct coexistence (DC) method and free energy calculations (see [10] for a detailed discussion) at selected state points. These results, summarized in table 3, show the good agreement between both routes.
The full phase diagram of the Z60 model is shown along with the corresponding data of the reference model in the bottom panel of figure 2. We note that for the Z60 model the zero temperature phase diagram was previously evaluated through an optimization technique relying on a minimization of the enthalpy with a strategy based on the ideas of evolutionary algorithms [23]. The extrapolation to zero temperature of the layer-FCC coexistence line calculated with the HGD method converges to the transition point obtained via the minimization of the enthalpy at zero temperature (see the open blue symbol in figure 2), proving the consistency of our calculations. As discussed in [10] this consistency also holds for the reference system (see the open black symbol in figure 2). The main feature of the phase diagram for the Z60 model is that the range of stability of the layered structure is significantly expanded with respect to the reference model. Interestingly this is mainly achieved by moving the layer-FCC coexistence line to higher pressures, as the fluid-layered coexistence line remains practically unchanged. The stabilization of the layered structure with respect to the FCC solid can be understood on the basis of energetic considerations: as discussed before, the Z60 model is characterized by stronger EE repulsions as compared to the reference system. While in the layered structure these EE repulsive interactions can be avoided by simply increasing the distance between the planes, these interactions are always present in the FCC solid. Thus a stronger charge imbalance leads to a considerable increase of the energy of the FCC lattice, whereas it only increases moderately the energy of the layered structure. The higher energy of the FCC structure explains also why the transition to the PC phase occurs at slightly lower temperatures and why this latter phase is destabilized with respect to the fluid phase.

Effect of the patch size
When the patch size increases with respect to the reference model (BP model), the PP and EP interactions extend over larger angles, the effect being more pronounced in the PP case; on the other hand, the EE repulsion is less extended in its angular dependence. Concomitantly, the relative strength of the PP repulsion is reduced with respect to the EP attraction, whereas the EE interaction is slightly enhanced (see the radial and orientational dependence of the potential in the top panels of figure 3). The reason of this trend is the fact that larger patches are obtained by locating the patch center of charge closer to the particle center; due to the larger distance between the patch charges of two interacting IPCs the strength of the PP repulsion is therefore reduced. These changes in the potential are reversed when the patch size is reduced (SP model). , are used to display the potentials; the color code refers to the representative configurations as labeled. The radial dependence of the potentials is reported for two particles in a fixed mutual orientation as a function of their interparticle distance r (the particle displacement being represented by a black arrow in the corresponding particle-particle configuration); the orientational dependence of the potentials is reported for two particles at contact as a function of the rotation angle θ rot around a given axis (the rotation axis being represented by a dot in the center of the rotating particle and the rotation direction represented by the dotted circle and the arrow in the corresponding particle-particle configuration). Bottom panel: phase diagram of the reference (black lines), BP (red lines) and SP (blue lines) model with phases as labeled. The black open circle shows the coexistence pressure at zero temperature between the FCC and the layers for the reference model; the solid circles the liquid-vapour critical point for the reference (black), BP (red) and SP (blue) models.
The effect of the patch size on the phase diagram is shown in the bottom panel of figure 3. Increasing the size of the patches destabilizes the layered structure with respect to the FCC solid, while it enhances its stability with respect to the fluid. Inspection of the coexistence densities reveals that a larger patch size does not lead to significant changes in the densities of none of the three phases. On the other hand, at constant pressure the energies decrease for all three phases as the patch size increases, although the magnitude of the change is quite different for each case: while for the fluid phase the energy remains essentially the same as in the reference model (it is only about 2% lower in the BP model), the energy of the layered structure is reduced by about 6% and that of the FCC solid by about 12% with respect to the corresponding structure in the reference model. The different energy changes for the different phases can be understood by separately considering the contributions to the interaction energy: the increased angular patch extension leads to a larger negative EP contribution and to a subsequent reduction of the positive EE contribution; on the other hand, the PP energy is reduced due to the larger distance between the patch charges of two neighbouring colloids. The reason why the lowering of the energy is most pronounced for the FCC lattice is probably related to the fact that the reduced EE repulsions are important contributions to the energy of this phase, whereas they are less relevant in the layered structure. On the other hand, the larger average distance between the particles in the fluid phase explains why the energy is less affected in this case. Based on these arguments, also the shift of the layer-FCC and layer-fluid coexistence lines to lower pressures and higher temperatures in the BP model can be understood. The orientational orderdisorder transition of the FCC solid to the PC phase occurs at higher temperatures in the BP model as compared to the reference system: a higher temperature is needed to break the bonds in the FCC and to allow the colloids to rotate freely about their lattice positions. The coexistence line between the fluid and the FCC PC phase, on the other hand, is fairly independent of the patch size, which is most likely related to the fact that particles can rotate in both phases leading thus to comparable energy changes in both phases.
The results for the SP model, characterized by a decreased patch width, show opposite trends with respect to the BP model, thus confirming the scenario already depicted (results are summarized in the bottom panel of figure 3). Specifically, a smaller patch size stabilizes the layered structure with respect to the FCC phase and destabilizes it with respect to the fluid; the order-disorder transition is shifted to lower temperatures. The coexistence lines between the FCC PC and the fluid are not at all affected by the patch size: this line is essentially the same for the reference system and for the BP and the SP models. Note that for the BP model the two triple points nearly coincide.

Effect of interaction range
As mentioned in section 2, an increase in the range of the interaction by simply reducing κδ induces in the model also a larger surface extension of the patches (see table 1). On increasing the interaction range, both the EP attraction and the PP repulsion extend over larger angles and consequently the EE repulsive term vanishes already at shorter angles. The opposite trend is observed for a system with a shorter interaction range. In the top panels of figure 4 the effect of the interaction range on the interparticle potentials is visualized for systems LRI and SRI.
The effect of the interaction range on the phase diagram is shown in the bottom panel of figure 4. For the LRI model, the region of stability of the layered structure in the phase diagram is substantially reduced in favour of a larger range of stability of the FCC structure. On further increasing the The data for the reference, the LRI and the SRI models are depicted by solid, dashed and dotted lines, respectively. Three representative configurations, referred to as PP, EE and EP and reported on the left side of the top panels (from top to bottom, respectively), are used to display the potentials; the color code refers to the representative configurations as labeled. The radial dependence of the potentials is reported for two particles in a fixed mutual orientation as a function of their interparticle distance r (the particle displacement being represented by a black arrow in the corresponding particle-particle configuration); the orientational dependence of the potentials is reported for two particles at contact as a function of the rotation angle θ rot around a given axis (the rotation axis being represented by a dot in the center of the rotating particle and the rotation direction represented by the dotted circle and the arrow in the corresponding particle-particle configuration). Bottom panel: phase diagram of the reference (black lines), LRI (red lines) and SRI (blue lines) model with phases as labeled. The black open dot marks the layer-FCC transition at zero temperature for the reference model. Finally the solid dots signal the location of the liquid-vapour critical point for the reference (black), LRI (red) and SRI (blue) models.
interaction range, i.e. when considering the system LRIb, the layered structure completely disappears from the phase diagram (results not shown). Concomitantly, for the LRI model the fluid looses stability with respect to all the ordered phases, namely the layered structure, the FCC and the PC. For the SRI model opposite trends are observed: the region of stability of the layered structure is considerably extended at the cost of the stability range of the FCC solid while it is destabilized at lower pressures by an extending stability region of the fluid phase. Curiously the layer-FCC coexistence line exhibits a reentrant behaviour (i.e. there is a point at which dp/dT = H /(T V ) = ∞); hence by increasing the pressure at some suitably chosen temperature, the FCC PC phase transforms into a layered structure and then back again to the PC phase.
As in the previous two cases energetic arguments help to understand the reasons why the layered structure is strongly destabilized for longer ranged models. On one side, the fact that the EP attraction extends to larger angles should favour the layered solid (as seen in 4.2), on the other hand, the EE repulsion extends over shorter angles, favouring thus the FCC crystal. Since in the latter phase particles have many nearest neighbours with EE interactions, a reduction in these energies leads to a considerably lower energy. Therefore when the range of interaction increases, the energy of the FCC decreases more rapidly than that of the layered structure. The observed increasing stability of all the ordered phases with respect to the fluid is also due to the stronger EP attractions and the weaker EE repulsions that lower the energy of the ordered phases with respect to that of the fluid. As a result of the lower energy of the FCC solid, its transition line to the PC is shifted to higher temperatures.
Since our procedure to increase the range of the potential implies a change of the patch size, it seems reasonable to infer that the trends observed in the phase diagram emerge as a result of both parameter changes. To decouple these trends, we decided to investigate another possible modification of the reference model that changes the range of the potential while keeping the angular extension of the patches invariant (see section 2); these models will be denoted as SRII and LRII. The radial and angular variations of the potentials of the three reference configurations (EE, EP and PP) for the SRII and LRII models are shown in the top panels of figure 5. By construction, the EE, EP and PP interactions of these two models are very similar to those of the reference model. Of course, the strength of the PP repulsion is considerably enhanced when the range of the interaction increases. This is the result of simultaneously increasing both the interaction range and the eccentricity. In particular, by moving the location of the patch charges closer to the particle surface, the distance between patches belonging to two interacting IPCs is decreased at a given interparticle distance; thus the PP repulsion becomes stronger as compared to its reference counterpart. In the new procedure to change δ, the effect of the interaction range will be thus mixed with the effect of changing the strength of the interactions.
Surprisingly the same qualitative trends in the phase diagram are observed, irrespective of whether the range is increased by keeping the same location of the patches (thus increasing their angular extension) or by simultaneously changing also the location of the patches to guarantee a similar patch size (see the bottom panels of figures 4 and 5). The phase diagrams for the LRI and LRII models are remarkably similar. In both cases the stability of the layered solid is reduced to very low pressures, this effect being more pronounced for the LRI model. This observation is simply reflecting the fact that a larger patch size also leads to a destabilization of the layered structure (as seen in section 4.2). Again, on further increasing the interaction range, i.e. when considering the system LRIIb, the layered structure completely disappears from the phase diagram (not shown here). Larger differences Figure 5. Effect of the interaction range on the phase diagram for fixed patch spatial extension. Top panels: translational and orientational dependence of the selected model potentials. Data for the reference, the LRII and the SRII models are depicted by solid, dashed and dotted lines, respectively. Three representative configurations, referred to as PP, EE and EP and reported on the left side of the top panel (from top to bottom, respectively), are used to display the potentials; the color code refers to the representative configurations as labeled. The radial dependence of the potentials is reported for two particles in a fixed mutual orientation as a function of the interparticle distance r (the particle displacement being represented by a black arrow in the corresponding particle-particle configuration); the orientational dependence of the potentials is reported for two particles at contact as a function of the rotation angle θ rot around a given axis (the rotation axis being represented by a dot in the center of the rotating particle and the rotation direction represented by the dotted circle and the arrow in the corresponding particle-particle configuration). Bottom panel: phase diagram of the reference (black lines), of the LRII (red lines) and of the SRII (blue lines) model with phases as labeled. The black open dot marks the layer-FCC transition at zero temperature for the reference model. The solid dots signal the location of the liquid-vapour critical point for the reference (black), LRII (red) and SRII (blue) models.
in the stability ranges are observed between the SRI and SRII models. In particular, we observe that for the SRI model the layered structure is more destabilized with respect to the fluid phase but more stabilized with respect to the FCC (again as a consequence of the smaller patch sizes, see section 4.2); in addition, the layered structure remains stable up to higher temperatures in the SRII model.

The liquid-vapour critical point
We investigate the effect of the changes in the model parameters specified in section 2 on the location of the liquidvapour critical point in the temperature versus density plane. Results for the values of the critical temperature, T c and the critical density, ρ c , for all models are summarized in figure 6, where the respective data for the reference system are also reported. As in this contribution we are more interested in the ordered phases, we have not evaluated the pressure at the critical point, p c . A rough estimate for p c obtained for the reference system suggested that its value is below p * = 0.2 [10]; the small perturbations of the model parameters implemented in this investigation should leave this scenario unchanged. Thus we expect for all systems investigated the liquid-vapour critical point to be located in the low pressure region of the pressure versus temperature plane.
The condensation of the liquid phase in patchy systems is mainly related to the available bonding volume that a particle has. For IPC systems, the number of configurations in which two particles can form a bond is imposed by the complex interplay of geometric terms (related to the patch size and to the particle interaction range) and energetic factors (related, in turn, to the competition between the attractive and repulsive parts of the pair potential). As a consequence of this complex scenario, we are not able to infer general considerations from the specific systems studied here, nonetheless a few trends can be extracted from our investigations.
First we observe that by increasing the charge imbalance with respect to the reference model (i.e. for model Z60), the liquid-vapour phase separation is disfavoured: both the critical temperature, T c and the critical density, ρ c , become smaller than the corresponding values in the reference system. We observe that while the effect on the temperature is rather moderate (T c decreases by ≈2.5%), the change in density is quite significant (ρ c decreases by ≈9%).
Similarly, when the patch size is decreased (SP model) the critical point moves towards lower temperature and density values; in this case variations are smaller: T c decreases by ≈1.6%, while ρ c is reduced by ≈3.4%. The opposite trend is observed for bigger patch size (BP model), with smaller variations in the temperature-T c increases by ≈0.8%-and larger changes in the density-ρ c increases of ≈4%. Note that in both models, the change in the critical temperature is smaller than the one in the critical density. Interestingly, figure 3 shows that the shift of the critical temperature to higher values with the patch size goes in the same direction as the shift of the fluid-layer coexistence on passing from system SP to system BP, although the magnitude of the change in the coexistence temperature seems to be larger in the case of the fluid-layer transition. As a consequence, the critical point, in spite of being located at higher temperatures, might be metastable with respect to solidification if the size of the patches is further increased from the BP model.
Finally, when changing the interaction range the effect on T c becomes more important: although the critical density shows substantial changes with δ, the critical temperature shows even stronger variations. We can summarize that short ranged potentials lower the critical temperature, while long ranged potentials increase it. In particular, when the range is changed together with the patch size, model SRI has a critical density that is practically unchanged as compared to the reference system, while its critical temperature is reduced by ≈8%; on the other hand, model LRI shows a strong variation both in T c and ρ c of ≈9%. When the range is changed together with the position of the patch center of charge, model SRII shows a variation of ≈5% for T c and ≈6% for ρ c ; in contrast, in model LRII the critical density is essentially unaffected while its critical temperature is ≈7% higher than the corresponding value in the reference system.
We note that data referring to systems I exhibit a minimum in density-located roughly at the position of the reference model-and that systems II show the same qualitative behavior; nonetheless the line connecting the critical points of systems I and that connecting the critical points of systems II are displaced and cross each other at the reference model; this demonstrates how the changes in the range and in the patch size add up in a consistent way for the set of data specified by index I. Indeed, models I (LRI and LRIb) are characterized by the same interaction range but bigger patches with respect to models II (LRII LRIIb); therefore, their critical densities move to higher values (in agreement with model BP); on the other hand, model SRI has a smaller patch size than model SRII and thus its critical density is lower (in agreement with the SP behavior).
Finally our results indicate that the liquid-vapour phase separation could go from being stable to metastable with respect to the solid-fluid phase separation when the range of the interaction increases. Similarly to the trends observed for the BP model, both coexistence lines shift to higher temperatures for larger interaction ranges, but the effect is bigger for the fluid-layer coexistence line. As a consequence the liquidvapour critical point might become metastable, specially for the LRI model (for which the larger range and bigger patches effects add up).

Conclusions
In the present contribution we numerically evaluated the phase diagram for a selection of inverse patchy colloids with two positively charged polar patches and an oppositely charged equatorial belt. We started from a system of slightly overall charged particles that are decorated by relatively extended and long-ranged patches. Previous studies have shown that the phase diagram of this system is characterized by a broad region where a crystal formed by parallel monolayers is stable [10]; in addition, at high pressure values and low temperatures the system forms an FCC solid which transforms with increasing temperature into a plastic FCC lattice. In the low pressurehigh temperature region the system forms a fluid phase with a liquid-vapour transition, characterized by a very small value for the critical pressure.
In the present work, we investigated the effect of changes in the model parameters on the stability of the layered phase with respect to the fluid, the orientationally ordered FCC and the plastic crystal. Additionally, we show how the liquidvapour phase separation (quantified by the location of its critical point) can be favoured or disfavoured on changing the model parameters. In particular, we focused on the impact (i) of the charge imbalance between the differently charged surfaces of the colloids (i.e., the poles and the equator), (ii) of the patch size and (iii) of the particle interaction range. The strategy of how to change the model parameters is rather delicate since the aforementioned parameters are related to each other within the coarse-grained description of the model. For instance, a change in the interaction range can be accompanied by a change in the patch extension or in the effective patch interaction strength. This complexity is also reflected in experimental systems of heterogeneously charged units: if considering complex units emerging from the adsorption of charged polyelectrolyte stars onto the surface of oppositely charged colloids [25], a change in the salt concentration of the solution surrounding the particles leads to a partial detachment of the polyelectrolyte stars, thus inducing a change in both the interaction range and the patch size; on the other hand when salt is added to a suspension of heterogeneously charged model colloids synthesized in the lab [26], the patch size is not affected, while the surface charge-and hence the patch effective interaction strength-can be.
We can summarize our observations as follows: • Overcharging the particles in favour of the bare colloid charge stabilizes the layered structures at the cost of the ordered FCC crystal. Although one might naively think that stronger interactions lead to a higher T c , overcharging the particles leads to a slight decrease in T c , which is due to the complicate interplay between the attractive and repulsion terms. A stronger effect is observed for ρ c , which moves to lower densities as the colloid charge increases. • Smaller patch widths stabilize the layered structure with respect to the FCC lattice, but destabilize the laminar phase with respect to the fluid; obviously, wider patches induce opposite trends. Larger patch sizes lead to higher T c and ρ c , the effect being more pronounced on ρ c rather than on T c , for which only slight changes are observed. As a consequence our results indicate that the liquid-vapour phase separation might become metastable with respect to the formation of the laminar phase as the patch size increases, although further calculations would be needed to confirm this scenario.
• The same qualitative trends in the phase diagram are observed irrespective of whether the range is increased by keeping the same location of the patches (thus increasing their angular extension) or by simultaneously changing the location of the patches (thus changing the patch effective strength) to achieve a similar patch size. For shorter interaction ranges the region of the phase diagram corresponding to the layered phase increases considerably due to its stabilization with respect to the FCC (it looses some stability with respect to the fluid phase). The stabilization of lower density ordered structures with respect to dense solids for lower interaction ranges has also been observed for other usual patchy models [27,28], which leads us to speculate that this might be a general trend irrespective of the form of the interparticle potential. As with regards to the liquid-vapour separation, T c increases with the interaction range, as expected. However ρ c shows a non-monotonic behaviour, whose origin is again probably due to the mixed effects of changing two features of the model. Indeed this nonmonotonic effect is more pronounced in models for which the effect of changing both the interaction range and the patch size add up. Even though the range of stability of the liquid-vapour separation decreases for larger interaction ranges, it might remain always stable as increasing further the interaction range leads to the complete destabilization of the layered solid.