Phase behaviour of colloidal superballs mixed with non-adsorbing polymers

Inspired by experimental work on colloidal cuboid-polymer dispersions (Rossi et al., Soft Matter, 7, 4139 (2011)) we have theoretically studied the phase behaviour of such mixtures. To that end, free volume theory (FVT) was applied to predict the phase behaviour of mixtures of superballs and non-adsorbing polymer chains in a common solvent. Closed expressions for the thermodynamic properties of a suspension of hard colloidal superballs have been derived, accounting for fluid (F), face-centred cubic (FCC) and simple cubic (SC) phase states. Even though the considered solid phases are approximate, the hard superballs phase diagram semi-quantitatively matches with more evolved methods. The theory developed for the cuboid-polymer mixture reveals a rich phase behaviour, which includes not only isostructural F1-F2 coexistence, but also SC1-SC2 coexistence, several triple coexistences, and even a quadruple-phase coexistence region (F1-F2-SC-FCC). The model proposed offers a tool to asses the stability of cuboid-polymer mixtures in terms of the colloid-to-polymer size ratio.


Introduction
While colloidal hard spheres only exhibit an equilibrium fluid-solid phase transition, the phase behaviour of anisotropic hard particles is more intricate. Directional excluded-volume interactions between anisotropic colloidal particles also give rise to entropy-driven phase transitions [1,2], which have been studied for lyotropic systems such as rod-like [3][4][5] and platelet-like particles [6], as well as for non-axisymmetric colloids [7,8]. The immersion of non-adsorbing free polymers (depletants) into an ensemble of hard particles leads to a region around the colloidal particles where the polymers have a lower configurational entropy than in bulk, namely the depletion zone (bounded by the grey curves around the black shapes in fig. 1). The excluded volume between the depletant and the colloidal particle of interest defines this depletion zone. Overlap of these depletion zones (grey areas in fig. 1) leads to an effective attraction between the colloidal particles [9][10][11]. Depletant addition to a system of anisotropic particles induces effective depletion attraction patches [12][13][14][15] because the depletion attraction is stronger for larger overlap of depletion zones (see bottom panels of fig. 1). The effects of entropic patchia e-mail: r.tuinier@tue.nl ness [12,13] in lyotropic systems have received increasing attention [14,16]. However, analytical treatments for entropic patchiness between non-axisymmetric particles are scarce.
Among the possible non-axisymmetric particles, colloidal cuboids are of interest due to their potential application as photonic crystals [17,18], their possible roles in emulsion stabilization [19] and anti-reflective coatings [20] and to prepare porous membranes [21]. Due to recent progress in colloidal synthesis, it is nowadays possible to prepare colloidal cuboids with a well-defined shape and size [22,23]. A commonly applied model to describe colloidal cuboids is the superball shape. Formally, superballs are a subset of a family of geometric shapes called superellipsoids, introduced by Barr [24]. The implicit equation describing the shape of a superball [25] reads where a is the radius of the superball (the shortest distance from the centre of the superball to its surface) and m is the shape parameter. The surface of the superball is described for f(x, y, z) = 1, whereas the locus of points inside the superball are retained for f(x, y, z) < 1. Here we focus on m ≥ 2: the shape of the superball lies in between a sphere (m = 2) and a cube (m = ∞) [26]. We depict in fig. 2 a collection of superballs in the range of m values investigated.
The phase behaviour of colloidal superballs has been studied both experimentally [18,27] and via computer simulations [26,28,29]. Some experimental studies on the effect of non-adsorbing polymers on the phase behaviour of colloidal superballs have also been conducted [22,30]. However, a (relatively) simple model for the thermodynamic properties of superballs (and superball-polymer mixtures) is not available. Hence, we first describe in this paper a theoretical framework both for the superball fluid and for two possible superball solid states. The tools required to calculate the phase diagrams of the cuboidpolymer mixture are subsequently presented. With all closed expressions at hand, we present a collection of phase diagrams as well as a phase stability overview, which reveals various rich multi-phase coexistence regions, including a possible four-phase coexistence. We finalize by summarizing the main conclusions of our findings. We provide a reproducible, closed model whose insights may provide a better understanding of experimentally observed phase behaviour of superball-polymer mixtures [22,30]. 2 as a function of the shape parameter m. Gray dots correspond to numerical solutions of the superball area and surface mean curvature (see appendix A for details). The black curve shows a fit through the data points, whose accumulated relative error is 5.15 Npoints is the number of points used to find the fit in eq. (3).

Theory
In this section we first present the canonical expressions developed for a suspension of colloidal superballs, both in the fluid and in the two solid phases considered. Secondly, we present the free-volume theory developed to account for the grand-canonical thermodynamic properties of a mixture of superballs plus non-adsorbing polymers.

