text:
# Critical Casimir Forces in Colloidal Suspensions

## I Introduction

In 1978 Fisher and de Gennes  predicted that the confinement of critical fluctuations of the order parameter in a binary liquid mixture near its critical demixing point $`T_c`$ gives rise to long-ranged forces between immersed plates or particles. In particular, they pointed out that these long-ranged forces would eventually lead to the flocculation of colloidal particles which are dissolved in a near-critical binary liquid mixture . A few years later, such a solvent-mediated flocculation was observed experimentally for silica spheres immersed in the binary liquid mixture of water and 2,6-lutidine . However, the precise interpretation of these experiments is still under debate; although fluctuation-induced forces as predicted by Fisher and de Gennes certainly play a major role, additional mechanisms, such as screening effects in the case of charged colloidal particles , may also contribute . In the meantime additional experimental evidence for this kind of flocculation phenomena has emerged for other binary mixtures acting as solvents .

From a colloid physics perspective effective, solvent-mediated interactions between dissolved colloidal particles are of basic importance . The richness of the physical properties of these systems is mainly based on the possibility to tune these effective interactions over wide ranges of strength and form of the interaction potential. Traditionally, this tuning is accomplished by changing the chemical composition of the solvent by adding salt, polymers, or other components . Compared with such modifications, changes of the temperature or pressure typically result only in minor changes of the effective interactions. On the other hand, however, effective interactions generated by bringing solvents close to a phase transition of their own are extremely sensitive to such changes. In this broader context, flocculation of colloidal particles induced by critical fluctuations near a second-order phase transition of the solvent as predicted by Fisher and de Gennes  provided the first instance in which effective interactions are tuned by changing the physical conditions rather than the chemical composition of the solvent. Later on, this idea has been extended to first-order transitions of the solvent, for which wetting transitions can occur which give rise to wetting films of the preferred phase coating the colloidal particles . In the meantime the study of effective interactions between dissolved particles generated by changing the physical conditions of the solvent, and the resulting phase behavior of the colloidal suspension, has emerged as one of the main fields of research in colloid physics .

The physical origin of the force predicted by Fisher and de Gennes is analogous to the Casimir force between conducting plates, which arises due to the confinement of quantum mechanical vacuum fluctuations of the electromagnetic field . (Related phenomena occur in string theory , liquid crystals (see, e.g., ref. ), and microemulsions . However, in contrast to the present so-called “critical Casimir effect”, in liquid crystals the director-fluctuation induced forces are notoriously difficult to separate from the background of dispersion forces due to the absence of a singular temperature dependence .) In recent years this “critical Casimir effect” has attracted increasing theoretical  and experimental  interest. In this case the fluctuations are provided by the critical order parameter fluctuations of a bona fide second-order phase transition in the corresponding bulk sample. For He4 wetting films close to the superfluid phase transition the theoretical predictions  for critical Casimir forces between parallel surfaces exhibiting Dirichlet boundary conditions have been confirmed quantitatively . These quantum fluids offer the opportunity to study critical Casimir forces in the absence of surface fields, i.e., the purely fluctuation induced forces.

In a classical binary liquid mixture near its critical demixing point, the order parameter is a suitable concentration difference between the two species forming the liquid. In this case the inevitable preference of confining boundaries for one of the two species results in the presence of effective surface fields leading to nonvanishing order parameter profiles even at $`TT_c`$ . (This is the case discussed by Fisher and de Gennes .) These so-called “critical adsorption” profiles become particularly long-ranged due to the correlation effects induced by the critical fluctuations of the order parameter of the solvent. In the case of a planar wall, critical adsorption has been studied in much detail . Asymptotically close to the critical point $`T_c`$ it is characterized by a surface field which is infinitely large so that the order parameter profile actually diverges upon approaching the wall up to atomic distances. Later on, critical adsorption on curved surfaces of single spherical and rodlike colloidal particles has been studied, where it turns out, in particular, that critical adsorption on a microscopically thin “needle” represents a distinct universality class of its own . Critical adsorption on the rough interface between the critical fluid and its noncritical vapor occurring at the critical end point of the binary liquid mixture has been addressed as well ; likewise the effect of quenched surface modulations and roughness . The interference of critical adsorption on neighboring spheres gives rise to the critical Casimir forces which have been argued to contribute to the occurrence of flocculation near $`T_c`$ . However, a quantitative understanding of these phenomena requires the knowledge of the critical adsorption profiles near the colloidal particles and the resulting effective free energy of interaction in the whole vicinity of the critical point, i.e., as function of both the reduced temperature $`t=(TT_c)/T_c`$ and the bulk field $`h`$ conjugate to the order parameter. So far, this ambitious goal has not yet been accomplished. Instead, the introduction of the additional complication of a surface curvature has limited the knowledge of the corresponding critical adsorption so far to the case of spheres for the particular thermodynamic state $`(t,h)=(0,0)`$ of the solvent .

In a previous work , at least the temperature dependence of the critical Casimir force between a sphere and a planar container wall has been discussed. It turns out that the force becomes maximal not at but above $`T_c`$. In ref.  also the force $`F`$ between a pair of spheres of radius $`R`$ at distance $`D+2R`$ between their centers (see Fig. 1) has been briefly analyzed. Close to the critical demixing point, this force separates into a regular background contribution and a singular contribution $`F_{\text{sing}}^\pm `$ of universal character, which is attractive if the same of the two coexisting bulk phases is enriched near both spheres. The results for $`F_{\text{sing}}^\pm `$ can be cast into the form

$$F_{\text{sing}}^\pm (D,R;t,h)=\frac{k_BT_c}{R}K_\pm \left(\mathrm{\Delta }=\frac{D}{R},\mathrm{\Theta }_\pm =\frac{D}{\xi _\pm },\mathrm{\Xi }=\text{sgn}(h)\frac{D}{l_H}\right),$$
(1)
where $`\xi =\xi _\pm `$ is the true correlation length governing the exponential decay of the order parameter correlation function in the bulk for $`t=(TT_c)/T_c0`$ at $`h=0`$; $`l_H`$ is the correlation length governing the exponential decay of the order parameter correlation function for $`h0`$ at $`t=0`$. $`K_\pm `$ are universal scaling functions. Note that surface fields $`h_1`$ break the symmetry $`hh`$, i.e., $`F_{\text{sing}}^\pm `$ depends on the sign of $`\mathrm{\Xi }`$. Here and in the following we consider only the fixed point value $`h_1+\mathrm{}`$.

The aim of this work is to study the critical Casimir force between two spherical colloidal particles as function of both relevant thermodynamic variables $`t`$ and $`h`$. Figure 2 shows the schematic phase diagram of the solvent, e.g., a binary liquid mixture. Different thermodynamic paths are indicated along which the dependence of the force on $`t`$ and $`h`$ is studied.

The remainder of this work is organized as follows: In Sec. II we define the field-theoretic model used to describe the near-critical solvent. In Sec. III the numerical calculation of the adsorption profiles in mean-field approximation is discussed. In Sec. IV we present the results for the force between two spherical particles and compare them with the limiting cases of small and large particle radii. The discussion of the results for two parallel plates, which inter alia are relevant for a pair of spheres via the Derjaguin approximation, is left to Sec. V. Finally, in Sec. VI, we summarize and discuss our results with regards to the phase behavior of a suspension of colloidal particles in a near-critical fluid.

## II Field-theoretic model

Near criticality the behavior of the system is governed by fluctuations of the order parameter $`\mathrm{\Phi }`$ on a large length scale such that the emerging universal properties are independent of microscopic details. A coarse-graining procedure leads to the continuum description in terms of the Ginzburg-Landau Hamiltonian

