Nonadditive interactions and phase transitions in strongly confined colloidal systems

The behaviour of colloids can be controlled effectively by tuning the solvent-mediated interactions among them. An extensively studied example is the temperature-induced aggregation of suspended colloids close to the consolute point of their binary solvent. Here, using mean field theory and Monte Carlo simulations, we study the behaviour of colloids confined to a narrow slit containing a nearly-critical binary liquid mixture. We found that the effective interactions in this system are highly non-additive. In particular, the effective interactions among the colloids can be a few times stronger than the corresponding sum of the effective pair potentials. Inter alia, this non-additivity manifests itself in the phase behaviour of confined colloids, which depends sensitively on the slit width and temperature. In addition, we demonstrate the possibility of a first-order bridging transition between colloids confined to a slit and suspended in a phase-separated fluid well below the critical point of the solvent and at its critical composition in the bulk. This transition is accompanied by a remarkably large hysteresis loop, in which the force between the colloids varies by two orders of magnitude.


Introduction
Colloids attract increasing attention from scientists and engineers as models for atoms and molecules 1 and due to numerous technological applications in the food industry, 2 medicine, 3 and other fields. [4][5][6][7] The behaviour of a colloidal suspension can be effectively controlled by tuning the properties of its solvent, such as the concentration, temperature, and composition. An outstanding example is the coagulation of colloids close to the critical point of the solvent, which gives rise to the so-called critical Casimir forces induced by order parameter fluctuations. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] Correspondingly, the phase behaviour of such a colloidal system exhibits well-defined colloidal, liquid and gas phases. Additionally, a number of crystalline phases 23,24 and glassy states 25 have been identified.
As a result of confinement, most physical processes deviate from their bulk behaviour, and may exhibit even new physical phenomena. The corresponding plethora of phenomena includes, e.g., unusual drying transitions 26 and freezing 27,28 of confined water, an anomalous increase of capacitance in subnanometer pores, 29,30 a remarkable long-ranged repulsion between confined van der Waals dimers, 31 etc. As we shall report, remarkable phenomena emerge also in strongly confined colloidal systems. Such systems have attracted much attention recently, [32][33][34][35] in particular as models to study the behaviour of particles in (quasi) two spatial dimensions. 36,37 Previous experimental and theoretical studies have focused on such systems with short-or long-ranged repulsive interactions in the context of the Kosterlitz-Thouless transition. [38][39][40][41][42][43] However, identical colloids in a narrow slit filled by a nearly critical solvent acquire strong long-ranged attractive interactions, which can be manipulated by minute changes of temperature or solvent composition. This allows one to study two-dimensional systems with tunable attractive interactions, and it may open up the door for novel types of twodimensional photonic crystals to serve as sensing devices. 44 The strong colloid-colloid interactions near the consolute point of a solvent are due to an increased correlation length of the solvent and due to modifications of its order parameter fluctuations caused by the presence of the colloids (for instance concentration fluctuations in a binary liquid mixture or density fluctuations in a simple fluid). 45,46 By analogy with the quantum Casimir effect, where the confinement of zero-point fluctuations of the electromagnetic fields leads to a measurable effective force between the confining surfaces, 47 the solvent-mediated force close to the consolute point of the solvent is called a critical Casimir force. Similar to their quantum-mechanical analogues, the critical Casimir interactions can be attractive or repulsive depending on the boundary conditions, viz. the force is attractive for like and repulsive for unlike boundaries. 48 For a binary liquid mixture, the two types of boundary conditions (say '+' and 'À') are determined by the preference of the surfaces for one of the two solvent species. Since confined near-critical binary liquid mixtures belong to the Ising bulk universality class and to the so-called normal surface universality class, such surface preferences map onto fixing the boundary spins to the values +1 or À1. 49,50 In the aforementioned context of the collective behaviour of the colloids, an important feature of critical Casimir forces is their non-additivity. Previous theoretical studies have focused on nonadditivity in systems consisting of three colloids 51 and of two colloids near a wall, 52 using a numerical mean field approach. 53,54 In a very recent experimental study 55 the non-additivity of critical Casimir interactions has been observed for three colloidal particles. The measured colloid-colloid Casimir potential was modified by the presence of a third colloid, as compared to the system with two colloids only. In these systems, however, the contribution due to non-additivity constitutes less than 20% of the total Casimir interaction, which is difficult to capture experimentally or numerically due to potentially significant inaccuracies.
Here, we shall demonstrate that the critical Casimir interactions in a confined colloidal suspension are highly non-additive. It is evident that if such interactions were additive, then the lateral force between the colloids would be unaffected by the presence of the slit walls. We shall show, however, that this is not only not the case, but that the non-additivity can strongly change the colloid-colloid interactions. This non-additivity gives rise to a sensitive slit-width dependent phase behaviour of the confined colloids, which can be used to probe and visualize the non-additivity of critical Casimir interactions experimentally.
While in the critical region of the solvent the colloid-colloid interactions are determined by strong long-ranged correlations of the fluid, a different physical mechanism takes place below the critical point of the solvent, where two distinct homogeneous phases can coexist. For example, in the case of a binary mixture of species A and B, these phases are rich in A and B; or they are the liquid and gas phases in the case of a simple fluid. In a system in which one of these phases is thermodynamically favoured in the bulk, it is possible that the other phase is preferred by the colloidal surfaces. In this case a wetting layer forms around a single colloid, which may extend to form a capillary bridge between two (or more) colloids. [56][57][58][59] Such capillary bridges may undergo so-called bridging transitions and give rise to strong capillary forces, which can be used to construct temperaturecontrolled colloidal networks and microstructures. 60 Here, we shall show that in a confined colloidal system the capillary forces and bridging transitions can be observed even in the case of no preference in the corresponding unconfined system and that they can be effectively controlled by the confinement.