Fluid state
We consider a collection of N c hard superballs in a volume V , each superball having a volume v sb , surface area s sb and surface integrated mean curvature [31] c sb . The second virial coefficient (B 2 ) for hard particles is given by the orientationally averaged excluded volume between two particles [32], and for a suspension of monodisperse convex particles (hence for superballs with m ≥ 2) in a fluid state reads [33,34] where γ is the so-called asphericity parameter. The numerically obtained normalized second virial coefficient B * 2 for hard superballs is shown in fig. 3. In appendix A we explain how γ can be computed numerically, yielding B * 2 using eq. (2). It follows that B * 2 smoothly increases with m from the sphere limit (m = 2, B * 2 = 4) to the cube limit (m = ∞, B * 2 = 5.5) due to the increase of the particle anisotropy. In this work, we use a closed expression for B 2 by fitting (solid curve) the calculated data (see appendix A for details) with the (inverse) equation of an ellipse, which accurately matches the data (accumulated error shown in fig. 3): An equation of state (EOS) for a fluid of hard convex particles was first derived by Gibbons using scaled particle theory [35], and a more accurate EOS was proposed by Boublík [36][37][38] taking into account virial coefficients higher than B 2 in a Carnahan-Starling-like fashion [39]: where Π o is the reduced osmotic pressure of the pure hard superball dispersion, β is 1/k B T (with k B the Boltzmann constant and T the absolute temperature), and where φ c is the volume fraction of hard superballs: The superscript "o" is used to indicate the (depletant-free) osmotic pressure and chemical potential of a suspension of pure superballs. Using the closed relation between B 2 , γ and Π o of eqs. (2)-(4), the EOS for a fluid of hard superballs as a function of m is completely defined. Computer simulations have shown the accuracy of the Boublík EOS for a wide range of m values [28,29]. Obviously, the Carnahan-Starling [39] EOS is recovered for m = 2. The chemical potential of the superballs is related to the osmotic pressure through the Gibbs-Duhem relation for a single-component system at constant temperature: with μ o = βμ o the reduced chemical potential. The chemical potential follows from eqs. (4) and (5) as with μ ref = ln(Λ 3 /v sb ) the reference chemical potential of a superball fluid and Λ the thermal wavelength. The free energy follows from the chemical potential and the osmotic pressure through the thermodynamic relations: with F = βF v sb /V the reduced free energy. The approximated free volume is illustrated as a grey region. Also indicated are the nearest-neighbour distance (r) and the maximum distance from the centre of the superball (r 2D max ) for m = 10. Gray shapes represent a particle that just touches a nearest neighbour. Bottom panels as top ones but for the simple cubic (SC) lattice. In this case the arrows hold for r and the superball radius (a) as these quantities defined the volume that the central particle explores without overlapping with its nearest neighbours. In all cases considered here the solids are 20% off their close-packing structure.

Solid states
For the colloidal solid phases we modify the cell theory proposed by Lennard-Jones and Devonshire (LJD) for hard spheres [40]. We consider each particle to be contained in a closed region whose shape is determined by its neighbouring particles, which are fixed at their lattice positions [41], see fig. 4. The free energy of the solid is calculated from the number of configurations determined from the volume V f that the centre of the particle explores without overlapping with its nearest neighbours. This leads to the following normalized free energy for a solid: The free volume V f depends on the shape parameter m and the volume fraction φ c , but also on the relative position of the nearest neighbours and hence on the structure of the solid. In this work we consider two crystal structures: facecentred cubic (FCC) and simple cubic (SC). A schematic view of the FCC and SC structures of superballs for several m values is shown in fig. 4. We focus first on the FCC crystal. The exact free volume depends on the shape of the Wigner-Seitz cell [41], which for an FCC crystal has a rather complicated geometry [42,43], but is usually approximated as a sphere [11]. The free volume considering this spherical approximation is given by where r is the distance between the centres of a superball and its nearest neighbours, and r cp is r at close packing. For the FCC crystal, we consider a "frozen" crystal where the particles are perfectly aligned. Hence, for the FCC lattice r cp is two times the distance between the edges of the superballs (2r 2D max , see appendix A). The distance r at a certain volume fraction can be determined from r cp as with φ cp c the close-packing fraction. Combining r cp = 2r 2D max = 2a √ 2(1/2) 1/m (see appendix A) with eqs. (8)-(10) provides the free energy for the FCC phase, using a Taylor expansion for the term (φ cp,FCC c /φ c ) 1/3 − 1 as in the original LJD approach for hard spheres [11,40]. The chemical potential and osmotic pressure are calculated via the thermodynamic relations given in eq. (7), leading to the following closed (yet approximate) thermodynamic expressions for the FCC phase as a function of m: The m-dependency of the close-packing volume fraction in an FCC crystal is provided as (see appendix B) where f (m) is (see appendix A) where Γ is the Euler Gamma function. For m = 2, eq. (11) recovers the free energy for hard spheres in the FCC phase [40]. The thermodynamic properties of the SC superball crystal is found using a similar approach as for the FCC crystal. For the SC structure, the free volume has the shape of a cube. To approximately take into account the effect of rotations of cuboidal particles on the SC free energy, the size of the free volume is chosen such that the central particle can rotate inside its unit cell. This effectively reduces the free volume by a factor of eight compared to the one obtained for perfectly parallel cuboids. This approach overestimates the free volume for particles with low m values. The free volume of the SC crystal is given by with r defined via eq. (10) and r cp = 2a. Following a similar procedure as for the FCC phase state, the thermodynamic functions of the SC phase read with the close-packing fraction in the SC phase given by We note here that effects of particle rotations are only qualitatively accounted for via our estimation of V f . For non-axisymmetric hard particles, Onsager-like theories [3] for crystalline phases are non-trivial due to the lack of a single reference axis for the inter-particle orientations. Already for biaxial hard particles no analytical solutions are found [44]. Further discussion on possible ways of incorporating rotational contributions is presented in appendix C.
In table 1 we provide the close-packing volume fractions for perfect spheres (m = 2) and perfect cubes (m = ∞) and for a limiting intermediate case. While FCC packings are more efficient for small m, SC arrangements can pack closer for large m. It follows from eqs. (12) and (16) that both the FCC and the SC phases have the same φ cp c at m = 3. Further details on the close-packing fractions are provided in appendix B.
Most of the limitations in our model are those inherent to the cell theory used in the calculation of the free energy of the crystalline phases. Cell theory is known to give accurate results for FCC and SC crystals of spherical colloids [42,43], but extending cell theory to other crystal structures is not straightforward due to the complex geometries of the space explored by the centres of mass of the particles. Already for a body-centred-cubic crystal of hard spheres, cell theory does not match with simulations [43]. Furthermore, the cell theory approach followed does not account for defects, which affect the fluid-crystalline phase transition of cuboidal hard particles [45,46]. The accuracy of cell theory for anisotropic particles is a matter of debate. Particularly, it is difficult to accurately account for rotational contributions for non-axisymmetric particles into the solid-phase partition function. More advanced theoretical approaches to properly describe the complex solid phases of superballs are rather involved and would make the description at hand less tractable (and most likely not involving a set of closed expressions).