$$\{\mathrm{\Phi }(𝐫)\}=\underset{V}{}d^dr\left\{\frac{1}{2}[\mathrm{\Phi }(𝐫)]^2+\frac{\tau }{2}\mathrm{\Phi }^2(𝐫)+\frac{u}{24}\mathrm{\Phi }^4(𝐫)h\mathrm{\Phi }(𝐫)\right\}$$
(2)
of an order parameter field $`\mathrm{\Phi }(𝐫)`$ in the bulk of volume $`V`$. For the liquids considered here $`\mathrm{\Phi }(𝐫)`$ is a scalar field. The Hamiltonian determines the statistical weight $`\mathrm{exp}[\{\mathrm{\Phi }(𝐫)\}]`$ of the configuration $`\{\mathrm{\Phi }(𝐫)\}`$. The external bulk field $`h(𝐫)`$ is conjugate to the order parameter and breaks its symmetry. The temperature dependence enters via

$$\tau t=\frac{TT_c}{T_c}.$$
(3)
The partition function of the system is given by

$$𝒵(\tau ,h;u)=𝒟\{\mathrm{\Phi }\}\mathrm{exp}[\{\mathrm{\Phi }(𝐫)\}],$$
(4)
where $`𝒟\{\mathrm{\Phi }\}`$ denotes an appropiately defined functional integration . The corresponding singular part of the free energy reads

$$_{\mathrm{sing}}(\tau ,h;u)=k_BT_c\mathrm{ln}𝒵(\tau ,h;u).$$
(5)
The coarse-graining procedure leading to a continuum model from a given lattice model can be generalized to systems with surfaces resulting in the Hamiltonian

$`\{\mathrm{\Phi }(𝐫)\}`$ $`=`$ $`{\displaystyle \underset{V}{}}d^dr\left\{{\displaystyle \frac{1}{2}}[\mathrm{\Phi }(𝐫)]^2+{\displaystyle \frac{\tau }{2}}\mathrm{\Phi }^2(𝐫)+{\displaystyle \frac{u}{24}}\mathrm{\Phi }^4(𝐫)h\mathrm{\Phi }(𝐫)\right\}`$ (6)
$`+`$ $`{\displaystyle \underset{S}{}}𝑑S\left\{{\displaystyle \frac{c}{2}}\mathrm{\Phi }^2(𝐫_S)h_1\mathrm{\Phi }(𝐫_S)\right\}.`$ (7)
The first integral runs over the volume $`V`$ available for the critical medium, which is bounded by the surface(s) $`S=V`$. Apart from the restricted space integration, the first integral in Eq. (7) has the same form as the Hamiltonian (2) for bulk systems disregarding surfaces. The second integral in Eq. (7) represents a surface contribution in which the integral runs over the boundary surface $`S=V`$; $`c`$ is related to the strength of the coupling of critical degrees of freedom near the surface and $`h_1`$ is the surface analogue of the bulk field $`h`$ .

The systematic treatment of the critical fluctuations as implied by Eq. (4) and a justification of the form of $`\{\mathrm{\Phi }(𝐫)\}`$ is provided by the field-theoretic renormalization group framework . The basic idea is to consider the physical system on larger and larger length scales and to single out those quantities which eventually remain unchanged under the repeated action of such scale transformations. Such quantities tend to their fixed point values right at the critical point . For bulk systems, there remain only two relevant adjustable parameters, corresponding to $`\tau `$ and $`h`$. The surface generates a subdivision of a given bulk universality class into different surface universality classes characterized by different fixed point values for $`c`$ and $`h_1`$. For $`h_1=0`$ the corresponding leading singular behaviors are determined by the fixed-point values $`c=+\mathrm{}`$ for the ordinary transition (O), $`c=0`$ for the special or surface-bulk transition, and $`c=\mathrm{}`$ for the extraordinary or normal transition (E) . For the O and E transitions only $`\tau `$, $`h`$, and $`h_1`$ remain as relevant variables. The three cases above were originally studied for magnetic systems. For liquids, however, E is the generic case and therefore it is called the normal transition. The case $`h_1=\mathrm{}`$ with $`c`$ arbitrary corresponds to critical adsorption, but the leading singular behavior is the same as for the E transition .

Neglecting fluctuations the Ginzburg-Landau Hamiltonian $`\{\mathrm{\Phi }(𝐫)\}`$ itself represents a free energy functional which upon minimization with respect to $`\mathrm{\Phi }(𝐫)`$ yields the Euler-Lagrange equations for order parameter profiles and correlation functions. The minimization is equivalent to the mean-field description of the system,

$$\frac{\delta \{\mathrm{\Phi }(𝐫)\}}{\delta \mathrm{\Phi }(𝐫)}|_{\mathrm{\Phi }=\mathrm{\Phi }}=0,$$
(8)
where the fluctuating order parameter $`\mathrm{\Phi }`$ is replaced by its mean value $`\mathrm{\Phi }`$. The universal critical behavior is captured correctly by mean-field theory if the space dimension $`d`$ exceeds the upper critical dimension $`d_{\text{uc}}=4`$.

Because of the relative complexity of the geometry considered here (see Fig. 1), it is not possible to carry out the complete renormalization procedure in closed form. Instead, available knowledge for bulk systems in $`d=3`$ is used in combination with the mean-field solution for the geometry considered here in order to estimate universal quantities in $`d=3`$. An alternative way to overcome mean-field theory is followed within the concept of the Local Functional Theory developed in refs. , in which the free energy is described by a functional depending locally on the order parameter and its derivative. With mean-field theory as a special case, it takes into account available knowledge in $`d=3`$, such as the correct values for the critical exponents.

## III Order parameter profiles

In this Section we describe some technical aspects of the tools and numerical methods used for obtaining the mean-field results, i.e., the determination of the order parameter profiles at the critical point, the short distance expansion of the profile close to a surface, the minimization procedure to determine the profiles away from the critical point, and the calculation of the force between the particles based on the stress tensor.

Minimization in Eq. (8) leads to the differential equation

$$\mathrm{\Delta }m+\tau m+m^3H=0,$$
(9)
where the coupling constant $`u`$ is absorbed in the order parameter and in the bulk magnetic field,

$$m=\sqrt{\frac{u}{6}}\mathrm{\Phi },H=\sqrt{\frac{u}{6}}h.$$
(10)
Note that we are interested in an estimate of the order parameter profile $`m`$ in $`d=3`$ where, at the fixed point, $`u`$ is a positive number. For $`d4`$ the mean-field solution is exact, but the fixed point value of the coupling constant $`u`$ vanishes in this limit.

The surface integrals in Eq. (7) lead to a boundary condition for the order parameter at the surfaces of the spherical particles. At the critical adsorption fixed point this boundary condition is given by

$$\mathrm{\Phi }_S=\mathrm{}.$$
(11)
### A Order parameter profiles at criticality

While our main interest is the dependence of the Casimir force on $`t`$ and $`h`$, the knowledge of the order parameter profile at the critical point, i.e., for $`(t,h)=(0,0)`$, is a useful starting point. In ref.  this latter profile is given for a critical system confined to a spherical shell between two concentric spheres with radii $`R_\pm `$ at which the order parameter diverges. In terms of elliptic functions it reads

$$m(r)=\frac{\sqrt{{\displaystyle \frac{1k^2}{k^2\frac{1}{2}}}}}{r\text{cn}\left(\sqrt{{\displaystyle \frac{1}{2k^21}}}\mathrm{ln}\left({\displaystyle \frac{r}{R_0}}\right)\right)},$$
(12)
where $`r`$ is the radial distance from the common center of the spheres and $`k`$ is the module of the elliptic function cn . The module $`k`$ must be adjusted such that the profile diverges at the given radii $`R_\pm `$ leading to the following implicit equation for $`k`$:

$$\frac{R_\pm }{R_0}=\left(\frac{R_+}{R_{}}\right)^{\pm \frac{1}{2}}=\mathrm{exp}\left(\pm \sqrt{2k^21}K(k)\right),$$
(13)
where $`R_0=\sqrt{R_+R_{}}`$ is the geometric mean of the radii and $`K`$ is the complete elliptic integral of the first kind.