Solvent-mediated interactions between confined colloids
Although we do not restrict our considerations to the critical region of the solvent, it will be convenient to discuss the solvent-mediated forces comprehensively in terms of scaling functions. For two identical colloids centered in the symmetry plane of a slit (see Fig. 1), we represent the force between the colloids as reported in ref. 61 and 62 (for the notations see Fig. 1) where T is the temperature, k B the Boltzmann constant, t = (T À T c )/T c the reduced temperature relative to an upper critical point T c , and AE corresponds to t _ 0. The bulk correlation length is given by x AE = x AE 0 |t| Àn , with the nonuniversal amplitudes x AE 0 forming a universal ratio x + 0 /x À 0 , and the critical exponent n = 0.63002(10) E 0.63 in spatial dimension d = 3 (see ref. 63) and n = 1/2 in d = 4, i.e., within mean field theory. In the critical region |t| { 1, F AE is a universal scaling function defined as the critical Casimir force expressed in terms of k B T per colloid radius R. Thus, we shall interchangeably call F AE a force or a scaling function, depending on the context. For notational simplicity, we shall also use F in place of F AE and Y in lieu of Y AE as long as it does not lead to confusion.
Integrating f over the colloid-colloid separation gives the critical Casimir potential f, which likewise can be presented in the scaling form The scaling function F AE is the critical Casimir potential in units of k B T. Similar to the force, in order to avoid a clumsy notation, we shall omit AE in F AE where appropriate.
In the limit O = W/R -N one might expect that the scaling functions F and F reduce to the scaling functions for two Fig. 1 Two identical spherical colloids with the same radius R confined to a narrow slit filled with a nearly critical fluid. W is the slit width and D is the surface-to-surface distance between the colloids. The boundary conditions on the slit walls and colloid surfaces can be '+' or 'À'. For a binary mixture, a '+' (or '-') surface means that one (or the other) of the two species is preferred by the surface. We consider generic symmetric configurations [a(b)a], where a = 'AE' and b = 'AE' correspond to the wall and colloid surface, respectively. In this case the colloids position themselves in the middle of the slit. Our main interest lies in the symmetric configuration with opposite boundaries between the walls and the colloids, as shown in the figure. Note that the configurations [+(À)+] and [À(+)À] are equivalent because the fluid is taken to be at its critical concentration. colloids in a bulk critical fluid, i.e., F(Y,D,O = N) = F bulk (Y,D) and F(Y,D,O = N) = F bulk (Y,D), which have been extensively studied in the literature. 19,61,62,64 Therefore, the bulk scaling functions lend themselves as a reference to compare with, and deviations from F bulk and F bulk serve as a measure of nonadditivity, as shown in previous analyses. 51,52,55 We shall use mean field theory (MFT) and Monte Carlo (MC) simulations (of the Ising model) in order to calculate the scaling functions F and F for the system shown in Fig. 1. We consider a generic symmetric configuration, [a(b)a], where a and b (= '+' or 'À') denote the boundary conditions at the surface of the slit walls and of the two identical colloids, respectively. Such symmetric configurations can be realized experimentally by placing colloids into a binary liquid mixture (e.g., water/lutidine) confined between parallel glass plates. [65][66][67] First, we shall study how such a narrow slit confinement alters the bulk critical Casimir interactions between the colloids, and thus manifests non-additivity. Then, we shall consider two important ramifications resulting from the confinement and from non-additivity.
Reduction, enhancement, and screening of colloidal interactions by confinement Mean field predictions. Within mean field theory, the equilibrium state of a system corresponds to the minimum of the standard Landau-Ginzburg-Wilson (LGW) Hamiltonian; here, the LGW Hamiltonian has been minimized numerically for the configurations of Fig. 1 by using F3DM 54 (see Methods and Section S1 in the ESI † for details). These calculations predict that, for all boundaries of the [+(+)+] configuration being the same, the strength of the colloid-colloid interactions is reduced by confinement (Fig. 2). This is understandable because the walls promote a fluid environment preferred by the colloids and thus the strength of the interactions, determined by the non-uniformity of the mean field profile, decreases. Furthermore, in this case the non-additivity contributes 20% or less to the critical Casimir force. This is similar to the results already reported for other geometries, 51,52,55 for which the effect of non-additivity is weak or at most moderate.
However, the critical Casimir force shows a qualitatively different behaviour for the opposite boundaries of the [+(À)+] configuration ( Fig. 3(a and b)). The corresponding scaling function starts to deviate from the scaling function F bulk in the bulk for Y = Y AE = R/x AE = sgn(t)|t| 1/2 R/x AE 0 t 5. This deviation increases strongly as the scaled temperature Y decreases towards negative values, i.e., T o T c (Fig. 3a). Remarkably, upon increasing the slit width, the scaling function F À does not approach F À bulk , unlike in all other cases (see the inset in Fig. 2 as well as Fig. 3b).
Although at first glance this behaviour might seem surprising, it has a simple explanation. Indeed, for T o T c (i.e., t o 0) there are two stable states with the order parameter values j AE ¼ AEA ffiffiffiffiffi ffi Àt p , both of which minimize the bulk free energy in the limits h -AE0, respectively, where h is the external bulk field conjugate to the order parameter (see h in eqn (3) in Methods), and A is the amplitude of the bulk order parameter. In the unconfined system (i.e., free boundary conditions at infinity), the value of the order parameter far from the colloids (j + or j À ) follows the preference of the colloid surfaces, while in confinement that value is set by the slit walls. This is borne out by Fig. 3c, where we plot the order parameter in the middle of the slit at a distance 2R outside the colloids (along the x axis, see Fig. 1). For Y t À3, the curves for O = W/R = 3 and O = 10 practically coincide, while for the unconfined (bulk) system the order parameter has the opposite sign and is consistent with the boundary conditions at the colloid surfaces. Correspondingly, the forces in wide slits (O -N) and in the unconfined bulk system differ substantially in that the wide slits generate the presence of a +/À interface. Fig. 3(a) and (b) show that for the [+(À)+] configuration the critical Casimir force between colloids is stronger for stronger confinement. However, this only holds if the distance D between the colloids is small. Indeed, Fig. 3d demonstrates that for large D (\ 2R) the force is significantly larger for wider slits. In other words, confinement screens the critical Casimir interactions at large colloid-colloid separations. This surprising behaviour can be related to the formation of an order-parameter 'bridge' linking two colloids at small distances ( Fig. S1 in the ESI, † see, cf., also Fig. 6a). [56][57][58]68 Pulling the colloids away from each other will eventually cause this bridge to break at a certain distance D bridge . Without the bridge the order parameter coronas around each colloid overlap barely and hence the interaction between the colloids is much weaker. D bridge decreases for narrowing the slit because the slit walls are opposite to the colloids and enforce the bridge to break. Therefore, for thin slits the decrease of the force occurs at small D, while the force remains strong for slits wide enough to keep the bridge intact. For instance, for Y À = À2 in Fig. 3d we find D bridge /R E 2 for O = W/R = 3, while D bridge /R \ 4 for O = 10 (not shown).
Interestingly, for wide slits and Y o 0 the force between the colloids practically coincides with the force obtained for the unconfined systems upon applying a weak external bulk field h 4 0 (O = 10 and O = N, respectively, in Fig. 3d). This is View Article Online because in this case the bulk field sets the value of the order parameter far from the colloids to be opposite in sign to the boundary conditions at the colloid surfaces. This leads to the formation of order-parameter bridges connecting the two colloids and hence to quantitatively comparable forces. We shall discuss 'bridging' and bridging transitions in more detail below. Before, however, we analyze how our MFT results compare with the corresponding Monte Carlo simulation data.
Monte Carlo simulations. The above mean field results for the universal scaling functions are quantitatively reliable for spatial dimensions d 4 d u , where the upper critical dimension for the Ising universality class considered here is d u = 4 (with logarithmic corrections in d = 4). 50,69 Moreover, the MFT results are the first term in a systematic expansion of universal quantities in terms of e = d u À d. In this case the MFT results provide an approximation for the critical behaviour in d = 3. Accordingly, in order to judge our mean field predictions, we have performed suitable Monte Carlo (MC) simulations of the lattice Ising model using the method developed in ref. 66 for quasi-spherical colloidal particles (see Methods and Section S2 in the ESI †).
We emphasize that due to computational limitations the systems studied via MC simulations had to be relatively small. For instance, the diameter of the colloids was 11 lattice constants (corresponding to 11 spins) and the distance between the surfaces of two colloids in some cases was as small as just two spins. (Larger distances up to 30 spins have also been considered.) Such small systems are expected to show some degree of non-universality, which is the case also for our MC results (unlike the present version of MFT, which focuses on the universal behaviour). It seems that there is no alternative for removing this non-universal contribution other than by simulating larger systems. Therefore, we present our unprocessed simulation data. However, by comparing the MC results for various system sizes, we have checked that these data are dominated by the universal contribution (see Fig. S3 and Section S2A in the ESI †).
In view of the above discussion, nevertheless the MFT and MC results are consistent with each other. Fig. 4(a) shows that the depth of the Casimir potential can reach values as large as À20k B T, which is more than 6 times stronger than the maximum of À3k B T obtained in the bulk system. Remarkably, for low temperatures (T o T c ) the Casimir interactions are stronger for wider slits, while narrow slits provide stronger interactions close to T c (Y { 1). Similarly, for small colloidcolloid separations (D t R), the interactions are stronger for thinner slits, but thin slits screen the Casimir interactions more effectively if D is large (Fig. 4b). This is in line with the corresponding MFT results (Fig. 3d) and can be related to the breaking of a bridge formed between colloids as the distance between them increases ( Fig. S1 and S4 in the ESI †). The order parameter j in the middle of the slit at a distance 2R outside the colloids (in the x direction, see Fig. 1) as a function of temperature for various slit widths; j is normalized to the amplitude A of the order parameter in the bulk system without colloids. In panels (a)-(c) the surface-to-surface distance between the colloids is D = D/R = 0.5. (d) Casimir force scaling functions versus surface-to-surface distance D = D/R between the colloids. The force is normalized to the amplitude C þ;þ of the critical Casimir force at T c for a slit with (+,+) boundary conditions. The case O = N has been obtained in a bulk system (i.e., no slit walls) by applying a weak external bulk field, with the sign opposite to the boundary conditions on the colloid surfaces, in order to mimic the presence of the slit walls infinitely far from the colloids (see the main text and Methods). This scaling function practically coincides with the scaling function calculated for O = 10. In all plots the centers of the colloids are located in the symmetry plane of the slit (Fig. 1).

Phase behaviour of confined colloids
In the following we use the critical Casimir potential obtained from our Monte Carlo simulations in order to study the phase behaviour of confined colloids. We consider only narrow confinements, in which colloids position themselves in the symmetry plane of a slit (Fig. S6 in the ESI †). In order to calculate the phase diagram of such an effectively two-dimensional system, we have used the random phase approximation (RPA) in two dimensions (see Methods and Section S3 in the ESI †). We note that within this approach we take into account the nonadditivity due to the slit walls, but we ignore the many-body Casimir interactions among the colloids. Although the latter contribution is expected to be small, 51,52,55 it will nevertheless be interesting to analyse such effects in future studies.
As an example, we consider colloids of radius R = 250 nm (ref. 17) suspended in a 3-methyl-pyridine (3MP)/heavy water mixture, which has the bulk critical temperature T c E 312 K and the correlation length amplitude x + 0 E 1.5 nm. 18,19 We consider charged colloids so that the effective hard disk radius is s E 1.05R due to a short-ranged electrostatic repulsion (Section S3 in the ESI †). The resulting phase diagram is shown in Fig. 5 in the plane of the two-dimensional colloidal packing fraction Z and the deviation DT = T À T c from the consolute point T c of the solvent in bulk. As in d = 3, there is a region in the phase diagram where the colloidal gas phase and the colloidal liquid phase coexist. This region shrinks as DT increases and ends in a colloidal critical point (Z cc , DT cc = T cc À T c ). The colloidal critical temperature deviation DT cc increases for narrower slits because confinement enhances the critical Casimir interactions (for small colloid-colloid separations, see Fig. 4b and Fig. S5 in the ESI †). Screening of interactions at large separations is not sufficient to significantly change this behaviour; however, it generates only a relatively moderate dependence of the colloidal phase diagram on the slit width. Nonetheless, taking into account that current cooler/heater systems can control temperature with a precision of about 2 mK (see, e.g., ref. 55), this slit-width dependence of the colloidal phase diagram appears to be experimentally well accessible.
It is instructive to estimate the effect of non-additivity on the colloidal phase diagram. To this end we recall that if the interactions were additive, the lateral interactions between colloids would not be altered by the presence of the slit walls. Thus, we take the critical Casimir potential of the bulk (unconfined) system and calculate the corresponding two-dimensional phase diagram. (Physically this means that the colloids should be confined to the midplane by alternative means, for instance by optical tweezers; to which extent this can be implemented experimentally has to be tested.) For the 3MP/heavy water mixture discussed above, we find that the colloidal critical point of this system would be at DT cc E 0.53 K, which is significantly closer to the bulk critical point of the solvent (i.e., DT = 0) than the system of colloids in the actual slit confinement (see O = 3.1 and 4.5 in Fig. 5).
We thus conclude that confined colloids present a simple and convenient system for probing experimentally the occurrence of  (Fig. 1). Note that F for the wide slit (O = 12) differs from F bulk , in accord with the MFT results ( Fig. 3b and d).  (Fig. S7 in the ESI †). The phase diagrams (symbols) have been calculated using the random phase approximation in two dimensions. The colloid-colloid interaction potentials have been obtained from the three-dimensional Monte Carlo simulations of the Ising model. The diagrams are calculated for colloids of radius R = 250 nm, while the hard disk radius is s = 1.05R, taking into account the short-ranged repulsive electrostatic interactions between actual charged colloids (see Methods and Section S3 in ESI †). The solid line and the diamonds denote the phase diagram of a hypothetical system with the bulk critical Casimir potential taken for the effective colloid-colloid interaction, which would be realized for the slit confinement if the interactions were additive. non-additivity of critical Casimir interactions via measuring the slit-width dependence of the colloidal phase diagram. Fig. 3d and 4b show that for T o T c the colloid-colloid interactions decay rapidly as the separation between them increases. We have related this behaviour to the breaking of a bridge formed between the colloids at small separations ( Fig. S1 and S4 in the ESI †). Such capillary bridges in near-critical fluids have been previously investigated both theoretically [56][57][58][59][70][71][72][73] and experimentally 74 for colloids in bulk systems [56][57][58][59][72][73][74] and for inhomogeneous slits without colloids. 70,71 Here, motivated by the remarkable behaviour of the force (Fig. 3d and 4b), we study the temperatureinduced bridge formation between two colloids in a strong slit confinement using mean field theory and Monte Carlo simulations of the Ising model (see Methods).

Bridging transitions for confined colloids
Our mean field calculations predict that the formation and the breaking of bridges can occur via first-order (discontinuous) phase transitions. Fig. 6(a) shows that the configurations with and without bridge can coexist, i.e., both configurations can be stable or metastable at the same temperature, slit width, and colloid-colloid separation.
In this context, a remarkably large hysteretic loop of the solvent-mediated force is observed as a function of temperature (Fig. 6b). This hysteresis occurs because the system follows different metastable branches, with and without bridges, as temperature increases and decreases. The force is very strong if there is a bridge between the colloids (almost two orders of magnitude stronger than the force at the critical point), but it practically vanishes if the bridge ruptures. Interestingly, the free energies of these two states are surprisingly close to each other (Fig. S2 in the ESI †), so that it is even difficult to determine the transition point (as the intersection of the free energy branches) with certainty, and hence it is not shown in Fig. 6b. However, this also means that the bridge and no-bridge configurations are de facto equally stable and exhibit quasicoexistence in a wide range of temperatures. Taking as an example the 3-methyl-pyridine/heavy water mixture (critical point T c E 312 K and correlation length amplitude x + 0 E 1.5 nm, see ref. 18 and 19), and colloids with 0.5 mm in diameter, from Fig. 6b we roughly estimate the size of this coexistence region to be 0.74 K, starting at approximately 0.36 K below T c . We regard this to be experimentally accessible. (Note, however, that this estimate follows from MFT, accounting for fluctuations is expected to modify the size and the location of the coexistence region.) Although we do observe the formation of bridges in our MC simulations ( Fig. 7a and Fig. S4 in ESI †), we have not found sufficiently strong arguments to identify the bridging transition. In order to check this possibility, we have studied -in Ising language -the magnetization m 0 at the mid-point between the two colloids, which can serve as an order parameter distinguishing the bridge and no-bridge states (m 0 4 0 in the no-bridge and m 0 o 0 in the bridge state for the [+(À)+] configuration). In the case of a first-order transition, m 0 should jump at the transition. However, our simulations suggest that m 0 varies continuously as a function of Y (Fig. 7b), which implies a continuous transformation between the bridge and no-bridge configurations. This is likely due to the fact that the colloids in our simulations are too small to accomplish the transition. Indeed, our colloids have a radius of 5.5 to 7.5 lattice constants, which roughly corresponds to a few nanometers and thus is comparable with the correlation length at the suspected transition. This means that even thermal fluctuations can render a metastable bridge state unstable. In ref. 72 and 73 it has been demonstrated that for a bulk system with a square-well fluid-fluid interaction potential (with a range of 1.5a, where a is the diameter of the molecules of the solvent) the bridging transitions can only take place if the colloidal radius is larger than E10a, which translates to roughly 10 lattice constants in the simulated Ising model. With Fig. 6 Bridging transition and hysteresis for two confined colloids corresponding to the [+(À)+] configuration obtained from mean field theory. (a) Normalized order parameter (OP) distribution for two coexisting configurations with and without a bridge between two identical colloids. The OP is normalized by its amplitude A in the bulk (unconfined) system with no colloids. The same parameters are used in the top and in the bottom panel, viz. the slit width W = 3R, where R is the colloidal radius, the surface-to-surface distance D = 0.8R between the colloids, and the temperature scaling variable Y = À20. (b) The normalized attractive force between two colloids as a function of scaled temperature Y = sgn(t)R/x(t). The arrows indicate the directions in the hysteresis loop, and the black diamonds denote the temperature used for the plots in (a). Concerning the formation and breaking of a bridge along the hysteresis loop, see Supplementary Video (ESI †). Note that within MFT the discontinuity of the force upon a bridging transition corresponds to a break in slope of the effective potential due to the crossing of the free energy branches for the bridge and no-bridge states. Concerning the locus of the bridging transition see the corresponding discussion in the main text.
the currently available computational resources we have not been able to simulate such large systems. We leave this task for future studies.
Finally, we point out that, according to the line of arguments provided by ref. 56, the bridging transition between spheres involves only a quasi-zero dimensional volume of ordering degrees of freedom, which smears out this first-order phase transition. This finite size effect has to be taken into account for the interpretation of MC data for large systems.

Conclusions
Using mean field theory and Monte Carlo simulations of the Ising model, we have found a remarkably strong non-additivity of critical Casimir interactions in confined colloidal suspensions. This follows from noting that if the assumption of additivity were correct, the lateral interactions between colloids in a slit would not be affected by the walls. In sharp contrast, however, we have found that for small colloid-colloid separations the critical Casimir interactions can be enhanced through confinement several times over. At larger separations, in contrast, the colloidcolloid interactions are strongly screened as the slit becomes narrower ( Fig. 3 and 4).
This non-additivity manifests itself in the phase behaviour of confined colloidal suspensions. Similar to a bulk (i.e., unconfined) system, there are liquid-like and gaseous colloidal phases, which can coexist in certain ranges of temperature and colloidal packing fraction. The coexistence region ends in a critical point, which depends sensitively on the slit width and is located far above the critical point of a hypothetical system with additive interactions (Fig. 5). This means that determining the phase diagram of a confined colloidal system provides an indirect but convenient method for probing experimentally the presence of non-additivity of critical Casimir interactions.
An important consequence of confinement and non-additivity is the formation and breaking of capillary bridges between colloids immersed into a solvent at its critical composition. Mean field theory predicts that this process can proceed via a first-order bridging transition, characterized by a remarkably extended hysteresis loop (Fig. 6). Experimentally, the formation of bridges has been observed, e.g., for glass colloids in a water/lutidine mixture close to its two-phase coexistence line. 74 However, experiments on bridging transitions are scarce. Our study proposes a convenient setup, with the slit width as a new control parameter, for analyzing such transitions and thereby providing a framework to observe them experimentally.
We have focused here only on a few aspects of the behaviour of confined colloidal systems. Such systems, however, can show a far richer behaviour. For instance, we have discussed colloidal liquid-gas transitions (Fig. 5), but two-dimensional crystalline structures can also form; 36,43 they may have potential applications in two-dimensional photonic sensors. 44 We have considered the temperature-dependent formation of bridges between confined colloids, but the slit width presents an additional parameter to control the bridges and the bridging transitions. Colloids can also help to detect interface localization/delocalization transitions, predicted for a noncritical fluid in the ordered phase confined between two walls with opposing boundary conditions. Although extensively studied theoretically, 75-80 such transitions have not yet been observed experimentally. Our preliminary calculations show that the force acting on a colloid confined in such a system experiences a vivid change when the interface localizes at one of the walls. Therefore this change can serve as an experimental indicator of the localization/delocalization transition.
We thus conclude that confined colloidal systems offer exciting opportunities to study a number of interesting phenomena and have the potential in technological applications. We expect that our results will stimulate future activities in this area, both in basic and applied research.

Mean field approach
The MFT equilibrium order-parameter profiles have been obtained by minimizing numerically the standard dimensionless Ginzburg-Landau-Wilson (LGW) Hamiltonian in spatial dimension d = 3 using the finite elements method with F3DM. 54 Here dx a is the d-dimensional volume element, j(r) is the order parameter (for a binary liquid mixture it is the deviation of the concentration of one species from its critical value), h is the external bulk field conjugate to j (in the present study mostly taken to be zero, see below), V is the volume accessible to the critical fluid, and the coupling constant g 4 0 stabilizes H below T c . Within MFT t = t/(x + 0 ) 2 where t = (T À T c )/T c is the reduced temperature, and x + 0 is the amplitude of the bulk correlation length in the disordered phase. Note that we consider here threedimensional spherical colloids, which form hyper-cylinders in d = 4 for which MFT predicts the correct critical behaviour.
The computational box for the numerical minimization was chosen sufficiently large in the lateral directions in order to ensure that the order parameter (OP) far from the colloids coincides with the OP profile dictated by the slit walls. A surface-to-surface distance of 5R from the colloids of radius R to the sides of the computational box turned out to be sufficient to meet this criterion. The values of the order parameter at these box sides were left free during minimization and there was no surface free energy associated with them due to choosing von Neumann boundary conditions in both lateral directions.
We consider the colloids and the slit walls to belong to the so-called extraordinary (or normal) surface universality class. Accordingly, the order parameter diverges upon approaching such surfaces (see, e.g., ref. 49, 50 and 81). In order to deal with this numerical challenge, we have used a short-distance expansion (see ref. 53 and 82) for calculating the value of the OP at some small distances d from the surface (we took d = 0.05R in all calculations). During minimization we kept the OP fixed on all mesh nodes belonging to that surface. This approach has proven to be successful in a number of previous studies. 19,[51][52][53]83,84 The critical Casimir force was calculated using the stress tensor method. 53 Each colloid was enclosed by an ellipsoidal surface (S i ) and the force acting on colloid i was calculated by integrating the stress tensor over S i (Section S1 in the ESI †).
For calculating the curve for O = N in Fig. 3d, we have used the 'À' boundary condition at the surfaces of both colloids and a weak opposing external bulk field h 4 0 (see eqn (3)), in order to mimic the effect of having '+' walls far from the colloids. We checked that the calculated forces are practically indistinguish- is the correlation length along the critical isotherm within mean field theory. Apart from this aspect we have set the bulk field to zero in all other calculations. (For further details concerning x h see eqn (4.2) in ref. 62 and Section 4.7 therein, as well as the appendix and the text after eqn (21) in ref. 85.) The hysteresis loop in Fig. 6b has been determined as follows. We first minimized the LGW Hamiltonian for a certain small t just below t = 0, starting from a random OP configuration. Then we used this solution as a starting configuration in a new minimization process for a lower, but still sufficiently small temperature t, as to ensure that the minimizer chooses the required (local) minimum. We continued this procedure until the lowest temperature presented in Fig. 6b was reached. The return path of the loop was obtained similarly by starting from low t o 0, at which there is no metastability; then we increased t stepwise and proceeded as described above.

Monte Carlo simulations
We have considered the standard Ising model with classical spins AE1 located on a cubic lattice, which mimics an incompressible binary mixture with spins corresponding to different species of the mixture. This spin model is a representative of the d = 3 Ising universality class, pertinent to the vicinity of the consolute critical point of the binary liquid mixture. For this model the value of the inverse critical temperature is b c = 1/k B T c = 0.2216544(3) and amplitudes of the correlation length for the low-and hightemperature regions are x À 0 = 0.243(1) and x + 0 = 0.501(2), respectively, in units of the lattice constant. 86 The critical exponent of the correlation length is n = 0.63002(10). 63 The system is limited by two confining surfaces along the z coordinate (Fig. 1) and is periodic in the x and y directions. We have considered the so-called normal surface universality class, 49,50,81 which corresponds to an infinitely strong surface field applied to the surface sites of the lattice (h a 1 = AEN, a A {top,bot} for the upper or lower wall). This implies that effectively the surface spins are frozen (i.e., they do not fluctuate) and their values are set to +1 or À1, depending on the sign of h a 1 . The colloids are defined by drawing a sphere of radius R with its center located at a lattice site, and all spins within this sphere are treated as being frozen (ref. 64 and Section S2 in the ESI †). Similar to the slit walls, we consider colloids the surfaces of which belong to the normal surface universality class, and hence we also fixed the spins on the lattice sites neighboring the sites of the colloid surfaces.
The critical Casimir potential was calculated by computing the difference between the insertion free energies of the system with two colloids (placed along the x direction, see Fig. 1) separated by a distance D and a system with two colloids a distance D max apart, where D max is a sufficiently large separation at which two colloids de facto do not interact with each other (Section S2 in the ESI †). In most of our calculations the system size in terms of lattice sites was N x Â N y = 100 Â 40 large in the x and y directions; the size in the z direction is determined by the slit width. We used the hybrid MC method, 87 in which each MC step consists of a combination of Wolff cluster update and N x Â N y Â N z /2 single spin Metropolis updates. For thermalization, 1-2 Â 10 5 MC steps were performed and thermal averaging was obtained from more than 1-2 Â 10 6 hybrid MC steps, which have been split into 10 series for the evaluation of the statistical errors. In order to check the presence of hysteresis loops we proceeded in the same way as for the MFT analysis. In these simulations, which probe metastable states, we did not use the cluster moves in order to prevent the system from jumping into the global minimum.

Calculation of the colloidal phase diagram
We have considered the case that the colloids position themselves in the mid-plane of a slit, as shown schematically in Fig. 1. For thin slits this is a minor restriction (Fig. S6 in the ESI †). Accordingly, we studied a two-dimensional system of colloids with its grand canonical thermodynamic potential J ¼ U À TS À mN as the sum of the contributions due to the excluded volume and the long-ranged critical Casimir interactions; here m is the chemical potential and N = rA is the number of colloids (r is the twodimensional density of colloids and A the surface area). Using the random phase approximation one readily obtains U % U 0 Z 2 ð2pÞ, where Z = prs 2 is the two-dimensional colloidal packing fraction, U 0 ¼ 2p s 2 À ÁÐ 1 2s FðrÞrdr is an interaction parameter, and s is the hard-disk radius. Considering charged colloids, we included the (short-ranged) screened electrostatic repulsion into the hard-disk reference system so that s 4 R, where R is the geometric radius of the colloidal particles. We used Rowlinson's description 88 of the hard sphere system in order to determine s. For the entropy S we used the analytical expression for hard discs following from scaled-particle theory. 89,90 In order to obtain an equilibrium packing fraction for a given chemical potential, we minimized J ðZÞ numerically. The colloidal phase diagram was constructed by identifying the region in which J has two minima corresponding to high and low packing fractions. For details of the calculations see Section S3 in the ESI. †

Conflicts of interest
There are no conflicts to declare.