Free-volume theory
We account for the mixtures of superballs plus nonadsorbing polymer in a semi-grand-canonical fashion via free-volume theory (FVT) [11,47]. Within FVT, the cuboid-polymer system (S) is considered to be in equilibrium with a reservoir (R) of polymers. In R and S the solvent is treated as background. The system and the reservoir are connected through a membrane permeable for the polymers and the common solvent, but impermeable for the cuboidal particles. A sketch of this osmotic equilibrium approach is given in fig. 5. The key quantity relating the polymer concentrations in R and S is the free-volume fraction available for depletants in the system α. Following the original ideas of FVT, we assume that α is independent of the chemical potential of the depletants in R. Furthermore, we take the simplest model for polymeric depletants, namely the penetrable hard-sphere (PHS) model [48]: depletants are treated as ghost-like spheres (with radius δ) that can freely interpenetrate each other but do not overlap with the superballs. PHSs mimic ideal polymer chains [11]. One last approximation made is that the ensemble-averaged free volume for depletants in the system, V free , is independent of the concentration of depletants: V free ≈ V free o . This results in the following (normalized) expression for the grand potential of the system [16]: with V the volume of the system, Since the depletants are considered to behave ideally, the (reduced) osmotic pressure in the R is simply given by Van 't Hoff's law: The depletant concentration in the system follows from From the semi-grand potential, the chemical potential of the superballs and the osmotic pressure of superball plus PHS mixtures are obtained through the standard thermodynamic relations: with N R d the number of depletants in R. Using these quantities, phase coexistences follow from where i and j denote the two (or more) coexisting phases (fluid, FCC or SC). If coexistence between three phases takes place a triple point (TP) arises, and we denote fourphase coexistence via a quadruple point (QP) [16,49]. Colloidal systems may exhibit isostructural phase coexistence (such as gas-liquid equilibrium [50,51]) when attractive interactions between particles are present [47,52,53]. In such a case, the low-density phase will be entropically favourable and the high-density phase is stabilized by attractive interactions between the particles. The limit of isostructural phase coexistence is defined via the critical point (CP), at which Whenever a phase state has a CP an isostructural phase coexistence can take place, which can be stable or metastable. The transition from a potentially stable to a metastable isostructural phase coexistence is defined by the critical end point (CEP), where the CP and the TP or QP of the corresponding isostructural coexistences merge [52,53]: These conditions enable to determine the topology of the phase diagrams as a function of the system parameters: the colloidal shape (through m) and the relative depletant size through