In ref.  it is outlined how the results for concentric spheres can be used for different geometries via a conformal transformation of the coordinates,

$$\frac{𝐫^{}}{|𝐫^{}|^2}=\frac{𝐫+𝐑}{|𝐫+𝐑|^2}\frac{𝐑}{2|𝐑|^2},$$
(14)
where $`𝐫^{}`$ and $`𝐫`$ are the new and the old coordinates, respectively, and $`𝐑`$ defines the symmetry axis and the eventual shape of the new geometry. For $`|𝐑|=R_0`$ the two concentric spheres with radii $`R_\pm `$ are mapped onto two spheres of equal size which is the generic case we are interested in. For $`|𝐑|=R_+`$ the two concentric spheres are mapped onto a single spherical particle in front of a planar wall for which results were obtained in ref. .

Under the conformal transformation (14), each local scaling field is multiplied by a scale factor

$$b(𝐫^{})=1+\frac{𝐑𝐫^{}}{|𝐑|^2}+\frac{|𝐫^{}|^2}{4|𝐑|^2},$$
(15)
which is raised to the power of the scaling dimension of this scaling field. In mean-field approximation and $`d=4`$, the scaling dimension of the order parameter $`\mathrm{\Phi }`$ is $`1`$, which implies for the order parameter profile

$$m(𝐫^{})=b(𝐫^{})^1m_{\text{conc}}(|𝐫(𝐫^{})|),$$
(16)
where $`m_{\text{conc}}(|𝐫|)`$ is the profile in the concentric geometry given by Eq. (12).

Let $`D`$ denote the surface-to-surface distance between the two spheres and $`R`$ their common radius; then the corresponding radii $`R_\pm `$ in the concentric geometry are given by

$$R_\pm =\frac{D}{4}A\left(\frac{A+1}{A1}\right)^{\pm 1}$$
(17)
with the abbreviation

$$A=\sqrt{1+4\frac{R}{D}}.$$
(18)
Putting the pieces together, Eqs. (10)-(17) yield the order parameter profile in the geometry of two spherical particles with radius $`R`$ a distance $`D`$ apart from each other right at the critical point $`(t,h)=(0,0)`$.

### B Short distance expansion

In the following we describe the methods for obtaining results off criticality, i.e., for $`t0`$ and $`h0`$, i.e., $`\tau 0`$ and $`H0`$. Equation (9) will then yield a different profile $`m(𝐫^{})`$. However, the profile at the critical point serves as a suitable starting point for the minimization procedure (8). Close to the surfaces it is even possible to calculate the deviation of the order parameter from the profile at the critical point analytically. In this context the so-called short distance expansion turns out to be useful, and has been carried out for $`\tau 0`$ in ref. . With $`H0`$ in addition, it reads:

$`m(s0;R,\tau ,H)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{s}}{\displaystyle \frac{\sqrt{2}}{6}}{\displaystyle \frac{d1}{R}}`$ (19)
$`+`$ $`{\displaystyle \frac{\sqrt{2}}{6}}s\left\{\tau +\left[{\displaystyle \frac{5}{6}}{\displaystyle \frac{(d1)^2}{R^2}}{\displaystyle \frac{(d1)(d2)}{R^2}}\right]\right\}`$ (20)
$`+`$ $`{\displaystyle \frac{\sqrt{2}}{4}}s^2\left({\displaystyle \frac{H}{\sqrt{2}}}{\displaystyle \frac{d1}{3R}}\tau {\displaystyle \frac{4d^339d^2+120d85}{54R^3}}\right)`$ (21)
$`+`$ $`𝒪(s^3),`$ (22)
where $`s`$ is the radial distance from the surface of a sphere with radius $`R`$. The difference $`\delta m=m(s0;R,\tau ,H)m(s0;R,\tau =0,H=0)`$ gives the deviation from the profile at the critical point close to the surfaces. While the profile itself diverges at the surfaces, the deviation $`\delta m`$ vanishes with the radial distance $`s`$ as

$$\delta m(s0;R,\tau ,H)=\frac{\sqrt{2}}{6}\tau s+\frac{1}{4}\left(H\frac{\sqrt{2}}{3}\frac{d1}{R}\tau \right)s^2+\mathrm{},$$
(23)
where the bulk field $`H`$ and the radius $`R`$ of the spherical particles appear only in $`𝒪(s^2)`$ and higher order. Equation (19) is only valid for small $`s`$ since the presence of the second particle is neglected. Corrections to Eqs. (19) and (23), however, are of order $`𝒪(s^3)`$ as derived for the case of parallel plates, where they are called distant wall corrections .

### C Numerical solution for the complete profiles off criticality

So far we have provided results for the order parameter profile at the critical point for the two sphere geometry \[Eq. (16)\] and off criticality but close to the surface of a single sphere \[Eq. (19)\]. In order to obtain the complete order parameter profile away from the critical point for the geometry consisting of two spheres with equal size (see Fig. 1) additional, numerical effort is required.

Upon discretization the order parameter profile $`m(𝐫^{};D,R;\tau ,H)`$ in the two sphere geometry (Fig. 1) is described by a set of parameter values $`a_{ij}`$ where $`i`$ and $`j`$ indicate the dependence on the two independent spatial variables required for the present cylindrical symmetry around the axis between the two centers of the spheres \[see, c.f., Eq. (26)\]:

$$m(𝐫^{};D,R;\tau ,H)=m(\{a_{ij}\}).$$
(24)
Within mean-field theory $`(\{a_{i,j}\})`$ has to be minimized with respect to $`\{a_{ij}\}`$ for given $`D`$, $`R`$, $`\tau `$, and $`H`$. To this end we apply the method of steepest descent ,

$$a_{ij}^{(n+1)}=a_{ij}^{(n)}\kappa \frac{(\{a_{ij}\})}{a_{ij}}|_{a_{ij}^{\left(n\right)}},$$
(25)
where $`\kappa `$ is a positive number and $`a_{ij}^{(n)}`$ is the $`n`$-th iteration of the parameter $`a_{ij}`$.

It is known that in general the method of steepest descent converges slowly. Nevertheless in the present case it is found to be useful since the dependence on the choice of the initial set of parameters $`\{a_{ij}^{(0)}\}`$ is weak. In fact, it turns out that the profile at the critical point is a rather useful starting profile. The convergence also depends on the value of $`\kappa `$: if it is chosen too large, the convergence breaks down, if it is chosen too small, the convergence becomes too slow. As we will discuss later, one advantage of the present minimization method is that the explicit evaluation of the Hamiltonian $``$ is not necessary, since only the gradient $`\frac{(\{a_{ij}\})}{a_{ij}}`$ is needed, which can be obtained independently. In order to keep this advantage, iterative methods to determine the best choice for the value of $`\kappa `$, which explicitly need the value of the functional that is to be minimized, are not applied. Nonetheless the best choice of $`\kappa `$ is the largest value for which the minimization still proceeds.

In particular, we use the parametrization

$$m(𝐫_{ij}^{};D,R;\tau ,H)=m(𝐫_{ij}^{};D,R;0,0)+a_{ij},$$
(26)
where $`𝐫_{ij}^{}`$ denotes discretized spatial positions in the two sphere geometry, $`m(𝐫_{ij}^{};D,R;0,0)`$ is the value of the order parameter profile for $`(\tau ,H)=(0,0)`$ at those positions, according to Eq. (16), and $`a_{ij}`$ denotes the deviation from this value at the position $`𝐫_{ij}^{}`$ used in the minimization procedure (25). Thus starting the procedure with the profile for $`(\tau ,H)=(0,0)`$ implies the choice

$$a_{ij}^{(0)}=0,i,j.$$
(27)
In addition, Eq. (23) provides the behavior of the $`a_{ij}`$ for $`𝐫_{ij}^{}`$ close to the surface, i.e., $`a_{ij}^{(n)}=0`$ for all $`n`$ if $`𝐫_{ij}^{}=𝐫_S^{}`$. From Eq. (23) even the slope of a function $`f(s)`$ used for interpolating the $`a_{ij}`$ as function of the radial distance $`s`$ from the surface of the sphere is known for $`s0`$.