Free-volume fraction
The only remaining unknown parameter in eq. (17) is the free-volume fraction for depletants in the system, α.
Widom's insertion theorem [54] relates the free-volume fraction α to the work (W ) required to bring a depletant from R to S via where V free o is the free volume for depletants in the undistorted (depletant-free) system. This work (W ) is approximated using Scaled Particle Theory (SPT) [55,56], by connecting the limits of inserting a very small depletant (up to second order) and a very big depletant in the system of interest, followed by scaling back to the actual size of the depletant. As the depletants are spherical, a single scaling factor (λ) enables to express this work as In the small depletant insertion limit (λ 1) no overlap of depletion zones is assumed. This allows writing W as a function of the excluded volume between a superball and a sphere (v exc ), which actually defines the depletion zone volume. This leads to where the scaled depletion volume is obtained by scaling the depletant thickness: δ → λδ. We use normalized units also in the big-depletant limit for convenience: By combining eqs. (24) to (27), a general expression for with v exc = v exc /v sb and A closed expression for v exc to compute α from eqs. (23) and (28) is convenient. The depletion zone due to spherical depletants around any hard, convex body is embedded within the Connolly surface [57] around the particle of interest with thickness δ. An apparently simple expression is available for the excluded volume between general, convex bodies and spheres [31]. The excluded volume between a superball and a sphere reads where s sb = s sb /d 2 and c sb = c sb /d. Equation (29) can however not be solved analytically (see appendix A for details on the calculation of s sb and c sb ). Due to the linear relation between δ and q (δ = qa), v exc (λ) is simply obtained by taking q → λq in eq. (29). Via interpolation of s sb and c sb it is possible to obtain a (non-closed) expression for eq. (28), hence completely defining the grand potential Ω (eq. (17)). Alternatively, we found that the depletion zone is accurately described by a tilted dampening sinus function: with p = 1 − 2/m. An advantage of eq. (30) over eq. (29) is that it has a closed form. Furthermore, the deviation between eqs. (29) and (30) is very small: (N points being the number of points used to fit eq. (30)). Thus, α is incorporated in our calculations via inserting eq. (30) into eq. (28). Moreover, the final thermodynamic properties of superball-polymer mixtures as considered here do not depend on the approach followed for calculating v exc (see appendix D). All results presented were computed using Mathematica [58].

Results and discussion
In the present section, we first provide the results for the ensemble of hard superballs in the absence of depletants. The free-volume fractions for depletants in this system is then briefly discussed, which provides all components to understand the phase diagrams of the superball-polymer mixtures of interest. Based upon these phase diagrams, an overview of the rich multi-phase coexistences exhibited is presented in a single, comprehensive plot.

Phase diagram of hard superballs
The calculated phase diagram for a suspension of pure hard superballs is presented in fig. 6 (left panel), and compared with more evolved simulation results (right panel). The well-known fluid-FCC coexistence for hard spheres [11,59] is recovered at m = 2, and it slightly shifts to higher densities with increasing m. As can be appreciated, this is in agreement with computer simulation results [26]. This shift can be attributed to a thermodynamically less favourable FCC state with increasing m. The forbidden region (in grey) simply identifies densities beyond the close packing of the considered phases. The discontinuity on our phase diagram corresponds to m = 3, which locates the transition between preferred FCC to more favourable SC. The thermodynamically preferred phase is roughly the one with the largest closepacking fraction, see table 1. However, at intermediate pressures, crystal structures with lower close-packing fraction may be also stable. We find a triple F-FCC-SC coexistence at m ≈ 3.71. Triple-phase coexistence for pure hard superball system has also been reported based upon Monte Carlo simulation results (right panel of fig. 6) [26]. Between m = 3 and m ≈ 3.71 we find SC-FCC coexistence, which arises due to the different φ c -dependencies of the solid equations of state (for more details, see appendix B). Above m ≈ 3.71, only F-SC coexistence is found, which shifts towards lower packing fractions with increasing m, also in qualitative agreement with simulations [26]. In the cube limit (m = ∞), we find φ F c ≈ 0.36 and φ SC c ≈ 0.54. Simulation studies [29] indicate that for perfect cubes phase coexistence between a fluid and a cubatic liquid crystal takes place at φ F c ≈ 0.47 (a phase showing high orientational order but no long-range translational order), while at higher densities a transition into an SC crystal occurs at φ SC c ≈ 0.58. This points towards a complex nature of the F-SC phase transition for perfect cubes that can not be accounted for with our simple theory.
The overall topology of the simple, theoretical phase diagram corresponds to the one found using more evolved computer simulations (right panel of fig. 6) [26,28,29]. The differences can be justified because we do not account for the same solid phases for superballs as in simulations. Jiao et al. [25] showed that the optimal packing for superballs are retained in the C 0 (low m values) and C 1 (high m values) crystalline phases. Both the C 0 and C 1 lattices are obtained via deformation of the FCC (m = 2) and SC lattices (m = ∞), respectively. In fact, these solid phases are accounted for in simulations [26,28]. Not surprisingly, the triple point from simulations is a fluid-FCC-C 1 [26]. Due to the limitations inherent to the simple model developed, the C 0 and C 1 phases are not accounted for. The FCC phase features (and its coexistences) roughly match those of the plastic FCC. The role played in simulations by the C 1 is mimicked by the SC phase in our simpler model.
When considering colloidal cuboid-polymer mixtures, the phase diagrams would only enrich upon refinements of the method. The liquid-crystalline and crystalline coexistence regions are found in simulations in a broader range of m values. Based on experimental observations [22,30], we expect the depletion attraction to enhance solid phases where overlap of the depletion zones is maximized (leaving space for the depletants to fit in the voids of the respective lattices). Note however that these experimental observations correspond to colloid-polymer mixtures confined at a surface, whereas the results here presented holds for bulk systems.

Free-volume fraction
We show examples of free-volume fractions for PHSs in a fluid state of hard superballs (α F ) using W as in eq. (30) in the left panel of fig. 7. It follows that α F only weakly depends on m. For m > 2 and low q values (q 0.05) the free-volume fraction is always slightly smaller as for the sphere case (m = 2). However, with increasing q the intricate shape of the depletion zone around a superball causes α F to be higher than for hard spheres.  table 1). For m = ∞, the solid phase with the higher α S is the SC one. The sudden vanishing of α S (clearly observed for all m = 2 and m = ∞ in fig. 8) corresponds to the close-packing fractions of the FCC and SC lattice at each m value. Within SPT, the work required to bring a depletant into a close-packing state is infinitely large.

Phase diagrams
Firstly, we consider a superball whose shape is still close to a sphere. We present phase diagrams of superballs with m = 2.5 and added depletants for three relative size ratios q in fig. 9. The depletant-free baselines (φ R d = 0) for the fluid-FCC coexistence correspond to the densities shown in fig. 6 (left panel). Upon addition of depletants, the FCC phase at coexistence gets denser and the coexisting fluid phase becomes more dilute in order to maximize the total free volume available for the depletants in the system. At sufficiently large q values (q = 0.4 and q = 0.6 in fig. 9), an isostructural colloidal F 1 -F 2 (also termed gas-liquid) coexistence appears (which is metastable for low q values). We do not further address metastable coexisting phases. This F 1 -F 2 coexistence spans until the F 1 -F 2 and the F-FCC coexistences match: at this φ R d value a triple line is found (upper panel for q = 0.4), which becomes a region in the system representation (lower panels of fig. 9). In fact, an isostructural phase coexistence is always connected to a triple coexistence when more than one phase identity is considered. When q increases it can be seen that the F-S coexistence narrows and the F 1 -F 2 critical point shifts to higher depletant volume fractions.
The coexistence regions in the system representation (bottom panels in fig. 9) show that the fluid phase with a low concentration of superballs has a high concentration of depletants, whereas the FCC phase has a high concentra- Fig. 9. Phase diagrams for a mixture of colloidal superballs for m = 2.5 and PHS depletants for several relative depletant sizes q as indicated in the reservoir representation (top panels) and system representation (bottom panels). The F 1-F2 critical point is indicated by a black dot. Triple-phase coexistences are shown as a horizontal black line in the reservoir representation and as a coloured area bounded by black lines in the system representation. All coexisting phases present are indicated in the reservoir representation. Insets in the system representation zoom in on the low depletant concentration region, and some of the coexistence regions are indicated. A few illustrative tie-lines are shown as dashed grey lines for q = 0.2. Above each phase diagram, a 2D illustration of the superball (black) and its depletion zone (grey) are shown, and the m and q values are indicated.
tion of superballs but a low concentration of depletants. The incorporation of partitioning of depletants over the different phases is one of the key elements of FVT [11]. The system representation also shows that for a superballdepletant mixture a single solid phase (without a coexisting fluid phase) only occurs at nearly imperceptible depletant concentrations. So far we observe no special features compared to FVT for hard spheres mixed with penetrable hard spheres [47], even though the colloidal shape considered is not perfectly spherical. For hard-spheres plus polymeric depletants, simulation results taking into account multi-body effects show that the trends predicted by the simple FVT hold [60].
Next, we show in fig. 10 phase diagrams for m = 3.33, where FCC-SC coexistence was found for hard superballs (see fig. 6, left panel). At sufficiently high depletant concentration, the FCC phase becomes metastable w.r.t. the SC phase, as expected because φ cp,SC c > φ cp,FCC c for this m value: the FCC phase completely disappears, and an F-FCC-SC triple point is always found. For sufficiently large q values (q = 0.4 and q = 0.6) two triple-phase coexistences are present whenever F 1 -F 2 coexistence is found: F 1 -F 2 -FCC and F-FCC-SC, with the corresponding triplepoint areas in the system representation. This illustrates the enrichment of the phase diagram topology due to depletion effects.  Phase diagrams for even more cubic particles, m = 5, are presented in fig. 11, for which only an SC state is present in the pure hard superballs system (see left panel of fig. 6). Similar qualitative trends as in fig. 9 are observed, but with the F-SC coexistence playing the role of the F-FCC equilibrium. For small q values the broadening of the coexistence lines occurs at lower depletant concentrations with respect to the superball-polymer mixture with m = 2.5 in fig. 9: the overlap of depletion zones is larger for particles with an increased cubicity, which results in a stronger depletion attraction. For m = 5 and q = 0.4, there is no F 1 -F 2 equilibrium phase coexistence, whereas F 1 -F 2 coexistence was found at this q value for superballs with m = 2 [47], m = 2.5 and m = 3.33. Due to the tendency of flat faces to align upon addition of depletant into the system, stable F 1 -F 2 coexistence shifts to higher q values: longer ranges of attraction are required for larger m values to induce stable F 1 -F 2 coexistence.