The cylindrical symmetry of the problem allows one to consider the planar geometry depicted in Fig. 1. Instead of discretizing the space by a square lattice we use a lattice which takes into account that there are areas with strong changes in the order parameter $`m`$, i.e., close to the surfaces, and areas with weak changes, i.e., far away from the spheres. In those areas with weak changes less lattice points are needed. For this purpose we use what we call the conformal lattice. The idea is to use the conformal mapping (14) in order to create this lattice from that one for the case of concentric spheres. In analogy to electrostatics the field lines and the lines of equal potential of a system with two concentric spheres forming the two electrodes of a capacitor are mapped onto the two sphere geometry such that the new lines are the field lines and the lines of equal potential of a system with two charged spheres. The intersections of these lines form the lattice sites of the conformal lattice. The advantage of this lattice is that the field lines approach the spheres perpendicularly to their surface, which allows one to use the short distance behavior \[Eqs. (19) and (23)\] as a boundary condition for the $`a_{ij}`$. Figure 3 shows the conformal lattice for the geometry of two spheres. In principle, it would be possible to use only one half of this lattice, since this system is symmetric with respect to the vertical midplane. However, this reduction leaves one with a disadvantage. The value of the order parameter at this midplane wall is not known, and thus it must be determined within the minimization process. The minimization turns out to be inefficient then. Therefore for the two sphere geometry we resort to the full conformal lattice enclosing both spheres as shown in Fig. 3.

The Hamiltonian $``$ \[Eq. (7)\] is given by

$$=\frac{6}{u}\underset{V}{}𝑑V(m)+\frac{6}{u}\underset{S}{}𝑑S_S(m),$$
(28)
with

$$(m)=\frac{1}{2}(m)^2+\frac{\tau }{2}m^2+\frac{1}{4}m^4Hm$$
(29)
and $`_S(m)=\frac{c}{2}m^2H_1m`$ with $`H_1=h_1\sqrt{u/6}`$. Performing the functional derivative with respect to the order parameter $`m`$ leads, via partial integration, to

$$\frac{\delta \{m(𝐫)\}}{\delta m(𝐫)}=\frac{6}{u}_V𝑑V\{\mathrm{\Delta }m+\tau m+m^3H\}+\frac{6}{u}_S𝑑S\{_nm+_S^{}(m)\},$$
(30)
where $`_nm`$ is the derivative of the order parameter in direction of the surface normal. We note that the surface integral vanishes, if the boundary conditions \[Eqs. (19) and (23)\] are fulfilled. Equation (25) contains the derivative of the Hamiltonian $``$ with respect to the parameters $`a_{ij}`$ which represent the order parameter $`m`$. With $`m=m(a_{ij})=const`$ in a small volume $`\delta V(𝐫_{ij}^{})`$ around the position $`𝐫_{ij}^{}`$ the gradient in Eq. (25) is approximated by

$$\frac{(\{a_{ij}\})}{a_{ij}}=\frac{6}{u}\delta V(𝐫^{})\{\mathrm{\Delta }m+\tau m+m^3H\}|_{𝐫^{}=𝐫_{ij}^{}}.$$
(31)
Note that a constant prefactor in Eq. (31) can be absorbed by a rescaling of $`\kappa `$ in Eq. (25). If Eq. (9) holds, i.e., if the expression in curly brackets in Eq. (31) vanishes, Eq. (25) reduces to

$$a_{ij}^{(n+1)}=a_{ij}^{(n)},$$
(32)
i.e., the minimization has led to a fixed point solution of the iteration.

For the minimization process the second derivatives of the order parameter profile are needed \[see Eq. (31)\]. Since the order parameter is calculated only at the discrete sites $`𝐫_{ij}^{}`$, numerical differentiation must be applied. For a numerical derivative, in general, it is preferable to have the sites on an equidistant square lattice. However, it is also possible to deal with derivatives on the conformal lattice. To this end we have used the spline approximation  in order to interpolate between all the $`a_{ij}`$ on each of the lines of the conformal lattice (see Fig. 3). This approximation directly yields the second derivative of the order parameter at the lattice sites with respect to the variables parametrizing the lines. Under a conformal mapping angles are conserved, so that the lines in the conformal lattice intersect perpendicularly which allows one to translate the derivatives obtained from the spline approximation into derivatives with respect to the cylindrical coordinates $`𝐫^{}=(\rho ,z)`$.

### D Stress tensor

With the tools now available to determine the order parameter profile in the two sphere geometry the next step is to calculate the effective forces between the two spheres. The straightforward way to determine such forces is to calculate the free energy $`k_BT\mathrm{ln}𝒵`$ as function of the distance $`D`$ between the spheres, and then to estimate the derivative with respect to the distance by considering finite differences. However, the total value of the free energy is expected to be large as compared with the differences, and thus one has to deal numerically with differences between large numbers, which causes high inaccuracies. Therefore we resort to an alternative method, in which the force is calculated directly using the so-called stress tensor $`T_{\mu \nu }`$, which is defined by the linear response

$$\delta =_Vd^dr\frac{b_\mu }{x_\nu }T_{\mu \nu }(𝐫),$$
(33)
where $`𝐛`$ is a non-conformal coordinate transformation, $`\mu `$ and $`\nu `$ index the spatial coordinates, and $`\delta `$ is the corresponding energy shift. With the Lagrangian $`(m)`$ in Eq. (29) it reads

$$T_{\mu \nu }=\frac{6}{u}\left[\frac{}{(_\nu m)}_\mu m\delta _{\mu \nu }\right].$$
(34)
The knowledge of the stress tensor on any infinitely extended surface which separates the two spheres from one another, or likewise on any closed, finite surface which contains only one of the spheres, allows one to calculate the force between the spheres. This is illustrated in terms of the spherical volume $`V_S`$ with radius $`\stackrel{~}{R}`$ containing the right sphere in Fig. 4. Using the coordinate transformation

$`b_z`$ $`=`$ $`\alpha ,𝐫V_S`$ (35)
$`=`$ $`0,𝐫V_S`$ (36)
$`b_{\nu z}`$ $`=`$ $`0,`$ (37)
the volume $`V_S`$ is shifted in $`z`$-direction by the amount $`\alpha `$. The derivative $`\frac{b_\mu }{x_\nu }`$ then leads to a delta function such that the integral over the whole volume $`V`$ reduces to an integral over the boundary of $`V_S`$. In this case the force $`F_{\text{sing}}=\frac{k_BT_c}{R}K`$ with the universal scaling function $`K`$ \[see Eq. (1)\] is given by

$`{\displaystyle \frac{\delta }{\alpha }}`$ $`=`$ $`{\displaystyle \frac{}{\alpha }}{\displaystyle _V}d^dr{\displaystyle \frac{b_\mu }{x_\nu }}T_{\mu \nu }`$ (38)
$`=`$ $`{\displaystyle _{V_S}}𝑑\mathrm{\Omega }n_\nu T_{z\nu }`$ (39)
$`=`$ $`4\pi \stackrel{~}{R}^3{\displaystyle _0^\pi }𝑑\phi \mathrm{sin}^2\phi [T_{zz}(\stackrel{~}{R},\phi )\mathrm{cos}\phi +T_{z\rho }(\stackrel{~}{R},\phi )\mathrm{sin}\phi ],`$ (40)
where $`𝐧`$ is the unit vector perpendicular to the boundary $`V_S`$. The quantities used in the parametrization in Eq. (39) are sketched in Fig. 4.

In the geometry with cylindrical symmetry \[$`m=m(\rho ,z)`$\] considered here, Eq. (34) leads to