Multi-phase coexistence overview
The rich multi-phase coexistence behaviour hinted at in the previous section is quantified by calculating the critical end points (CEP) of all possible isostructural coexistences as a function of the system parameters m and q. The calculated CEP curves are summarized in fig. 12, which constitute the main result of our investigations on the phase behaviour of superball-depletant mixtures. The limiting values of the FCC phase (m 3.71) and the SC phase (m 3) are indicated as vertical dashed lines in fig. 12: for m ∈ {3, 3.71}, F-FCC-SC coexistence always takes place. To the left of this m interval the only solid phase found is FCC, and to the right the only solid state found is the SC. In fig. 12, the solid curves hold for the CEPs defining the limiting q values at which (stable) isostructural phase coexistences are found. For sufficiently high q values, F 1 -F 2 isostructural coexistence is expected for all m values, which spans from an F 1 -F 2 -FCC triple region or from an F 1 -F 2 -SC triple coexistence as indicated.
To gain insight about the F 1 -F 2 isostructural coexistence and the connection with a triple F 1 -F 2 -FCC or F 1 -F 2 -SC coexistence we show in fig. 13 phase diagrams for q = 0.4 and a few selected m values (moving along a horizontal line in fig. 12). For m = 3.4, an F-FCC-SC triple coexistence occurs at a higher φ R d than the F 1 -F 2 -FCC triple line (F 1 -F 2 -SC coexistence is metastable). On the other hand, for m = 3.65 the F 1 -F 2 -FCC triple point becomes metastable and an F-FCC-SC triple point arises at lower φ R d than the F 1 -F 2 -SC isostructural coexistence. This is explained by the fact that the stability of the FCC decreases as m increases. The condition at which the F 1 -F 2 -FCC and the F 1 -F 2 -SC coexistences merge results in a quadruple coexistence (F 1 -F 2 -FCC-SC). This four-phase coexistence is present for a range of m values, and shows an asymptotic behaviour from the CEP of the quadruple coexistence towards the pure hard superball triple point (see fig. 12). The shift of the QP from the TP of the depletant-free system towards lower m values reflects the patchiness of the depletion attraction between superballs. As a consequence of the enhanced alignment of the flat faces upon addition of depletants, F-SC coexistence takes place at m values below those of the depletant-free system. Hence, depletion-mediated entropic patchiness promotes the appearance of the SC phase. We conclude that quadruple coexistence arising from merging two isostructural triple-phase coexistences is possible for superball-polymer mixtures. Recently, we found that similar quadruple coexistences can appear in plateletdepletant mixtures [16,61]. In fact, quadruple coexistence containing isotropic-isotropic coexistence has been experimentally reported for platelet-polymer mixtures [62] and for a hydrate forming system [49]. As an F-FCC-C 1 triple point is present from simulations for the depletant-free superball system [26], the corresponding quadruple-phase coexistence may be found from simulations with the C 1 phase instead of the SC phase used here. However, hinting at this rich phase behaviour directly from simulations may be too demanding without the results provided by the simple model presented. Four-phase coexistences are possible in effective two-component systems provided an extra field variable [61]: in this case, the superball's shape parameter m.
A remarkable finding is the appearance of an SC 1 -SC 2 isostructural coexistence found at low q and high m values (see lower right part of fig. 12). We depict a few illustrative phase diagrams in fig. 14, where small isostructural SC 1 -SC 2 coexistence regions appear. The single-phase fluid and simple cubic regions get smaller upon decreasing q. For m = 10 the binodals shift towards lower φ R d values with decreasing q, following the same trend as observed for the F 1 -F 2 , F-FCC and F-SC coexistences (see figs. [9][10][11]. This solid-solid coexistence is driven by the entropic gain for depletants upon phase separation of the colloids into a dense solid SC phase and a more dilute SC phase. The low q values at which this coexistence takes place are related to the low (yet non-zero for cubic enough superballs) freevolume fraction for depletants available in solid phases at high colloid concentrations (see fig. 8). As can be observed in the rightmost panel of fig. 14, the m value tunes the depletant concentration at which SC 1 -SC 2 coexistence is found: SC 1 -SC 2 equilibria are driven by the alignment of the flat faces, and thus for more curved particles (decreasing m) the SC 1 -SC 2 coexistence requires a higher depletant concentration. This induces the crystal state to demix into an attractive solid (depletion zones optimally overlapping) coexisting with a repulsive one, as expected for short-ranged attractions in colloidal systems [63,64]. With more accurate models or in experimental systems, these SC 1 -SC 2 coexistences may be replaced for example by a C 1 coexisting with an SC phase. The absence of a stable FCC 1 -FCC 2 coexistence can be rationalized by the non-optimal overlap of depletion zones between the flat faces of the superballs in an FCC state.

Concluding remarks
A simple model for the thermodynamic properties of superballs was presented, where all thermodynamic functions required for the phase diagram calculation are expressed in closed form. We account for three phase states: fluid (F), face-centred cubic (FCC) and simple cubic (SC). Some of the closed expressions provided, such as the accurate fit for the second virial coefficient for hard superballs, may be of direct application not only in theoretical, but also in experimental studies. Despite the assumptions made, the found phase behaviour of pure hard superballs semi-quantitatively recovers the trends observed as compared with Monte Carlo (MC) computer simulations. Further improvements of the theory developed may be possible, particularly in the solid phases, but most likely lacking the simple, closed expressions reported here. The increase of the excluded volume with increasing particle anisotropy as well as the tendency for flat faces to align allowed us to rationalize the phase diagram obtained in terms of the colloidal packing fraction and the particle shape parameter.
The addition of free, non-adsorbing polymers to a collection of hard superballs induces effective attractive patches: alignment of the less-curved areas of the superballs is enhanced. This is clearly reflected in the presence of F-SC phase coexistences for shape parameters where it was not stable in the depletant-free system, and also in the manifestation of SC 1 -SC 2 phase coexistences for sufficiently cubic colloids upon addition of small depletants. Such solid-solid coexistences may be not only of fundamental relevance, but also of relevance for the designing novel photonic crystals. When the depletion attraction is sufficiently long-ranged, isostructural F 1 -F 2 coexistence is found for all shape parameters, which can coexist either with an FCC or an SC state depending on the superball shape. The boundary between these two triple points that included isostructural fluid phases (F 1 -F 2 -FCC and F 1 -F 2 -SC) defines a window for quadruple-phase coexistence (F 1 -F 2 -FCC-SC). The main trends were collected in a simple, comprehensive plot ( fig. 12) containing the system parameters (colloidal shape and relative depletant size) that summarizes the effectively patchy nature of the depletion attraction in colloidal suspensions of cuboidal particles. The system parameters at which rich multi-phase behaviour of hard-superballs polymer mixtures is revealed with our simple model may guide more accurate computer simulations, and may serve as a first qualitative guide to interpret experimental results.

Author contribution statement
RT proposed the research and did some preliminary calculations on the superball system. AGG developed the model for superball-polymer mixtures, which was revised and improved by JO. All authors jointly wrote and discussed the manuscript.

Appendix A. Superball properties: calculation
Superball properties required for the calculation of the second virial coefficient and for the free energy of the solid phases are detailed in this appendix. The distance r between the centre of a superball and an arbitrary point on the surface of the superball is given by [65] r(θ, φ) = a (| cos φ| m | sin θ| m + | sin φ| m | sin θ| m where θ and φ are the polar angle and the azimuthal angle, respectively. The maximum distance (r max ) between the centre and the surface of the superball is the distance from the centre to the corner, as shown in fig. 2 for the 2D superball projection. The angles corresponding to this maximum distance are θ = π/4 and φ = 0 for a 2D superball and θ = arccos( 2/3) and φ = π/4 for a 3D superball, which leads to a maximum distance given by The volume of the superball V sb is obtained by integration of eq. (A.1) [65]: where integration is performed over one octant due to symmetry. Equation (A.4) can be solved analytically, resulting in [21,26] with d the diameter of the superball (d = 2a) and Γ the Euler Gamma-function. Exact equations for the surface area s sb and for the mean curvature c sb of a superball are not known, but they can be calculated numerically using the surface integral and the integral of mean curvature [65]: where subscripts denote partial derivatives and x represents a vector from the centre of the superball to the surface, given by

Appendix B. Close packing and free volume of FCC and SC crystals
In this appendix, clarification on the close-packing fraction of the two solid states considered is provided. The general equation for the close-packing fraction of a superball crystal is given by with N uc the number of superballs in the crystal unit cell and V uc,cp the volume of the unit cell at the close-packing fraction.
For the FCC crystal, the number of particles inside the unit cell is 4 and the volume of the unit cell at the close-packing fraction is given by which, combined with eq. (B.1), gives the close-packing fraction of superballs in the FCC crystal (eq. (12)). For the SC crystal, there is only a single particle inside the unit cell and the volume of the unit cell at the closepacking fraction is simply given by which results in the close-packing fraction of superballs in an SC crystal given by eq. (16). In fig. 16 the close-packing fractions for the FCC and SC superball crystal are shown as a function of the shape parameter m. As the particles become more cubic, the close-packing fraction of the FCC lattice decreases, whereas the close-packing fraction of the SC lattice increases. We also show in fig. 16 the normalized free volume for the FCC and SC crystal as a function of the packing fraction. For the FCC crystal, the relative free volume decreases when m increases, which seems counterintuitive since the distance between particles r increases as a function of m for a constant packing fraction due to the fact that the size of the particles increases as a function of m. However, due to the diagonal stacking of the particles, the free volume per particle decreases as a result of the increase of r 2D max as a function of m. The opposite trend is observed for the SC crystal, which is as expected because the free volume of the SC lattice does not depend on r 2D max , but decreases as a function of the particle radius a which does not depend on the m value. These trends in the free volume are also visible in the 2D projections of the free volume in the FCC and SC crystals shown in fig. 4.

Appendix C. SC crystal for perfect cubes with approximated contribution for rotations
Here, an approximation for the rotational contribution of cubes (m = ∞) to their free volume in an SC crystal is discussed beyond the particular choice of the unit cell. Instead of reducing the free-volume cage by a factor 8, Fig. 17. Sketches of the rotation approximation for cubic particles in an SC crystal in 2D (left) and 3D (right). Dashed lines indicate the SC unit cells, two nearest neighbours are shown in black and the free volume is indicated in grey. In the 2D scheme the minimum distance between two nearest neighbours (grey line), the angles θ and θ/2 and the distances 2a sin θ and 2a cos θ are also indicated.
as was done in the main text, we rescale the free volume based on the maximum angle that a cube can rotate without entering the unit cell of a neighbouring particle (see fig. 17).
The maximum angle that a cube can rotate in one direction inside its unit cell θ can be found by solving The average angle of a cube inside its unit cell is approximated by θ/2. Note that cuboids can rotate in three dimensions: here we only consider rotations around one of their symmetry axes. The free volume is approximated as two times the minimum distance between two particles, rotated in one direction with opposite orientation, to the power of three. Following this approach, the free volume of a cube in an SC crystal V SC The maximum angle, the free volume, and the free energy as a function of the packing fraction are shown in fig. 18. The angle θ goes from zero at the maximum packing fraction to π/4 at the situation where cubes freely rotate in one direction without touching each other (at φ c = 1/(2 √ 2)). Note that this is not the same maximum packing fraction at which cubes in an SC crystal can freely rotate in all directions, namely φ c = 1/(3 √ 3). Taking into account this approximate orientation of the cubes leads Fig. 18. Left: maximum angle a cube can rotate in one direction in an SC unit cell θ, centre: free volume V free and right: normalized free energy e F of an SC crystal of cubes, with the rotation approximation shown in fig. 17 and for an SC state of perfectly aligned cubes. to a reduction of the free volume for all volume fractions. This free volume is not a monotonically decreasing function with φ c , which shows the intricate of accounting for rotations even qualitatively. Following this rough rotation approximation, the F-SC coexistence densities for hard cubes are found at φ F c ≈ 0.37 and φ SC c ≈ 0.58. These densities are very similar to those found with the appropriate choice of the unit cell made in the main text: φ F c ≈ 0.36 and φ SC c ≈ 0.54. Note that even this simple approach, that only takes an approximated average particle orientation in two dimensions into account, may only be analytically resolved for perfect cubes. However, this simple derivation allowed us to estimate how much the free-volume cage is reduced due to particle rotations in an SC state: namely about 1/8 of the free volume for perfectly aligned cubes.

Appendix D. Effect of different depletion zones on phase stability
Naively, the depletion zone may be considered as having the same shape as the hard superball but whose radius is a + δ: v eq.shape exc = (1 + q) 3 . (D.1) Derivation of the free-volume fraction from eqs. (28) and (23) is straightforward. In table 2 we show the differences for the CEP of the quadruple F 1 -F 2 -FCC-SC coexistence. Differences between eq. (29) and eq. (30) are only in the third decimal place. When the depletion zone is overestimated (eq. (D.1)) the CEP of the quadruple coexistence occurs at a lower range of the depletion attraction. We finally present in fig. 19 the phase stability overview. A nearly imperceptible difference is found when applying eq. (29) or eq. (30), whereas the F 1 -F 2 coexistence occurs at lower q when applying eq. (D.1). As the Fig. 19. Same as fig. 12, but comparing different methods for the depletion zone calculation. Black curves corresponds to the interpolation of the numerically calculated superball properties (eq. (29)), the orange dashed curves that practically coincide with the black ones correspond to the fitting used (eq. (30)) and the gray solid curve correspond to the equal shape approximation for the depletion zone (eq. (D.1)).
SC 1 -SC 2 coexistence only occurs at rather low q values, it is less sensitive to the particular shape of the depletion zone considered.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.