$$T_{zz}=\frac{1}{u}\left[3\left(\frac{m}{z}\right)^2+4\frac{m}{\rho }\frac{m}{\rho }+2m\frac{^2m}{\rho ^2}\left(\frac{m}{\rho }\right)^23\tau m^2\frac{3}{2}m^4\right]$$
(41)
and

$$T_{z\rho }=\frac{1}{u}\left[4\left(\frac{m}{z}\right)\left(\frac{m}{\rho }\right)2m\frac{^2m}{\rho z}\right],$$
(42)
where the so-called improvement term $`\frac{1}{u}(_\mu _\nu m^2\delta _{\mu \nu }\mathrm{\Delta }m^2)`$ has been added .

In our calculation, we choose $`V_S`$ as indicated in Fig. 3, determine the stress tensor at the lattice points, and use spline interpolations between them. Note that the possibility of choosing different spherical surfaces enclosing one of the spheres, which via Eq. (39) should all result in the same force $`F_{\text{sing}}`$, provides a valuable and sensitive check for both the validity and the accuracy of the numerical results.

## IV Results for the two-particle interaction

The goal of this work is to provide a major step towards the understanding of aggregation of colloidal particles dissolved in a binary liquid mixture close to criticality. In the previous section we explained the various methods for obtaining the critical Casimir force between two spherical particles of equal size. The results are presented in this section.

The aggregation phenomenon is driven by the dependences on the temperature $`t=(TT_c)/T_c`$ and on the bulk composition $`c_A`$ of the solvent. While the former is related to the temperature parameter $`\tau `$ the latter is related to the bulk field $`H`$. Thus also the Casimir force between two particles must be determined as a function of those two parameters. We present results for the scaling function $`K_\pm (\mathrm{\Delta }=\frac{D}{R},\mathrm{\Theta }_\pm =\frac{D}{\xi _\pm },\mathrm{\Xi }=\text{sgn}(h)\frac{D}{l_H})`$ which is related to the force $`F_{\mathrm{sing}}^\pm `$ between the particles via Eq. (1). The correlation length $`\xi _\pm `$ for $`H=0`$ depends on the temperature $`T`$ as \[see also Eq. (3)\]

$$\xi _\pm =\xi _0^\pm |t|^\nu =\{\begin{array}{cc}\hfill \tau ^{1/2},& \tau >0\hfill \\ \hfill |2\tau |^{1/2},& \tau <0\hfill \end{array},H=0,$$
(43)
with $`\nu =\frac{1}{2}`$ within mean-field theory while the correlation length $`l_H`$ at $`T_c`$ depends on the bulk field $`H`$ as

$$l_H=l_0\left|\frac{c_Ac_{A,c}}{c_{A,c}}\right|^{\nu /\beta }=\frac{1}{\sqrt{3}}|H|^{1/3},\tau =0,$$
(44)
where $`\nu /\beta =1`$ within mean-field theory; $`\xi _0^\pm `$ and $`l_0`$ are nonuniversal amplitudes. The first parts of Eqs. (43) and (44) hold in general whereas the second parts correspond to the present mean-field description. The scaling variable $`\mathrm{\Xi }=\text{sgn}(H)D/l_H`$ depends on the sign of $`H`$ \[Eq. (1)\]. The relation between the bulk field $`H`$ and the concentration $`c_A`$ will be discussed further in Subsec. IV G, but we keep in mind that a negative value of $`H`$ means that the fluid component, which is preferred by the particles, is the one with a low concentration in the bulk (see Fig. 2).

### A Dependence on temperature at the critical composition

The variation of the temperature leads to an interesting behavior of the force . According to Fig. 5 the force exhibits a maximum for a temperature above the critical point. The position of the maximum of the force as function of the scaling variable $`\mathrm{\Theta }_+`$ depends on the distance $`D`$ between the particles. In the limit of large distances the position $`\mathrm{\Theta }_{\mathrm{max}}=\mathrm{\Theta }_{\mathrm{max}}(\mathrm{\Delta },\mathrm{\Xi }=0)`$ of the maximum approaches that one of the critical point, i.e.,

$$\mathrm{\Theta }_{\mathrm{max}}\left(\mathrm{\Delta }=\frac{D}{R}\mathrm{}\right)0.$$
(45)
On the other hand, for small distances the position of the maximum approaches a finite value above $`T_c`$

$$\mathrm{\Theta }_{\mathrm{max}}\left(\mathrm{\Delta }0\right)\mathrm{\Theta }_{\mathrm{max}}^{(0)}.$$
(46)
These two limiting behaviors are shown in Fig. 5. The disappearance of the maximum for large distances along with the derivation of the small radius expansion will be discussed in more detail later on in this section.

### B Dependence on composition at the critical temperature

The above results show that there are rich structures in the force curves off criticality. Therefore it is important to monitor the force not only at the critical point, but also in its vicinity. This becomes even more evident if not the temperature but the composition of the solvent is varied. Along the dashed path in Fig. 2 and for small distances $`D`$ the force exhibits a pronounced maximum as can be seen in Fig. 6. For large distances there is only a maximum at the critical point, i.e., $`H=0`$. For $`\mathrm{\Xi }>0`$ the force decays monotonously for all distances $`D`$.

Figure 7 shows the maximum values of the force for a given distance $`D`$ as function of the temperature (solid line in Fig. 2) or the composition (dashed line in Fig. 2). For comparison, also the force at the critical point for the same particle distance $`D`$ is plotted (dash-dotted line in Fig. 7). For intermediate and larger distances the three lines merge, which is due to the fact that for large distances the position of the maximum of the force approaches that of the critical point. For small distances the differences are pronounced. Obviously, the composition of the solvent plays an important role.

### C Dependence on temperature off the critical composition

So far the thermodynamic paths were chosen such that either the temperature or the solvent composition were fixed at their critical point value. However, it is experimentally relevant to consider also thermodynamic paths along which the temperature is varied at fixed compositions $`c_Ac_{A,c}`$. The corresponding results are given in Fig. 8 for a fixed distance $`D`$. For small deviations $`|c_Ac_{A,c}|c_{A,c}`$ the position of the maximum of the Casimir force is more or less unchanged while its absolute value is changed considerably. However, the main features of the temperature scan do not differ much from that at $`c_{A,c}`$ so that they are sufficiently robust in view of possible experimental tests with finite resolution for $`c_A`$.

If the dotted path in Fig. 2 is extended to temperatures below $`T_c`$ the force is analytic as function of the temperature until the path intersects with the coexistence line.

### D Derjaguin approximation

In the limit that the radius $`R`$ of the spherical particles is much larger than both the correlation length $`\xi `$ and the distance $`D`$ of closest approach surface-to-surface between the two particles, the particles can be regarded as composed of a pile of fringes. Each fringe builds a fringe-like slit with distance $`L=D+r_{||}^2/R`$, where $`r_{||}`$ is the radius of the fringe. Within this Derjaguin approximation  and in $`d=4`$, the force between the particles in units of $`\frac{k_BT_c}{R}`$ \[see Eq. (1)\] is given by

$$K_\pm (\mathrm{\Delta },\mathrm{\Theta }_\pm ,\mathrm{\Xi })=4\pi \mathrm{\Delta }^{5/2}_0^{\mathrm{}}𝑑\upsilon \upsilon ^2\frac{K_\pm ^{||}(\mathrm{\Theta }_\pm (1+\upsilon ^2),\mathrm{\Xi }(1+\upsilon ^2))}{(1+\upsilon ^2)^4},$$
(47)
where $`\upsilon `$ is a dimensionless integration variable and $`K_\pm ^{||}(L/\xi _\pm ,\text{sgn}(H)L/l_H)`$ is the analogous scaling function for the Casimir force in the slit geometry with parallel walls at distance $`L`$ \[see, c.f., Sec. V\].

### E Small radius expansion

If on the other hand the radius $`R`$ $``$ albeit large on the microscopic scale $``$ is much smaller than $`D`$ and $`\xi `$, the statistical Boltzmann weight $`e^{\delta _S}`$ characterizing the presence of the sphere centerd at $`𝐫_S`$ can be systematically expanded in terms of increasing powers of $`R`$ , i.e.,

$$e^{\delta _S}/d_{}=1+c_{}^\mathrm{\Phi }R^{x_\mathrm{\Phi }}\mathrm{\Phi }(𝐫_S)+c_{}^{\mathrm{\Phi }^2}R^{x_{\mathrm{\Phi }^2}}\mathrm{\Phi }^2(𝐫_S)+\mathrm{}$$
(48)
where $`d_{}`$ is an amplitude and where $`x_\mathrm{\Phi }=\beta /\nu `$ and $`x_{\mathrm{\Phi }^2}=d\nu ^1`$ (with the standard bulk critical exponents $`\beta `$ and $`\nu `$) are the scaling dimensions of $`\mathrm{\Psi }=\mathrm{\Phi },\mathrm{\Phi }^2`$. The ellipses stand for contributions which vanish more rapidly for $`R0`$. The coefficients $`c_{}^\mathrm{\Psi }`$ are fixed by $`c_{}^\mathrm{\Psi }=A_{}^\mathrm{\Psi }/B_\mathrm{\Psi }`$, where $`A_{}^\mathrm{\Psi }`$ and $`B_\mathrm{\Psi }`$, respectively, are amplitudes of the half-space (hs) profile $`\mathrm{\Psi }(z)_{\mathrm{hs},T=T_c}^{}=A_{}^\mathrm{\Psi }(2z)^{x_\mathrm{\Psi }}`$ at the critical point of the fluid for the boundary condition $``$ corresponding to critical adsorption, and of the bulk two-point correlation function $`\mathrm{\Psi }(𝐫)\mathrm{\Psi }(0)_{b,T=T_c}=B_\mathrm{\Psi }r^{2x_\mathrm{\Psi }}`$ . Since both spheres are equally small the free energy $``$ follows from the bulk average of two statistical weights of the type given in Eq. (48):

$$\frac{(D)(\mathrm{})}{k_BT_c}=\mathrm{ln}\left[\frac{e^{\delta _{S_1}}e^{\delta _{S_2}}_\mathrm{b}}{e^{\delta _{S_1}}_\mathrm{b}e^{\delta _{S_2}}_\mathrm{b}}\right],$$
(49)
where the spheres $`S_1`$ and $`S_2`$ are a distance $`D`$ apart from each other. Using Eq. (48) we find

$$\frac{(D)(\mathrm{})}{k_BT_c}=(c_{}^\mathrm{\Phi }R^{x_\mathrm{\Phi }})^2\mathrm{\Phi }(𝐫_{S_1})\mathrm{\Phi }(𝐫_{S_2})_\mathrm{b}[12c_{}^\mathrm{\Phi }R^{x_\mathrm{\Phi }}\mathrm{\Phi }_\mathrm{b}+\mathrm{}],$$
(50)
where $`\mathrm{\Phi }(𝐫_{S_1})\mathrm{\Phi }(𝐫_{S_2})_\mathrm{b}`$ is the bulk two point correlation function and $`\mathrm{\Phi }_\mathrm{b}`$ provides a measure of the deviation of the bulk composition from its value at the critical point. The Casimir force is the derivative of the free energy with respect to the distance $`D`$ between the spheres. Note that in Eq. (50) the first term in the square brackets gives rise to an attractive force which is increased by the second term if $`\mathrm{\Phi }_\mathrm{b}`$ is negative, i.e., if the binary liquid mixture is poor in the component preferred by the colloids (compare Subsecs. IV B and IV C).

Within mean-field approximation, the correlation function is given by (see Eq. (B5) in ref. )

$`\mathrm{\Phi }(𝐫_{S_1})\mathrm{\Phi }(𝐫_{S_2})_b`$ $`=`$ $`G_b^+(|𝐫_1𝐫_2|=D;\xi _+)`$ (51)
$`=`$ $`B_\mathrm{\Phi }{\displaystyle \frac{1}{D^{d2}}}\underset{}{{\displaystyle \frac{2^{2d/2}}{\mathrm{\Gamma }(\frac{d2}{2})}}\left({\displaystyle \frac{D}{\xi _+}}\right)^{\frac{d2}{2}}\text{K}_{\frac{d2}{2}}(D/\xi _+)},`$ (53)
$`1,TT_c`$
where $`\text{K}_{(d2)/2}`$ is a modified Bessel function and where $`\xi _+=l_H`$ for the case $`\tau =0`$ and $`H0`$. Thus as expected two point-like perturbations of the bulk do not exhibit a rich interaction structure but a monotonous decay. If a small sphere is immersed at a distance $`D`$ from another small sphere with a given adsorption profile, the resulting perturbation of this profile leads to an effective interaction. At those distances, at which the small radius expansion can be applied the profile of the distant sphere is already close to the bulk value. This differs from the system studied in ref.  where a small sphere in front of a planar wall is considered. In that case the sphere is exposed to the half-space adsorption profile which is a much stronger perturbation of the bulk order parameter than that of a sphere (compare refs.  and ).

### F Comparison with dispersion forces

It has been stressed above that the Casimir forces depend sensitively on both temperature and composition of the solvent. The Casimir forces add to the omnipresent dispersion forces, which also depend on the thermodynamic state of the solvent. However, this latter dependence is smooth and - within the critical region of interest here - it is weak. Therefore the dependence of the Casimir force on temperature and concentration allows one to extract it from the actual total force which is the only one which is experimentally accessible. As an example for a binary liquid mixture we consider a water-lutidine mixture for which the critical temperature is $`T_c=307.1\text{K}`$ . The strength of the dispersion forces is characterized by the Hamaker constant, which for a system of silica spheres immersed in such a mixture is approximately $`A=10^{20}\text{J}`$. Figure 9 shows the comparison of the dispersion forces and of the Casimir forces at the critical point for spheres with a typical radius $`R=100\text{nm}`$ which add up to the total force

$$F_{\text{tot}}(\mathrm{\Delta },0,0)=\frac{k_BT_c}{R}K_+(\mathrm{\Delta },0,0)+F_{\text{disp}}(\mathrm{\Delta }).$$
(54)
The dispersion forces are given by

$$F_{\text{disp}}(\mathrm{\Delta })=\frac{A}{R}\frac{32}{3}\frac{1}{(\mathrm{\Delta }+2)(\mathrm{\Delta }^2+4\mathrm{\Delta }+4)(\mathrm{\Delta }+4)^2\mathrm{\Delta }^2},$$
(55)
where $`A`$ is the aforementioned Hamaker constant. In the limit of small distances the Casimir force is known from the Derjaguin approximation

$$K_+(\mathrm{\Delta }0,0,0)=\pi \mathrm{\Delta }_{}\mathrm{\Delta }^2,$$
(56)
where $`\mathrm{\Delta }_{}=0.326`$ is the estimate of the Casimir amplitude in $`d=3`$ , and in the limit of large distances the Casimir force can be obtained from the small radius expansion

$$K_+(\mathrm{\Delta }\mathrm{},0,0)=2\frac{\beta }{\nu }𝒜_{}(\mathrm{\Delta }+2)^{2\beta /\nu 1},$$
(57)
where $`\beta /\nu =0.52`$ and $`𝒜_{}=\frac{(A_{}^\mathrm{\Phi })^2}{B_\mathrm{\Phi }}=7.73`$ .

In Fig. 9 also the forces between a single sphere and a planar wall are shown, as considered in ref. . The dispersion forces then read $`F_{\text{disp}}(\mathrm{\Delta })=\frac{A}{R}\frac{2}{3}\frac{1}{(\mathrm{\Delta }+2)^2\mathrm{\Delta }^2}`$, where $`\mathrm{\Delta }`$ is the surface-to-surface distance in units of the radius $`R`$. For small distances the amplitude of the Casimir force is two times bigger than in the case of two spheres while for large distances the small radius expansion leads to

$$K_+(\mathrm{\Delta }\mathrm{},0,0)=2\frac{\beta }{\nu }𝒜_{}(2\mathrm{\Delta })^{\beta /\nu 1}.$$
(58)
For small distances all forces increase as $`\mathrm{\Delta }^2`$, only the amplitudes differ from each other. The values considered for $`A`$ and $`T_c`$ lead to Casimir forces that are about ten times larger than the dispersion forces. For large distances, however, the (unretarded) dispersion forces $`(\mathrm{\Delta }^7)`$ decay faster than the Casimir forces $`(\mathrm{\Delta }^{2.04})`$. (Ultimately retardation of the dispersion forces leads to a decay $`\mathrm{\Delta }^8`$ and $`\mathrm{\Delta }^5`$ for SS and SW, respectively. For $`\mathrm{\Delta }=1`$ the Casimir forces between the two spheres are about fifty times larger than the dispersion forces such that the latter contribution to the total force \[Eq. (54)\] can be neglected. For a sphere in front of a planar wall and for large distances the Casimir forces $`(\mathrm{\Delta }^{1.52})`$ decay even slower than they increase for small distances $`(\mathrm{\Delta }^2)`$.

### G Equation of state

For the demixing transition of a binary liquid mixture the difference $`(\mu _A\mu _B)(\mu _A\mu _B)_c`$ of the chemical potentials of the two components is proportional to the bulk field $`h`$. The potential difference leads to a concentration $`c_Ac_{A,c}`$ of species $`A`$ in the bulk. The bulk order parameter $`M_b`$ is proportional to the concentration difference $`M_b=𝒜(c_Ac_{A,c})`$, where the amplitude $`𝒜`$ depends on the actual definition of the order parameter. For such a system in the vicinity of its critical point the equation of state is given by

$$h=𝒟M_b^\delta f\left(t\left|\frac{B}{M_b}\right|^{1/\beta }\right),$$
(59)
where $`\beta `$ and $`\delta `$ are the standard critical exponents, $`𝒟`$ and $`B`$ are nonuniversal amplitudes, and the scaling function $`f(x)=1+x`$ describes the crossover between the critical behavior at $`t=0`$ and $`h=0`$, respectively. Eq. (59) yields $`M_b=\pm B|t|^\beta `$ for $`h0`$ and $`M_b=(\frac{h}{𝒟})^{1/\delta }`$ for $`t0`$ as limiting cases. The choice of the amplitude $`𝒟`$ defines a field $`h`$. Since, however, the combination $`e^{{\scriptscriptstyle 𝑑VhM_b}}`$ is the field contribution of the statistical weight of the system, $`𝒟`$ is fixed by the amplitude of the chosen order parameter.

Away from the critical point the two point correlation function decays for large distances $`r`$ exponentially as $`G_b(r)\mathrm{exp}(r/\xi )/r^{(d2)/2}`$, with the true correlation length $`\xi `$. Within mean-field theory, where one has $`\beta =\frac{1}{2}`$ and $`\delta =3`$, this correlation length is equivalent to that one derived from the second moment of the correlation function, which in momentum space is given by $`\stackrel{~}{G}_b(q)\frac{1}{q^2+\xi ^2}`$. From this relation one infers

$$\xi =\left(\tau +3m_b^2\right)^{1/2},$$
(60)
with $`m_b`$ obtained from the mean-field equation of state:

$$H=\tau m_b+m_b^3.$$
(61)
Equations (60) and (61) lead to the mean-field expressions of the correlation lengths $`\xi _\pm `$ and $`l_H`$ in Eqs. (43) and (44), respectively. The correlation lengths are experimentally accessible by measuring the scattering structure factor. This allows one to determine the nonuniversal amplitudes $`\xi _0^\pm `$ and $`l_0`$ in Eqs. (43) and (44). By measuring the coexistence curve $`M_b=\pm B|t|^\beta `$ the nonuniversal amplitude $`B`$ is obtained.

However, the amplitude $`l_0`$ in Eq. (44) is fixed once $`\xi _0^+`$ and the nonuniversal amplitude $`B`$ are known:

$$l_0=\xi _0^+\left(\frac{B}{𝒜c_{A,c}}\right)^{\nu /\beta }\left(\frac{Q_2}{\delta R_\chi }\right)^{\nu /\gamma },$$
(62)
where the first term in parantheses is nonuniversal \- but it does not depend on the definition of the order parameter - and $`Q_2`$ and $`R_\chi `$ are universal amplitude ratios  leading to $`\left(Q_2/\delta R_\chi \right)^{\nu /\gamma }0.38`$ in $`d=3`$.

Using Eq. (43) and $`m_b=\sqrt{|\tau |}`$ one finds a relation between the order parameter $`m_b`$ in the present model and the order parameter $`M_b`$ defined in an actual experiment:

$$m_b=\frac{1}{(\sqrt{2}\xi _0^{})^{\beta /\nu }}\frac{M_b}{B}.$$
(63)
We note that in a magnetic system the bulk field $`H`$ is directly accessible, while for liquids we use the correlation length $`l_H`$ as a measure of the bulk field $`H`$ according to Eq. (44). So far we have assumed that the coexistence curve close to the critical demixing point is symmetric with respect to the temperature axis. However, this might not be the case in a real liquid mixture, so that one would then have to use appropriate linear combinations of the thermodynamic variables.

## V Results for the film geometry

In the previous section \[see Eq. (47) in Subsec. IV D\], within the Derjaguin approximation, we have used results obtained for the film geometry in order to analyze the limiting case of the colloidal particles being very close to each other. Besides that the film geometry is of considerable interest in its own right as the generic case for simulations  and for experimental studies involving wetting films  or force microscopes with crossed cylinders of large radii of curvature. The Casimir forces in the film geometry have been considered for walls that do not break the symmetry of the order parameter  and also for symmetry breaking walls . However, most studies focus on the temperature dependence only, while in the following we consider also the case $`H0`$.

For the film geometry the same numerical methods can be used as described in Sec. III. However, since the system is effectively one-dimensional there are additional methods available. For instance, the first integral of the second-order differential equation (9) can be found analytically . The integration constant, which follows from fulfilling the boundary conditions is proportional to the force (see App. A in ref. ). For $`H=0`$ the second integration leads to elliptic integrals, while for $`H0`$ it must be carried out numerically.

First we obtain the scaling function $`K_\pm ^{||}(\frac{L}{\xi },\frac{L}{l_H})`$ of the singular contribution to the force between two parallel plates at distance $`L`$ for the various thermodynamic paths shown in Fig. 2. The temperature dependence (solid path) of the force is discussed in refs.  and  (see in particular Fig. 3 in ref. ) and reveals a maximum of the force at $`\mathrm{\Theta }_{\text{max}}=L/\xi _+=1.94`$. As a new result the dependence of the force on the bulk field $`H`$ (dashed path) is shown in Fig. 10 in terms of the scaling variable $`\mathrm{\Xi }=\text{sgn}(H)L/l_H`$.

It is remarkable that the force maximum is about ten times larger than the force at the critical point. Since this scaling function enters into the expression for the force between two spheres this peak structure carries over to that case, too. Indeed, Fig. 7 shows that the force maximum in the limit of two very close spheres is about five times higher than the force at the critical point.

Figure 11 shows the dependence of the force on the bulk field for $`TT_c`$. Above $`T_c`$ one finds a similar structure like at $`T_c`$ except that the maximum is less pronounced. (Here we refer to the maximum value of the modulus of the force. However, in order to facilitate an easy comparison with the results of ref.  in Fig. 11 the sign of the force is chosen such that it reveals its attractive nature.) Slightly below $`T_c`$ the maximum is more pronounced although the form of the curve is still the same. But at still lower temperatures the character of the curve changes. Instead of a smooth increase upon increasing the bulk field towards $`H=0`$ there is a jump from a small value of the force to a high one. This is accompanied by a first-order transition of the corresponding order parameter profile.

Figure 12 shows the two coexisting profiles at the critical field of this first-order phase transition for various temperatures $`T<T_c`$. The profiles, which are negative in the middle of the film, correspond to the weaker of the two respective forces. This means that the forces are weaker if the solvent component preferred by the surfaces is the one in which the bulk solution is poor. These first-order phase transitions correspond to capillary condensation in the slit which can occur for fixed temperature as function of $`H`$ as shown in Figs. 11 and 12 or as function of temperature for fixed values of $`H`$. The loci of these capillary condensation phase transitions are shown in Fig. 13. The transition line ends in the capillary critical point ($`\mathrm{\Theta }_{}=7.17,\mathrm{\Xi }=9.54`$) which is in quantitative accordance with the findings for the critical point shift .

In ref.  the solvation forces between two parallel plates confining a two-dimensional Ising spin system have been calculated. They exhibit the same features as those found in the present mean-field calculations valid for spatial dimensions $`d4`$. For temperatures far below $`T_c`$ upon increasing the field the force jumps from a small to a large absolute value and the varies linearly as function of the bulk field (compare Fig. 11). In the linear regime the solvation force is given approximately by $`K_{}^{||}/L^4=2\overline{m}(T)H`$ , where the bulk spontaneous magnetization $`\overline{m}(T)`$, i.e., the magnetization at $`H=0`$ is $`\overline{m}=|\tau |^{1/2}`$ within mean-field approximation which indeed gives the slope of the linear parts in Fig. 11. We note that in $`d=2`$ capillary condensation can be identified only as a quasi-phase-transition whereas in $`d=4`$ there is indeed a bona fide first-order phase transition. Figure 14 shows the difference $`\delta K_{}^{||}`$ between the two corresponding force values at the transition. Upon approaching the capillary critical point this difference vanishes. Along the coexistence line $`H(T)`$ of capillary condensation the quantity $`2\overline{m}(T)H(T)`$ provides a good account of the force difference (see Fig. 14).

The findings in $`d=2`$ and $`d=4`$ agree to a large extent qualitatively. Quantitatively, however, the results in $`d=2`$ and $`d=4`$ differ significantly. For example the ratio of the maximum force and of the force at the critical point as function of the bulk field and for $`T=T_c`$ is $`K_+^{||}(0,\mathrm{\Xi }_{\text{max}})/K_+^{||}(0,0)100`$ in $`d=2`$ and $`K_+^{||}(0,\mathrm{\Xi }_{\text{max}})/K_+^{||}(0,0)10.1`$ in $`d=4`$. A simple linear interpolation yields an estimated value of $`K_+^{||}(0,\mathrm{\Xi }_{\text{max}})/K_+^{||}(0,0)55`$ for this ratio in $`d=3`$. Thus the studies in ref.  and the present ones can be used to obtain quantitative estimates for the actual behavior in $`d=3`$.

## VI Concluding remarks and summary

We have studied the critical adsorption profiles and the effective free energy of interaction for a pair of colloidal particles which are immersed in a binary liquid mixture near its critical demixing point. Our results pertain to the whole vicinity of the critical point of the binary liquid mixture, i.e., they are given as function of both the reduced temperature $`t=(TT_c)/T_c`$ and the field $`h`$ conjugate to the order parameter (see Fig. 2 and Sec. IV).

In an ensemble of colloidal particles dissolved in a near-critical solvent, the critical Casimir forces between them are expected to lead to flocculation of the particles. This holds even in cases where dispersion forces alone are not strong enough to produce such a phase transition (see Subsec. IV F). Indeed, such a flocculation of colloidal particles has been observed experimentally for various binary liquid mixtures . In particular, these experiments exhibit an asymmetry in the shape of the observed flocculation diagrams in that flocculation occurs on that side of the phase diagram where the binary liquid mixture is poor in the component preferred by the colloids. Our results are consistent with this observation, since we indeed find that the critical Casimir forces are strongest in this region of the phase diagram (see Subsecs. IV B, IV C, and IV E).

In the following we summarize the main results of this work:

1. For two spherical particles of radius $`R`$ at a distance $`D`$ (see Fig. 1) we have investigated the effective interaction mediated by the solvent, such as a binary liquid mixture. This effective interaction depends on the thermodynamic state of the solvent, for which the phase diagram is sketched in Fig. 2.
2. The solvent is described by a field-theoretic model, provided by the Ginzburg-Landau Hamiltonian \[Eq. (2)\], with additional terms representing the solute spheres. We solve this model within mean-field theory as described in Sec. III.
3. (III A) The profile at the critical point is obtained by a conformal mapping of the corresponding profiles between concentric spheres. (III B) Close to the surface of either sphere the short distance expansion of the profile is valid, which is carried out off the critical point. (III C) For the determination of the complete profiles we first discretize the space using the conformal lattice shown in Fig. 3 and then apply the method of steepest descent for the order parameter at the lattice points. (III D) From the order parameter profiles the force between the spheres is calculated by integrating the stress tensor. The choice of this integration path is indicated in Fig. 4.
4. The results for the two-particle interaction are presented in Sec. IV. (IV A) The dependence of the force on the temperature at the critical composition is shown in Fig. 5. It exhibits a maximum above $`T_c`$ which approaches that of the critical point for large distances. (IV B) The dependence of the force on the composition at the critical temperature is shown in Fig. 6. It exhibits a pronounced maximum whose position also approaches that of the critical point for large distances. The values of the maximum of the force determined for various distances are shown in Fig. 7. For small distances these values differ significantly from the value at the critical point. (IV C) In Fig. 8 the dependence of the force on the temperature for fixed compositions off the critical one is shown. (IV D) Using the knowledge of the force between parallel plates the Derjaguin approximation can be applied for small particle distances. (IV E) The small radius expansion can be applied for large particle distances regarding the spheres as perturbation of the bulk phase. This expansion yields nonperturbative results for $`d=3`$. In leading order the small radius expansion leads to forces which exhibit a maximum only at the critical point. (IV F) In Fig. 9 a comparison of the Casimir forces and the dispersion forces at $`T_c`$ and in $`d=3`$ is shown. For large distances the dispersion forces decay faster while for small distances the amplitudes are smaller than those for the Casimir forces. Moreover, the temperature and composition dependences of the Casimir forces allow one to distinguish them from the dispersion forces. (IV G) Relations between thermodynamic quantities as they are used in the present model and their experimental counterparts are discussed.
5. In Sec. V results for the case of parallel plates are presented, which are used via the Derjaguin approximation for the limiting case of small particle distances. Moreover, the present mean-field results are in qualitative agreement with corresponding results for a two-dimensional Ising spin system obtained in ref. . Figure 10 describes the dependence of the force on the bulk field $`H`$ at the critical temperature. In Fig. 11 the dependence of the force on the bulk field is shown also for temperatures considerably below $`T_c`$. As $`H`$ is increased the force jumps from a small absolute value to a large one. This is accompanied by a first-order phase transition of the corresponding order parameter profiles as shown in Fig. 12. These first-order phase transitions correspond to capillary condensation in the slit and their loci are displayed in Fig. 13. The transition line ends in the capillary critical point. In Fig. 14 the size of the jump of the force at capillary condensation (see Fig. 11) is plotted as function of the bulk field. Away from the capillary critical point this difference is well approximated by an estimate in which only the order parameter value for $`H=0`$ at the transition temperature and the corresponding bulk field enter.

## Acknowledgments

We thank E. Eisenriegler and A. Maciołek for helpful and stimulating discussions. We acknowledge partial financial support by the German Science Foundation through Sonderforschungsbereich 237 (FS) and grant No. HA3030/1-2 (AH). AH also acknowledges financial support by the National Science Foundation through grant 6892372 and by the Engineering and Physical Sciences Research Council through grant GR/J78327.