Effective Pair Interaction of Patchy Particles in Critical Fluids

We study the critical Casimir interaction between two spherical colloids immersed in a binary liquid mixture close to its critical demixing point. The surface of each colloid prefers one species of the mixture with the exception of a circular patch of arbitrary size, where the other species is preferred. For such objects we calculate, within the Derjaguin approximation, the scaling function describing the critical Casimir potential, and we use it to derive the scaling functions for all components of the forces and torques acting on both colloids. The results are compared with available experimental data. Moreover, the general relation between the scaling function for the potential and the scaling functions for the force and the torque is derived.


I. INTRODUCTION
Colloids have been the subject of research for centuries [1][2][3][4]. The initial studies were mostly concerned with the observation and explanation of the behavior of naturally occurring colloids in suspensions, as they are small enough to exhibit certain properties typical for molecular systems and, at the same time, they are big enough to be directly observable by using a microscope. With the development of the corresponding theoretical description [4][5][6] and methods of synthesis [5,7] it has become possible to design colloidal particles exhibiting desired properties; this has found applications in various areas like pattern formations [5,7], drug delivery [8,9], phoretic motors [10], in the oil industry [11], and numerous others [12]. In many cases, the description of the system can be simplified by introducing an effective interaction between the colloids, which is mediated by the solvent and is present in addition to their direct interaction. Different types of such effective interactions have been proposed [4,6] like London-van der Waals forces, screened electrostatic repulsion, steric, and depletion forces. One of the ways of inducing and controlling the interaction between colloids is the use of the critical Casimir effect [13,14].
The critical Casimir force is one of the manifestations of effective interactions induced by fluctuations. Historically, the first example of recognizing such a force was the electromagnetic Casimir effect [15], according to which the force acting between two conductors originates from the confinement induced restrictions of the quantum fluctuations of the electromagnetic field. In the case of the critical Casimir effect, two or more objects are immersed in a medium, in which the thermodynamic state is tuned to be close to its critical point. For these latter systems, the order parameter fluctuations are large, and the * nimabafi@is.mpg.de † pionow@is.mpg.de resulting effective interaction between colloids is longranged [13]. One of the most interesting features of critical phenomena is the concept of universality [16], which stipulates that in the vicinity of the critical point certain properties of the system (such as critical indices, ratios of amplitudes, and scaling functions) depend only on certain general features, like the spatial dimension and the number of components of the order parameter, but not on the microscopic details of the system. This allows one to group specific critical systems into various so-called bulk universality classes, which provides a convenient means of theoretical analysis: instead of modeling a complicated system one can study a much simpler model which can act as a representative of the corresponding universality class. One of the prominent examples is the 3D Ising universality class which contains a simple fluid close to its critical point, a uniaxial ferromagnet in the vicinity of the Curie point, and a binary liquid mixture close to its critical demixing point [17,18].
In the case of semi-infinite systems, each bulk universality class splits into several surface universality classes, depending on some general properties of the interaction between the critical system and the confining wall. Similarly, for systems with a slab geometry there emerge various film universality classes which can typically be characterized by pairs of surface universality classes of the two confining walls. The universality of the critical Casimir force manifests itself via the occurrence of scaling laws: close to the critical point the force can be expressed in terms of a power law times a universal scaling function. The latter depends only on dimensionless ratios of geometrical parameters describing the system, the bulk correlation length, and suitable scaled bulk fields. Universality of the scaling functions implies that they are identical for systems from the same bulk and film universality classes [17,18].
Recent advances in synthesis have facilitated the fabrication of colloidal particles in a controlled way with spatially varying surface properties. Since the effective interactions between such particles are anisotropic, the resulting self-assembly patterns can be much more complex than in the case of chemically homogeneous particles [5,7,14,61,62]. It was observed that the properties of critical Casimir interactions allow one to change the self-assembly structure of certain colloids by varying the thermodynamic parameters of the solvent such as temperature and concentration [59,60]. These observations have also been confirmed by various numerical simulations of chemically inhomogeneous particles [14,35].
A full understanding of the relation between the properties of a single colloidal particle and the pattern formed by a large number of such particles can potentially provide a useful tool to create any kind of three-dimensional microstructures. In this respect, one of the necessary steps is to investigate the critical Casimir pair interaction for inhomogeneous particles. Some theoretical studies [34,40,[63][64][65][66] have already addressed this issue by using mean field theory [67], the Derjaguin approximation [68], and the exact two-dimensional solution [69]. So far, these studies were devoted to either patterned surfaces or particles with distinct chemical properties on half of their surface (so-called Janus particles [61]).
In the present contribution, we extend the studies of equilibrium critical Casimir interactions to the case of two identical spherical colloids with a circular cap forming a chemical patch of arbitrary size on their surfaces. We use the Derjaguin approximation in order to determine all components of the force and the torque acting on them. These calculations can straightforwardly be generalized to more complicated chemical patterns. We compare our results with available experimental data.
The paper is organized as follows: In Sec. II we introduce spherical colloids with chemically inhomogeneous surfaces. Section III is devoted to critical Casimir effect; we recall the pertinent results for the slab geometry and for two spheres. Moreover we discuss the Derjaguin approximation used in our calculations. In Sec. IV the procedure of calculating forces and torques acting on the colloidal particles is described. Afterwords, in Sec. V we comment on the validity and accuracy of the approximation used in our calculations. Next, in Sec. VI we present our numerical results and compare them with available experimental data. Finally, the summary of our research is presented in Sec. VII. Our presentation is supplemented by three appendices. In Appendix A we derive the formulae relating the scaling function for the critical Casimir potential with the scaling functions for the com-ponents of the critical Casimir force and torque; in Appendix B we recall the scaling functions for the critical Casimir force in the slab geometry; and in Appendix C we discuss nonanalyticities of the scaling functions.

II. PATCHY PARTICLES
We consider a system of two spherical colloidal particles immersed in a binary liquid mixture close to its critical demixing point. We assume that the mixture consists of two species A and B, the composition of which is equal to the one of the critical point, and that the temperature T is close to the critical one T c . The two colloidal particles are spheres of the same radius R and the surfaceto-surface distance between them is denoted as D. In this study we assume that the forces and torques acting between the colloids are balanced by external forces such that the particles are kept in fixed positions and the system is in thermodynamic equilibrium. We do not consider any dynamic effects.
The surface of the particles is inhomogeneous, i.e., the interaction with the two components A and B of the mixture depends on the position on the surface. We study the critical Casimir interaction in the scaling limit (cf. Sec. III), in which only general properties of the wall-fluid interaction are relevant. In order to describe the surface it is sufficient to specify at each point which component of the mixture is preferred. In all figures, the regions where the component A of the mixture is preferred are denoted by '+' and plotted in red while the preference for component B is denoted by '−' and plotted in blue. The detailed interaction between the mixture and the surface of the spheres gives rise to subdominant terms in the scaling limit, i.e., corrections to scaling; studying them is beyond the scope of the present analysis.
We note that the approach used here renders all crossovers between regions of different affinity of the surface to be sharp. A more realistic, gradual description of such interfaces calls for a separate study.
In order to fully describe the configuration of the system, we assume that the center of the first colloid is located at the origin of the laboratory reference frame O, and the center of the second one is at the point defined by a vector r of length D + 2R. The rotational configuration of each colloid is determined by three angles α, β, and γ in accordance with the following procedure: The colloid is initially put with its center at the origin of O in a predefined initial configuration. First, it is rotated around the z axis by the angle α. The second rotation is by the angle β around the x axis, and the third one is by the angle γ around the y axis. Finally, in order to obtain the desired configuration, the second colloid is shifted by the vector r. All three rotations are active and in the direction determined by the right-hand rule. The procedure is shown schematically in Fig. 1. All applied rotations are represented by the matrix R (α, β, γ) =    cos α cos γ + sin α sin β sin γ . . . − sin α cos γ + cos α sin β sin γ .  Schematic plot of the three rotations which define the rotational configuration of the colloid. The particle with an arbitrary pattern (a) is, firstly, rotated by the angle α around the z axis (b); secondly, it is rotated by the angle β around the x axis (c); and, thirdly, it is rotated by the angle γ around the y axis (d). The pattern shown in this picture is chosen in order to illustrate all the transformations, but it is not studied in the present analysis.
In order to obtain any possible rotational configuration, it is sufficient to consider α ∈ [0 • , 360 • ), β ∈ (−90 • , 90 • ], and γ ∈ [0 • , 360 • ). The case of β = 90 • requires special care, because in this case the rotations around the z axis and the y axis are, in fact, the same rotation and thus one can assume that γ = 0. We note that it is often helpful to consider α, β, or γ beyond the domains given above; in such cases the resulting rotations can always be replaced by the ones which fulfill the constraints.
It is convenient to introduce in order to denote the rotational configuration of the colloids. Here, α i , β i , and γ i describe the configuration of the first (i = 1) and of the second (i = 2) particle. Accordingly, the system is described by four quantities: temperature T , radius R of the particles, relative position r of the colloids, and the rotational configuration Ω.
The surface-to-surface distance follows from the relation (a) Schematic plot of two spherical colloids in a special configuration (see the main text). In this configuration the center of the first particle is located at the origin and the center of the second particle is on the positive part of the y axis. The surface-to-surface distance between the colloids is denoted as D. (b) Schematic plot of the spherical particle with a single circular patch. This type of particles is the one considered in all subsequent calculations. The angle θp defines the size of the patch.
In order to study the system of two colloids it is not necessary to consider all possible configurations, because the critical Casimir interaction is invariant under rotations. (The translational symmetry has already been utilized by keeping the first particle at the origin.) We introduce the rotation T which moves the center of the second particle to the positive y semi-axis, i.e., where e y is the unit vector in y direction. With such a rotation, not only the vector r is transformed, but also the rotational configuration Ω; we denote the new configuration by We note that T is not defined uniquely by the relation in Eq. (3); composing T with any rotation around the y axis preserves Eq. (3). Since the additional rotation of the colloids around the y axis is only increasing γ T 1 and γ T 2 by the angle of the rotation, it is possible to choose T in such a way, that γ T 2 = 0. This additional condition renders T unique. Concerning an example of constructing the matrix T see the calculation presented in Appendix A.
We call the configurations of colloids, for which r = r e y and γ 2 = 0, special configurations. As we have shown above, any general configuration of two patchy particles can be transformed to the special one by the rotation T constructed above. Therefore, it is sufficient to calculate the critical Casimir forces and torques only for the special configurations. An example of such a configuration is presented in Fig. 2(a). We note that, because forces and torques are vectors, in order to calculate them for an arbitrary configuration one must transform the results obtained for the special configuration by the rotation T −1 .
For any configuration of the colloids we additionally introduce the relative rotational configuration as the configuration the particles would have if they were transformed to the special configuration. All angles in Ω * are (rather complicated) functions of r and Ω. Wherever possible, we introduce the laboratory reference frame O such that the particles already assume the special configuration. This way, we do not have to determine the rotation matrix T which simplifies the calculation.
The above discussion is valid for an arbitrary pattern on the surfaces of the colloids (and can easily be generalized to nonspherical particles). Here, we consider only a simple pattern with one circular patch on the surface of each sphere: in the initial position of the particle (before applying any of the rotations), the component A (B) of the mixture is preferred for θ θ p (θ > θ p ) where θ is the polar angle of the spherical coordinates in O, and the opening angle θ p is the parameter describing the angular size of the circular patch. A schematic plot of the colloid is presented in Fig. 2(b). For such a pattern the first rotation (by an angle α around the z axis) does not change the configuration and thus Ω is defined completely by β 1 , γ 1 , β 2 , and γ 2 .
The critical Casimir interaction in the special case θ p = 90 • has already been addressed in Ref. [34]. This corresponds to a so-called Janus particle [70,71], in which the particle consists of two hemispheres preferring opposite components of the binary liquid mixture.

III. CRITICAL CASIMIR FORCE
In order to calculate the critical Casimir interaction between colloids we use the Derjaguin approximation which is based on the slab geometry. In this section, we recall all pertinent results for the slab and the spherical geometries and adapt them to the present case of chemically inhomogeneous surfaces.

A. Slab geometry
We start from the description of the thermodynamic state of the binary liquid mixture. In general, such a liquid is fully described by three intensive parameters. In the current study we assume that the pressure p is fixed and the concentrations (molar fractions x A and x B = 1 − x A ) of the components of the mixture are tuned to be equal to the critical ones. The temperature T , as the third parameter, is free to change. We assume that it is close to the critical temperature T c . (The values of both the critical concentration and the critical temperature depend on the pressure p.) If such a liquid is confined by two macroscopically large, parallel walls, the resulting slab system is described by only two macroscopic control parameters: the temperature T and the distance L between the walls.
Close to the critical point, the critical fluctuations of the concentration of the fluid lead to the effective critical Casimir force [13] acting between the walls. In spatial dimension d = 3 this force can be described by the universal scaling formula where k B is the Boltzmann constant, F slab c /A is the critical Casimir force per area (i.e., excess pressure); ϑ s (ω) is a scaling function which is universal, i.e., it is independent of the microscopic details of the system. The scaling function depends only on the bulk universality class of the critical point of the fluid, and on the interaction between the fluid and the two walls, encoded into the corresponding film universality class [18], denoted by the index 's' (see below). The scaling variable is where in the case of an upper critical point of a binary liquid mixture, t = (T − T c ) /T c is the reduced temperature, and ν is the critical exponent of the correlation length.
The expression in Eq. (6) is exact in the scaling limit, i.e., for T → T c and L → ∞ with ω fixed. If T is fixed and close to T c , and L is large but finite, Eq. (6) provides an approximation of the actual critical Casimir pressure; in order to improve the result, one has to include corrections to scaling (such as higher order terms in 1/L).
In the present context we are interested only in the surface universality class of a symmetry breaking surface field in which the surface prefers either A ('+') or B ('−') chemical species of the binary liquid mixture. For our study only two film universality classes are relevant: If both walls prefer the same component of the mixture the scaling function in Eq. (6) is ϑ sm (ω) (s = sm, same boundary conditions: '++' or '−−') and if the walls prefer different components of the binary mixture it is ϑ op (ω) (s = op, opposite boundary conditions: '+−' or '−+'). Since the analytical forms of these scaling functions are not known, we use their numerical estimates in spatial dimension d = 3 (see, c.f., Sec. IV A and Fig. 4).
Finally, we recall the scaling formula for the potential of the critical Casimir force: where ϕ s (ω) is a universal scaling function with s = sm or s = op; the form of this function can be determined from ϑ s (ω) by using the relation which follows directly from the relation between the critical Casimir force and its potential.

B. Spherical objects
When two spherical colloids are immersed into a critical fluid, like in the slab geometry, fluctuations induce a critical Casimir interaction between them. For homogeneous spheres, this effect has been studied theoretically [21,23], numerically [73], and experimentally [53]. There are three macroscopic parameters describing the system: the temperature T , the radius R of the colloids, and the surface-to-surface distance D between them. Following the literature, these variables are combined into the following two scaling variables: If the spheres are chemically homogeneous, due to symmetry the critical Casimir force acts in radial direction only. Close to the critical point, it is given by the scaling law where the index 'H' indicates that the considered quantity is evaluated for homogeneous spheres, F H s is a universal scaling function where the index s = sm or s = op denotes the same or opposite affinities of the two spheres, respectively. The additional factor ∆ −2 has been introduced in order to render the scaling function finite in the limit ∆ → 0 [21]. The potential of the critical Casimir force in the homogeneous case is given by where U H s is another universal scaling function. Like in the slab geometry, the scaling functions F H s and U H s are related. Equations (12) and (13) are valid in the scaling limit T → T c , D → ∞, R → ∞ with ∆ and Θ fixed.
We now turn to the case of inhomogeneous colloids studied here. If the preferences for the two components of the binary mixture vary along the colloid surface, the critical Casimir interaction is modified relative to the homogeneous case; in general, the force becomes non-radial and a torque appears.
In the special configuration (see Sec. II), Eq. (12) can be generalized to the following scaling formulae: where F (i) and T (i) denote the vector scaling functions for the force F (i) c and the torque T (i) c acting on the inhomogeneous spheres; i = 1, 2 labels the particles; Ω * denotes the relative orientation of the colloids (see Eq. (5)); and the scaling variables ∆ and Θ are given by Eq. (11). Note that we consider here the torques to be acting on the center of each particle.
Since the system is in thermal equilibrium, the total force and torque must vanish. This allows one to relate the above scaling functions: The formulae show that it is sufficient to calculate the force and the torque acting on one particle; the other quantities follow from Eq. (15). The critical Casimir potential of interaction between the colloids can also be determined in terms of an appropriate scaling law: where U is the scaling function and the factor ∆ −1 has been split off in order to keep the scaling function finite in the limit ∆ → 0. Unlike force and torque, the potential is a scalar quantity and thus the scaling formula in Eq. (16) holds even if the second colloid is not located on the y axis (because Ω * is the relative configuration of the spheres, it is invariant under the rotations of the system). The relation between the scaling function for the potential and the scaling functions for the forces and torques is derived in Appendix A and is given by Eq. (A8). Finally, we note that the scaling laws in Eqs. (14) and (16), together with the relation between the scaling functions studied in Appendix A, are very general and can be used for particles of arbitrary shapes and surface patterns.

C. Derjaguin approximation
In order to obtain the scaling functions for the critical Casimir interaction one can use several distinct tech-(a) Schematic plot of the (grayish) projection plane used in the Derjaguin approximation for an exemplary configuration of particles. Within this approximation the interaction between spheres depends only on points from the right hemisphere of the first particle and the left hemisphere of the second particle. Panel (a) illustrates the process of the orthogonal projection of the interacting surfaces (facing each other) to a plane parallel to the x and z axes which is the projection plane. The horizontal lines illustrate the projection for three points on the projection plane serving as examples. Panel (b) presents the resulting pattern on the projection plane. The projection of both spheres gives the same circle Λ of radius R, which we divide into four disjoint (nonoverlapping) regions Λ++, Λ+−, Λ−+, and Λ−−. The first sign in the index of Λ denotes the affinity of the surface of the first particle and the second sign of the second particle. The symbol '+' means that there is a patch at an appropriate point of the surface of the particle (denoted with red color) and '−' denotes that there is no patch (blue color).
niques like mean field theory [67], Monte Carlo simulations [74], or the so-called Derjaguin approximation [68,75]. All of these techniques provide only an approximation to the actual scaling functions: mean field theory is exact only for spatial dimension d 4 (with logarithmic corrections in d = 4); within the Monte Carlo simulations the size of the lattice is limited and therefore it is challenging to extract from the numerical data reliable results for the scaling limit; and within the Derjaguin approximation geometrical aspects of the system are not captured accurately. Out of these available methods, the Derjaguin approximation is the most straightforward scheme and requires the least numerical effort. Therefore, it provides an appropriate starting point for investigating the critical Casimir interaction in our system.
Within the Derjaguin approximation any available result for the planar geometry (typically the best one) can be used in order to estimate the effective interaction between more complicated objects. This approximation has been applied for many problems [75][76][77][78][79] and, in the case of the critical Casimir force, the results typically agree qualitatively (and under favorable circumstances even quantitatively) with the proper ones as far as they are available [23,30,34]. Here, we describe briefly the concept of this approximation (mostly in order to introduce those objects and quantities which turn out to be useful for our analysis); concerning the discussion of the validity of this approximation in our present case see Sec. V.
We assume that the particles are in the special configuration (i.e., the center of the first sphere is at the origin and the center of the second one is on the positive y semi-axis at the point (r x = 0, r y = 2R + D, r z = 0)). We introduce the projection plane, parallel to the x and z axes, and project orthogonally onto it the patterns on the hemispheres of both colloids facing each other (i.e., the right hemisphere of the first particle and the left hemisphere of the second one). The resulting figure is a circle Λ of radius R which we separate into four disjoint regions (sets) Λ ++ , Λ +− , Λ −+ , and Λ −− . The first and second sign in the index of Λ denotes the preference of the surface of the first and second particle, respectively. For example, for every point P ∈ Λ +− on the projection plane there are two points, one on each sphere, which are projected onto P ; the point on the first particle is in the patch while the point on the second particle is not. A complete example of the construction scheme described above is presented in Fig. 3. It is convenient to additionally define two sets Λ sm and Λ op , where the properties of the two points on the surfaces of both particles are the same and opposite, respectively: In order to make the definitions mathematically complete, it is necessary to define the surface affinity also at the edge of the patches. We assume that in these points the surface exhibits the same preference as the inside of the patch, i.e., the patch on the sphere is a closed set. As expected, our results do not depend on this convention. Within the Derjaguin approximation, for each point P on the projection plane, the distance between those two points which are projected onto P is calculated, and the contribution to the force of the surface element dA around the point P is estimated via Eq. (6) to be where both the distance and the scaling variable ω depend on P , and 's' denotes the pair of boundary condi-tions for the points projected onto P . The total force is obtained as the integral of Eq. (18) over the whole circle Λ on the projection plane: It is convenient to parametrize the circle Λ on the projection plane by using the spherical coordinates of the first particle 0 θ π and 0 φ π, where the range of φ is restricted because only the facing hemisphere of the first particle is used in the parametrization (see Fig. 3). For this choice of integral variables we have determined where R is the radius of the spherical colloids, D is their surface-to-surface distance, and ∆ and Θ are the scaling variables, given in Eq. (11). In order to simplify the notation, we have denoted the expression in Eq. (20c) byω; it is an argument of the scaling functions for the slab geometry when they are used to calculate the scaling functions for two spheres within the Derjaguin approximation. Using Eqs. (19), (20), and (14a), we have calculated the formula for the scaling function for the critical Casimir force acting on the first particle: where the variables θ and φ are the spherical angular coordinates on the first particle, s (θ, φ) denotes the same or opposite boundary conditions in the point parametrized by θ and φ (s depends on the configuration Ω * of the particles), and the argumentω of the scaling function is given by Eq. (20c). The force as given by Eq. (21) cannot be considered as a reliable approximation of the actual critical Casimir force acting between colloids with inhomogeneous surfaces. First, it always acts in the direction of the line connecting the centers of particles (i.e., it is a radial force). Second, within the Derjaguin approximation there is no straightforward way to determine the torques present in the system (besides the slab geometry in which there are none). Third, the interaction described by Eq. (21) is not conservative (because a radial force is conservative if and only if it does not depend on the angles).
In order to overcome the above problems, we use the Derjaguin approximation for deriving the potential of interaction instead of deriving the force. Straightforward calculation leads to whereω is given by Eq. (20c). The critical Casimir forces and torques are calculated as derivatives of the potential. Accordingly, by construction the obtained interaction is conservative. The formulae for the scaling functions F (i) and T (i) are provided in Appendix A. It is reassuring that the Derjaguin approximation for the force (Eq. (21)) renders the same expression for the radial component of the force as the corresponding derivative of the potential given in Eq. (22). The difference between these two approaches is that in addition the potential gives the torques and the non-radial components of the force in such a way that the interaction is conservative.

IV. METHOD OF CALCULATION
A. Scaling functions for the slab geometry In order to be able to calculate the scaling functions within the Derjaguin approximation, it is necessary to know the scaling functions for the slab geometry (see Sec. III A). They can be estimated, e.g., via Monte Carlo simulations or, alternatively, by using the extended de Gennes-Fisher local-functional method [80,81].
Here we use the data from the corresponding Monte Carlo simulations [73]. This technique allows one to estimate the scaling functions only for a limited number of values of the scaling variable ω. In order to obtain the full scaling functions ϑ sm (ω) and ϑ op (ω) it is necessary to interpolate and extrapolate the available data; the details of this procedure are described in Appendix B. We plot the resulting scaling functions in Fig. 4.
We note that the function ϑ sm (ω) is negative throughout (i.e., if both walls prefer the same component of the mixture, there is a critical Casimir attraction) whereas ϑ op (ω) is positive throughout (i.e., there is critical Casimir repulsion of walls preferring different liquid components). Additionally, the absolute value of the scaling function is larger in the case of opposing surface affinities.
Because of the unknown form of the leading corrections to scaling (see Ref. [73]), the interpolation of the Monte Carlo data and thus, the construction of the above scaling functions suffer from numerical errors. The systematic error is estimated to be up to 20% [34], which translates directly to all of our numerical results. In order to increase the precision, one can normalize all results by dividing them by the critical Casimir amplitude ϑ sm (0); such normalized functions have a numerical error of up to 5%. We note that in our calculations this systematic error, together with the inaccuracies of the Derjaguin approximation, is the main source of error; all other numerical inaccuracies present in our calculations are much smaller and can be neglected. On the other hand there is good reason to be confident about the reliability of the scaling functions shown in Fig. 4, because they agree excellently with high resolution experimental  6)) for the 3D Ising universality class, for opposite (ϑop) and same (ϑsm) boundary conditions as functions of the scaling variable ω (Eq. (7)). For further details see Appendix B.

B. Calculation of the scaling function for the interaction potential
The scaling functions ϕ s (see Eq. (9)) are prerequisites for calculating the critical Casimir potential of interaction for two spherical colloids (see Eq. (22)).
The calculation of the corresponding integral is based on the adaptive quadrature algorithm implemented in the GNU Scientific Library (GSL) [82]. In the course of carrying out the integral in Eq. (22), we first fix the value of θ and calculate the integral over φ: where I 0 (θ, φ) denotes the integrand in Eq. (22). We note that the function I 0 is discontinuous at all points where the surface boundary conditions change. Moreover, for certain values of θ it can happen that two points of discontinuity are located very close to each other and the change of integrand is easy to miss in the numerical integration. In order to avoid this problem, for each θ we calculate analytically all those values of φ, for which I 0 is discontinues and subdivide the integral as follows: where 0 < φ 1 < φ 2 < . . . < φ k < π denotes all points where the integrand has a discontinuity. This way all discontinuities of I 0 are properly taken into account in the course of the integration. Additionally, as all integrands on the right-hand side of Eq. (24) are now continuous functions of φ, this subdivision is reducing the time required for the numerical calculation.
Finally, we calculate the integral over θ in order to obtain the scaling function Here, we locate all values of θ, at which the patches start or end, and split up the integral accordingly. The resulting scaling function U (∆, Θ, Ω * ) is evaluated with a relative or an absolute error of 10 −6 , whichever is attained first.
The program performing the above algorithm of evaluation of the scaling function was written in C++. Further processing of the data was carried out using Mathematica [83].

C. Calculation of the forces and torques
In order to obtain the forces and torques acting on the colloids we use Eq. (A8). This implies that we have to calculate numerically the derivatives of the scaling function U (∆, Θ, Ω * ) for the potential. In this section we describe the corresponding procedure.
First, we note that it is not necessary to calculate the derivatives ∂U/∂∆ and ∂U/∂Θ. In the derivation of forces and torques both the radius R of the particles and the temperature T are fixed, and only the surfaceto-surface distance D can change (see Eq. (11)). Second, a close inspection of Eq. (A8) shows that the derivative ∂U/∂D appears only in the formula for the radial component of the force (see Eq. (A8b)), which can be calculated directly from the Derjaguin approximation for the forces (see Eq. (21)).
It remains to determine the derivatives of U with respect to the angles appearing in the relative configuration Ω * . In order to simplify the notation, in the following we discuss the calculation of ∂U/∂δ, where δ is one of the angles β * 1 , γ * 1 , or β * 2 . The value of δ, at which the derivative is calculated, is denoted as δ 0 . We also do not consider the dependence of U on all other variables as they are fixed.
The general procedure of calculating ∂U/∂δ| δ=δ0 is as follows: First, we fix a small positive numberε and a positive integer n. Second, we calculate the values of the function U for n values of δ uniformly distributed in the interval I = [δ 0 −ε, δ 0 +ε]. We denote these points in I by δ i for i = 1, 2, 3, . . . , n. Third, we fit the quadratic function to the points P i = (δ i , U (δ i )) (obtained in the second step) by using the method of least squares. This way we obtain the values of the coefficients a, b, and c. With them the estimate of the derivative is The quadratic term in the fitting function in Eq. (26) has been included in order to somehow account for the fact that in general the function U is not linear within the interval I. The method of least squares was used in order to reduce the error stemming from the chosen numerical integration scheme for calculating U; this error is similar to random noise. We note that other sources of error, like the inaccuracies of the scaling functions ϕ sm and ϕ op (see Eq. (9) and Appendix B), or inaccuracies of the Derjaguin approximation, are of a more systematic nature, and the fitting procedure leaves these errors unaltered.
The above general procedure needs to be adjusted near certain special values of δ 0 . One of these problems arises, if δ 0 is located close to the boundary of the domain of the function U, and certain values U (δ i ) cannot be calculated. In this case, for the fitting we use only those points δ i which are inside the domain. This way the number of points is reduced but it is still not smaller than n/2.
If there is a point of nonanalyticity of the function U inside the interval I, the fitting function in Eq. (26) might be not a good approximation. Around such a point, the general procedure fails and needs to be corrected. This is the reason why, in order to calculate the derivatives, it is necessary to investigate nonanalyticities of U. Below, in Sec. IV D, we present the observed types of singularities of the interaction and discuss how to adjust the general procedure described above.
If the point δ 0 is far away from the points of nonanalyticity, the above procedure gives, as we have checked, reliable results forε = 10 −3 and n = 21. These values have been used in order to calculate all results presented in Sec. VI.

D. Nonanalyticities of scaling functions
Upon changing the relative orientation Ω * of the colloids, we have observed three main types of singularities in the interaction of the colloids. In this subsection, we briefly characterize them and describe how the derivatives of the scaling functions around them can be calculated numerically. A more detailed mathematical analysis is presented in Appendix C.
For reasons of simplicity, we consider the particles in a special configuration (see Sec. II). In Fig. 5(a) we present the typical behavior of the scaling function for the potential U when one of the particles is rotated around the x axis; it exemplifies the three main types of singularities. In Fig. 5(b) we plot the derivative of the scaling function for the interaction potential, which is actually proportional to the x component T (2) x of the vectorial scaling function for the torque (see Eq. (A8h)). We note that in this figure we allow for β 2 > 90 • -in this case the rotation is equivalent to the one withβ 2 = 180 • − β 2 andγ 2 = 180 • . (Here the superscript ' * ' can be omitted because the initial configuration is already a special one.) The singularity of type I (type one) appears if the patches on both colloids form a so-called mirrorsymmetric configuration, i.e., for β 2 = −β 1 > −θ p and γ 1 = 0, like configuration E in Fig. 5(c). In this case, Λ sm = Λ and (within the Derjaguin approximation) the function U has the same value as for homogeneous spheres. If in this case any of the colloids is rotated, Λ op becomes a nonempty set and U increases. This produces a characteristic 'V' shaped cusp of the scaling function for the potential and a discontinuity of its first derivative (see Fig. 5 for β 2 = 90 • ). In Appendix C we show that the left-and the right-side derivatives of the function U at a point of type I singularity have the same absolute value but opposite signs.
In order to calculate numerically the derivative of U close to the nonanalyticity of type I, the procedure described in Sec. IV C has to be modified. For the fitting, instead of Eq. (26), we use where δ I denotes the value of the angle, for which the singularity occurs; a I , b I , and c I are fitting parameters. Exactly at δ I the derivative does not exist. However, because at δ I the potential reaches its minimum value, one can assume that the derivative at δ I is equal to zero, i.e., in the configuration with patches in mirror-symmetric configuration, there are no non-radial forces and torques. The occurrence of singularities of type I in the scaling function of the critical Casimir interaction, which is calculated within the Derjaguin approximation for Janus particles, has been reported in Ref. [34].
The singularity of type II occurs if the projections of the two patches on the projection plane are tangent, i.e., if there is only a single point in the region Λ ++ ; see the configurations C and G in Fig. 5(c). In order to analyze this case it is useful to introduce the overlap angle where ζ pp is the angular distance between the midpoints of the patches, i.e., the angle between n 1 and n 2 , where n i , for i = 1, 2, is a unit vector starting from the center of the i-th particle and ending at the surface of the same particle, in the midpoint of the patch. If ζ II overlap > 0, the patches overlap and the region Λ ++ is nonempty. When ζ II overlap < 0 there is no overlap and, on the projection plane, the two regions Λ +− and Λ −+ are separated by Λ −− . The nonanalyticity of U occurs for ζ II overlap = 0. In Appendix C we show that the nonanalytic part of U is proportional to with the corrections of the order of ζ II . This behavior makes the nonanalycities of type II hard to notice in the plots of the scaling function for the potential Up to the factor ∆, this quantity equals T (2) x , which is the x component of the scaling function for the torque acting on the second colloid. (c) The configurations of particles at special points. The letters A-H, denoting the configurations, are also marked at the top of the plots (a) and (b). In the configurations B-H the scaling function U is not analytic. In Sec. IV D we classify these nonanalyticities into three types: type I singularity (configuration E in the plots) where there is a 'V' shaped cusp nonanalyticity of U and where there is a jump of the derivative of U; type II singularity (configurations C and G) for which U exhibits a bump and for which the derivative of U exists but has an infinite slope; and type III singularity (configurations B, D, F, and H) where the second derivative of U has an infinite slope. See the main text for further details.
(especially if the amplitude of the nonanalytic term is small). On the other hand, the derivative of U behaves like a square root (i.e., a cusp with infinite slope) and thus in the plots of the forces and of the torques these singularities are well visible; see In order to calculate the derivative of U with respect to a certain angle δ in the region close to a singularity of type II, the procedure described in Sec. IV C has to be modified by replacing the fitting function in Eq. (26) with where a II , b II , c II , and d II are fitting parameters; δ II is the value of δ for which there is a nonanalyticity; Θ H denotes the Heaviside step function; and κ = +1 if ζ II overlap increases upon increasing δ, and κ = −1 otherwise. The additional term accounts for the nonanalytic part of U. We note that, unlike for the singularity of type I, at δ 0 = δ II the first derivative of U exists but the second derivative diverges from one side. Up to our knowledge, this is the first report of this type of nonanalyticity for the critical Casimir interaction calculated within the Derjaguin approximation.
Within the Derjaguin approximation, only half of each sphere takes part in the interaction, i.e., points on the right hemisphere of the first colloid interact with points on the left hemisphere of the second colloid. If the patch is fully located within one of the two other hemispheres, the energy of interaction does not depend on its precise position and the derivatives with respect to some of the angles in Ω * are zero. Accordingly, a nonanalyti-city of the function U is expected to occur if the patch passes through the great circle separating the two hemispheres. We call this type of singularity of the interaction as type III. It is expected to occur if the edge of the patch is tangent to the great circle separating the interacting and noninteracting hemispheres, i.e., if on the projection plane one of the sets Λ ++ ∪ Λ +− or Λ ++ ∪ Λ −+ have exactly one point in common with the edge of the circle Λ; this situation takes place if |β 1 | = θ p or |β 2 | = θ p . This nonanalyticity occurs for the configurations B, D, F, and H in Fig. 5(c).
Also in this case it is useful to define the two overlap angles where the patch-hemisphere angle ζ ph,1 (ζ ph,2 ) denotes the angle between the unit vector n 1 (n 2 ) pointing to the midpoint of the patch on the first (second) colloid (see the discussion after Eq. (29)), and the vector n h,2 = e y (n h,1 = −e y ) pointing to the midpoint of the noninteracting hemisphere on the second (first) colloid. The singularity of type III occurs in configurations for which ζ III overlap,i = 0 or ζ III overlap,i = 2θ p for i = 1 or i = 2. In Appendix C we argue that singularities of type III manifest themselves via a nonanalytic term in U of the form where the first term is relevant for |ζ III overlap,i | 1, while the second term holds for |ζ III overlap,i − 2θ p | 1. This means that the derivatives of U with respect to the angles (i.e., forces and torques) exhibit a nonanalyticity of the order of δ 3/2 .
For |ζ III overlap,i | 1 the nonanalytic term in the formulae for the force and torque is the leading order term and therefore it is well visible. In contrast, for ζ III overlap,i − 2θ p 1 the nonanalytic term is only a small correction and other analytic terms dominate. In the latter case, within our accessible numerical precision, we have not been able to verify that the function U contains this term.
Since the singularity of type III for the potential manifests itself through a nonanalyticity of the order of 5/2 -which is higher than the order of the polynomial in Eq. (26) used for the fitting (see Sec. IV C) -there is no need to adjust the general procedure in this case.
The singularities discussed above are relevant for calculating the derivatives with respect to angles and, via Eq. (A8), for all components of the scaling functions for the force and torque, except for the radial components of the force F (i) y . The resulting nonanalyticities of F (i) and T (i) can take three forms: jumps of the value, ∝ δ 1/2 , and ∝ δ 3/2 for the singularities of type I, II, and III, respectively (see Fig. 5(b)).
In the case of the radial components of the scaling functions for the force, the situation is different. They are calculated from Eq. (A8b), where there are no derivatives with respect to the angles. Therefore we expect that for these components the singularities of interaction should manifest themselves in the same way as they do for the scaling function U for the potential -the nonanalyticities of the form of a 'V' shaped cusp, ∝ δ 3/2 , and ∝ δ 5/2 for the singularities of type I, II, and III, respectively. In practice, these components are calculated directly from Eq. (21). The properties of the radial component of the scaling function for the force are presented in Sec. VI A.
Finally we emphasize that all three types of singularities discussed above follow from using the Derjaguin approximation in its present form; the scaling functions for the force and the torque onto objects of finite size are expected to be analytic functions of their rotational configuration and scaled distance. Nonetheless, it is useful to have an overview of the properties of the present, widely used, approximation scheme.

V. VALIDITY OF THE DERJAGUIN APPROXIMATION
Before presenting our numerical results it is necessary to discuss the reliability of the Derjaguin approximation as formulated in Sec. III C. Unfortunately, there is no systematic way to estimate the inherent error of this approximation. Moreover, up to our knowledge, so far the patchy particles considered here have not been studied by different techniques. Therefore we cannot estimate the accuracy of our results by making a simple comparison with other independent results. Instead, we look at similar models in order to identify possible shortcomings of the Derjaguin approximation.
We start with the case of homogeneous spherical particles of radius R separated by a surface-to-surface distance D. At T c for such a system, the potential following from using the Derjaguin approximation diverges ∝ D −1 for D → 0 and vanishes exponentially for D → ∞. General arguments based on conformal invariance [21] predict that for homogeneous spheres the critical Casimir potential diverges ∝ D −1 for small D, and decays ∝ D −β/ν for large D. This means that in this case the Derjaguin approximation fails to predict the long-ranged, algebraic decay. The detailed analysis shows that the Derjaguin approximation provides reliable results if D R [23,24,29]. In the scaling limit, the approximation is exact for ∆ = D/R → 0 [23,29] but it fails to reproduce the correct behavior for large ∆.
If the surfaces of the immersed objects are not homogeneous, the situation is different in that the Derjaguin approximation can give wrong results even for ∆ → 0.
So far, results are available for inhomogeneous wall-wall [64,66], sphere-wall [65], and cylinder-wall [34,40] geometries. Within mean field theory they are based on the numerical minimization of the corresponding Hamiltonian, which gives correct results (up to logarithmic corrections) in d = 4 spatial dimensions. Beside that, the approximation can be tested against Monte Carlo simulations, an exact solution in d = 2 [66], and, to some extent, experimental results in d = 3 [65].
The analysis of available data allows us to identify two general situations in which the Derjaguin approximation for the critical Casimir force is expected to give wrong results: (i) The characteristic length-scale of the pattern on the surface is much smaller than either the correlation length or the distance between walls. (ii) There is a region between interacting objects, where the gradient of the order parameter Φ (r) has a large component in the direction perpendicular to the direction n D used for the projections as they are applied within the Derjaguin approximation (in our case n D = e y ), i.e., in regions with Since obtaining the order parameter profile Φ (r) is typically at least as difficult as calculating the critical Casimir force, the above criterion is not useful for an exact analysis. Nevertheless, it is usually possible to roughly estimate the order parameter profile for a given system and then use Eq. (34) in order to check the validity of the Derjaguin approximation.
In situation (i), typically, the pattern on the surface is averaged and yields a certain effective homogeneous surface [43]; this aspect is not captured correctly by the Derjaguin approximation. Since in our present study we consider only colloids with a single patch, here this case is not relevant unless the patch is very small (θ p 1). In situation (ii), the free energy associated with the rapid change of the order parameter in all but one direction is completely missed by the approximation. The most prominent example of such a case is the capillary bridge formation at T < T c between patches on the two spheres. In this case, the force in the normal direction is increased by the contribution stemming from the surface tension of the interface present in the system [84][85][86]. This effect is not captured within the Derjaguin approximation. For T > T c the situation (ii) can occur if D is small, i.e., in the region where the Derjaguin approximation is actually expected to work quite well.
We note that the above criteria are only conjectures based on a limited amount of data available in the literature, and they do not provide a strict answer to which extent the Derjaguin approximation is reliable or not; they are supposed to identify regions, where the approximation can fail. This issue definitely deserves more research. An example for the case in which, despite of the large perpendicular gradients of the order parameter, the approximation can provide quite accurate results, concerns the lateral component of the force in a system with a capillary bridge. It can be correctly described within the Derjaguin approximation, except in the vicinity of the breaking transition [66].
It is also important to mention, that the nonanalyticities of the critical Casimir force reported above in Sec. IV D are all related to the discontinuous variation of the chemical surface properties. Such changes are known to not propagate in this form into the fluctuating medium [63]. Moreover, such nonanalyticities, up to our knowledge, have not been reported in calculations which are not based on the Derjaguin approximation. Therefore, we strongly expect that these singularities are an artifact of the Derjaguin approximation.
Finally, we note that both Monte Carlo simulations and mean field calculations for systems with two inhomogeneous spheres are numerically challenging. Up to our knowledge, such studies have not yet been reported, and most probably will give estimates of the scaling functions only for a rather limited number of points. The present results can provide guidance for how to interpolate these data points correctly, even in regions where the Derjaguin approximation is not working well.

VI. RESULTS
In this section we present our results for the critical Casimir interaction between patchy particles. Using the method described in Sec. IV, we have been able to calculate, within the Derjaguin approximation, the potential and all components of the forces and torques arising from the critical Casimir interaction between two spherical colloids with chemically inhomogeneous surfaces. In the special case of Janus particles (θ p = 90 • , i.e., the patch is covering half of the particle surface) some results have already been reported [34]; our analysis is in full agreement with them. Moreover, we have been able to extend these results by providing non-radial components of the critical Casimir force as well as of the critical Casimir torque.
In this section, we first discuss our results for the radial critical Casimir force and the potential, and compare them with those for the special case of Janus particles. Then, we present our results for the non-radial components of the critical Casimir force and for the torque. Finally, we present a comparison with corresponding experimental data.

A. Radial component of critical Casimir force
We start the discussion of our results by analyzing F (i) r for i = 1, 2, i.e., the radial component of the vectorial scaling function of the critical Casimir force acting on the first and second particle, respectively. Since F   Fig. 6. For all plotted curves we have chosen ∆ = D/R = 1, β 1 = γ 1 = 0, and θ p = 30 • . If β 2 = 0 (i.e., if the patch on the second particle is in the topmost position), the radial component is negative (i.e., attractive) for all values of Θ. Upon increasing β 2 (i.e., rotating the second particle such that its patch is moved towards the first particle), for any fixed Θ the radial component changes sign. This change occurs first for large negative values of Θ. If β 2 = 90 • (i.e., the patch is facing the first particle), F (2) r is fully positive (i.e., repulsive). Upon increasing β 2 further, the strength of the repulsion decreases and, eventually, attraction is recovered. This change starts at large positive values of Θ and moves towards negative values of Θ.
The observed behavior of the radial component of the scaling function of the critical Casimir force can be easily understood. If β 2 = 0, within the Derjaguin approximation, only the '++' and the '−−' boundary conditions are active, so that due to ϑ sm (ω) < 0 the resulting net radial component is negative. If β 2 is increased, the region Λ −+ on the projection plane (see Sec. III C) emerges (i.e., the set Λ −+ becomes nonempty) and, because ϑ op (ω) > 0, the force becomes less negative; if the patch is sufficiently large, the sign of the force eventually changes. The area of the region Λ ++ decreases upon increasing β 2 , and Λ ++ becomes empty for β 2 > 2θ p . Starting from there, the dependence of the radial component of the scaling function on β 2 is solely determined by the location of the patch on the second particle; therefore we focus on the region Λ −+ . For β 2 = 90 • , the mean distance between the points on the two spheres projected onto the region Λ −+ (see Sec. III C) is smallest, and, moreover, the area of Λ −+ is largest; thus the repulsive contribution is strongest. Increasing β 2 further moves the patch to the bottom, where the mean distance is large and, concomitantly, the area of Λ −+ shrinks. This explains the recovery of the attraction in this regime. Finally, we note that if the size of the patch θ p is not large enough, the repulsive effect may not be sufficiently strong in order to change the sign of the radial component of the scaling function. Figure 7 shows the dependence of F (2) r on the scaling variable ∆. Therein, both plots correspond to θ p = 30 • , β 1 = −90 • , and γ 1 = 0. For β 2 = 50 • (see Fig. 7(a)) the radial component of the scaling function is always positive (i.e., there is repulsion) and has a maximum for Θ < 0. Upon increasing the scaled distance ∆ between the particles, the strength of the interaction decreases together with a shift of the position of the maximum towards more negative values of Θ.
The situation is quite different for β 2 = 70 • (see Fig. 7(b)). In this case, for ∆ = 0, the function F (2) r is negative (corresponding to attraction) for all values of Θ. Upon increasing the scaled distance ∆ = D/R between the particles, the value of the scaling function grows for all values of Θ. This leads to a change of sign of F (2) r for Θ < 0; this change sets in for large negative values of Θ and, upon increasing ∆, it propagates towards larger values of Θ. Starting from ∆ ≈ 0.05 the particles repel each other for all Θ < 0 (i.e., in the demixed region of the binary solvent), and, upon a further increase of ∆, we observe the repulsion even for small positive values of Θ. This behavior continues until ∆ ≈ 0.5 is reached. Upon further increase of the scaled distance between the particles, the magnitude of the function F (2) r starts to decay. For ∆ = 0 and Θ < 0 we observe a maximum of the radial component. If ∆ is very large or very small the maximum is very broad and located at a very large negative value of Θ. For Θ > 0 (i.e., in the mixed region of the binary liquid mixture) the radial component of the scaling function has a minimum. Upon increasing ∆ the minimum becomes less deep and moves towards higher values of Θ. This behavior of the minimum is slightly altered for ∆ ≈ 2, where, upon increasing ∆, the minimum moves towards smaller values of Θ and becomes deeper. This anomaly appears in the region where the Derjaguin approximation is expected to be unreliable so that it is physically irrelevant; we refrain from a more detailed discussion of this phenomenon.
The behavior described above can be understood by referring to the properties of the critical Casimir force in the slab geometry. If the scaled distance ∆ between the particles is very small, the region where the surfaces are closest to each other contributes the most to the mutual interaction (due to the prefactor L −3 in Eq. (6)). For β 2 = 50 • the boundary condition at the point of smallest distance is '+−'; thus for small ∆ the particles repel each other. On the other hand, if β 2 = 70 • , the boundary condition at the same point is '++' and the particles attract each other. If ∆ is increased, regions with larger values of become relevant. In the case of β 2 = 50 • , the attractive contribution stemming from the region Λ sm is not sufficiently strong to dominate the radial component of the critical Casimir force. This is not surprising because the magnitude of ϑ op (ω) is larger than the magnitude of ϑ sm (ω) (see Fig. 4). Additionally, below the critical point, ϑ sm (ω) is very small while ϑ op (ω) has a maximum. This is the reason why for Θ < 0 the value of F (2) r grows rapidly upon increasing ∆. Above T c the absolute value of the function ϑ sm has a maximum, while ϑ op is relatively small. This is the reason why, upon increasing ∆, the radial component of the scaling function does not change sign for large values of Θ.
The characteristic narrow plateau of the scaling function F (i) r around Θ = 0, which is visible for all curves in Figs. 6 and 7, is inherited from the scaling functions ϑ sm (ω) and ϑ op (ω) for the slab geometry. For small |ω| these functions behave like A 1 + A 2 |ω| 1/ν [18], where A 1 and A 2 are constants.
In Fig. 8(a), for θ p 45 • the maximum of the radial component of the scaling function of the critical Casimir force is exactly at β 2 = 90 • , in which case the center of the patch on the second particle is in the closest possible position to the first particle. In this position the area of Λ −+ is largest. If θ p > 45 • , for β 2 = 90 • the region Λ ++ becomes a nonempty set and, as a result, the radial force decreases because for elements of Λ ++ one has a negative contribution in Eq. (21). This effect shifts the position of the maximum towards higher values of β 2 . Upon increasing θ p , starting from 45 • , the position of the maximum increases towards 180 • , which is attained for θ p = 90 • . Upon this increase the maximum becomes sharper. In Fig. 8(a), for θ p < 90 • , we observe a singu- 3/2 of type II located at the values β (i) 2,sing. of β 2 for which the patches are tangent (these points are marked by circles in the plot). These singularities are located to the left of the maximum for θ p < 45 • , coincide with the maximum for θ p = 45 • and θ p = 90 • , and for 45 • < θ p < 90 • they are shifted slightly to the right side of the maximum. We note that for θ p = 90 • the radial force has an inverse 'V' shape around β 2 = 180 • . If β 1 = −90 • , as shown in Fig. 8(b), for any size of the patch θ p the radial component of the scaling function is symmetric around β 2 = 90 • (and around β 2 = −90 • ). For β 2 = 90 • the patches are facing each other and (within the Derjaguin approximation) F (2) r does not depend on the size θ p of the patch. Upon decreasing β 2 , a region Λ op with opposing boundary conditions emerges and thus contributing to repulsion so that the radial component grows. Since the magnitude of ϑ op is larger than the magnitude of ϑ sm , we observe a change from attraction to repulsion even for relatively small patches. Upon further decreasing β 2 , the radial component of the scaling −2 function for the force reaches a maximum. This occurs if the overlap of the two patches is small. A further decrease of β 2 reduces F (2) r . This can be understood by noting the fact that moving the patch on the second particle upwards reduces the area of the projection of the patch. Reducing β 2 even further moves the patch to the right hemisphere of the second colloid, where it does not participate in the mutual interaction between particles. This occurs in a region around β 2 = −90 • , where the radial component is constant.
Like in the previous case, we have observed several singularities of F (2) r . Around the minimum at β 2 = 90 • , if the patches face each other, there is a singularity of type I and the scaling function for the force exhibits a 'V' shaped cusp with an opening angle which increases upon increasing θ p . In the case of Janus particles (θ p = 90 • ) the opening angle reaches 180 • . This property of the force follows directly from the Derjaguin approximation: If β 2 is slightly shifted away from 90 • , around the circumference of the patches a region Λ op emerges. If the size of the patches increases, so does within this region. Additionally, we have noticed singularities ∝ β 2 − β (i) 2,sing.

3/2
of type II which occur if the patches are tangent. In Fig. 8 these points are marked with circles.
Finally, we comment on the nature of singularities of the radial component of the scaling function for the critical Casimir force. As we have reported above, the nonanalyticities of F (i) r are similar to those of the scaling function U for the potential rather than to nonanalyticities of the scaling functions for other components of the forces. On one hand this can be explained as a simple consequence of the similarity of the formulae from which F (i) r and U are calculated (see Eqs. (21) and (22)). On the other hand, all singularities of U discussed in Sec. IV D manifest themselves upon changing the rotational configuration Ω of the system. In contrast to all other components of the scaling function for the forces, the calculation of F (i) r does not require derivatives with respect to angles (see Eq. (A8)); therefore the character of the nonanalyticities remains unchanged.

B. Critical Casimir potential
In this subsection we present our results for the scaling function U of the critical Casimir potential (see Eq. (16)) calculated within the Derjaguin approximation.
In Fig. 9 the scaling function U as a function of Θ is shown for various values of β 2 with all other parameters fixed (∆ = 1, β 1 = γ 1 = 0, θ p = 30 • ). If β 2 = 0 the patches on both spheres are in the topmost position and, within the Derjaguin approximation, the boundary conditions are the same everywhere. Thus for any value of Θ the function U attains its smallest value via this configuration. Upon increasing β 2 , with all the other pa- rameters fixed, a region Λ −+ emerges (i.e., it becomes a nonempty set) and, as a result, U increases. This growth is more pronounced for Θ < 0 (where the scaling function ϕ op relevant for the slab geometry is the largest), and for β 2 > 0 the scaling function as a function of Θ has a maximum at Θ < 0. Correspondingly, upon increasing β 2 the minimum, visible in Fig. 9 for Θ > 0, becomes less deep, more shallow, and, eventually, somewhere between β 2 = 60 • and β 2 = 90 • , it disappears. The scaling function reaches its maximum for β 2 = 90 • . In this special configuration the area of Λ op is largest and, concomitantly, the mean value of the distance , between pairs of points on the surfaces projected onto the region, is smallest. A further increase of β 2 reduces U and drives it negative again. For β 2 = 180 • the potential is only slightly larger than the lower bound corresponding to β 2 = 0. Finally, we note that the reported behavior of the scaling function U strongly depends on the values of the fixed parameters. If the size θ p of the patch is small and the scaled distance ∆ is sufficiently large, U can be negative for all values of Θ and β 2 .
In Fig. 10 we present the plot of the scaling function for the critical Casimir potential U for fixed values of ∆ = Θ = 0.5, γ 1 = 0, and θ p = 60 • . The value of the function for the whole ranges of β 1 and β 2 is presented by resorting to the color code. We note that, due to reflection symmetry, the scaling function in Fig. 10 is invariant under the transformations The minimum of the function U, with U min ≈ −1.1552, is attained (within the Derjaguin approximation) if the boundary conditions are '++' or '−−' for every point of Λ. This occurs if the patches are in the mirror-symmetric configuration (along the diagonal line β 1 = −β 2 ) or if both patches are on those hemispheres which do not participate in the interaction (rectangle θ p < β 1 < 180 • − θ p and −180 • + θ p < β 2 < −θ p ) (see Fig. 10(b)). In the latter case, the size of the rectangular region depends on the size θ p of the patch in that it shrinks upon increasing θ p and disappears completely for θ p = 90 • . As can be inferred from the cuts in Figs. 10(c) and 10(d), crossing the region in which U is minimal always reveals nonanalyticities: 'V' shaped cusps (i.e., a singularity of type I) for the line β 1 = −β 2 (outside the rectangle) and of type III for the edges of the rectangle. We note that all properties reported here have been obtained within the Derjaguin approximation and we expect them to be absent beyond this approximation.
The maximal value of the scaling function U, with U max ≈ 4.3949, is attained at four isolated points: (β 1 ≈ 19 • , β 2 ≈ 99 • ) and at three points which can be generated from this first one by exploiting the symmetries described in Eq. (35). At these four special points various influences counter each other. Changing β 1 or β 2 increases the area of Λ ++ ; it moves parts of the patch to the far side hemispheres, where it does not participate in the interaction; or it moves the patch up or down which reduces the area of Λ op . Unlike the minima, the precise positions of the maxima depend sensitively on the choices of ∆, Θ, and θ p .
The dependence of U on β 1 for fixed β 2 = 90 • is presented in Fig. 10(c). This plot is symmetric around β 1 = −90 • , where it has a minimum and a nonanalyticity of type I. In this special configuration the two patches are facing each other. Upon increasing β 1 from −90 • , there emerges a region Λ op and, as a result, U increases. For β 1 > −60 • , increasing β 1 moves some part of the patch on the first particle to the left hemisphere, where (within the Derjaguin approximation) it is not participating in the interaction. At first, this effect is not very pronounced, but, upon increasing β 1 from −60 • , it becomes more relevant. First, the growth rate of U is reduced and finally, for β 1 ≈ 28 • , the effect becomes dominant and the scaling function starts to decrease. If 60 • < β 1 < 120 • , the patch on the first particle is located fully on the left hemisphere and the function U becomes constant. Since U is symmetric around β 1 = −90 • , a further increase of β 1 does not yield any new phenomena.
If β 2 = −90 • , the patch on the second colloid is fully located on the right hemisphere, and it does not influence the interaction. The cut of the scaling function U in this special case is presented in Fig. 10(d). There is a maximum for β 1 = −90 • , which is attained if the patch on the first particle is in the closest possible position to the second particle. In this configuration, changing β 1 increases the distance between the points on the patch of the first particle and the points on the surface of the second sphere. Simultaneously, it reduces the area of Λ +− and thus the repulsion; accordingly the scaling function for the potential decreases. This decrease continues for 60 • < β 1 < 120 • until the patch on the first particle is fully located on the left hemisphere; within this range of values for β 1 the function U is the same as in the case of homogeneous spheres and thus has a flat minimum.
Next, we study the dependence of the critical Casimir potential on the distance between the particles. In Fig. 11(a) we present the plot of the potential as a function of the surface-to-surface distance D at a fixed supercritical temperature at which ξ b (T ) = R for two particles with patches of size θ p = 30 • , both in the leftmost position (i.e., β 1 = β 2 = 90 • and γ 1 = 0). Accordingly, the patch on the first particle is located on the left hemisphere and, thus, it does not influence the interaction. The potential has been calculated from the scaling func-tion U by using Eq. (16).
Within the Derjaguin approximation, the distance between pairs of interacting points varies between D and D+2R. Since the potential in the slab geometry increases rapidly upon decreasing the distance between walls (see Eq. (9)), for D R the closest pairs of points (for which the boundary conditions are '−+') dominate the whole interaction. This explains why for small D the potential is positive and diverges for D → 0. If D R, the effect of different distances for different pairs of points is negligible, and the sign of the potential is determined by the size of the patch. If the patch is sufficiently small, the area of Λ −+ is not large enough to facilitate a change of sign, so that for D → ∞ the potential approaches 0 from below. This leads to the conclusion that for small values of θ p the potential as a function of D has a minimum at a certain distance D 0 . We expect that this holds also beyond the Derjaguin approximation.
In Fig. 11(b), we plot the dependence of the position of the minimum D 0 of the critical Casimir potential U c (see Eq. (16) We note that the minimum of the potential is observed only for shifting the second particle in radial direction. The configuration corresponding to the minimum is not stable with respect to rotations, and thus for D = D 0 the system occupies a saddle point.

C. Non-radial components of the critical Casimir force
In this subsection, we present our results for the critical Casimir force and potential, when the second particle is moved in z direction (i.e., in a direction perpendicular to the line connecting the centers of the spheres in the initial position). For such a setup, even though the spheres are not rotated in the laboratory reference frame, their relative orientations change upon moving the particles. In the reference frame associated with Ω * (see Sec. II) both spheres rotate around the x axis and, at the same time, the distance between them increases. This way, the effect of rotations can be studied without actually performing any rotations, which makes the setup a good candidate for possible lattice-based simulations. We note that all calculations presented here have been carried out for a system in equilibrium; we do not consider any dynamic effects.
In Fig. 12, we plot the results for the potential and the z component of the critical Casimir force acting on the second particle for various configurations and sizes of the patches. The results have been calculated from the scaling functions by using Eqs. (14a) and (16). In all plots r x = 0, r y = 2.5R, and the temperature is chosen such that ξ b (T ) = R and T > T c .
We first consider spheres without any patches. In Fig. 12(a), the potential and the force for the case of homogeneous particles is plotted. The potential exhibits a minimum at r z = 0, which is the position where the colloids are in the closest possible position. Upon increasing r z , the negative potential is gradually increasing towards 0. If r z = 0, there is no z component of the critical Casimir force. Upon increasing r z a force appears which acts in the negative z direction. For small, increasing values of r z the absolute value of the z component of the force increases and reaches its maximum value for r z ≈ 0.9R, and, upon a further increase of r z , it decays to 0. This behavior can be understood by noting that the force acting between two homogeneous spheres is radial. For small values of r z the total force is large, but it is almost perpendicular to the z direction, so that the z component is small. For larger values of r z , the angle between radial and z direction decreases, but simultaneously the magnitude of the force decreases; the interplay between these two effects produces the maximum of the absolute value of the z component of the force.
The critical Casimir potential and the z component of the force in the presence of patches, initially located directly opposite to each other (β 1 = −90 • and β 2 = For all plots the temperature T > Tc is chosen such that ξ b (T ) = R, and the initial distance between the colloids is ry = 2.5 R (i.e., initial surface-to-surface distance D = 0.5 R). In the columns (b) and (c) the interaction exhibits a singularity of type I at rz = 0 and of type II for rz = ry tan θp; the latter is beyond the range of the plots except for the case of θp = 30 • in column (b). We note that all the results have been obtained for systems in equilibrium; in this study we do not consider any dynamic effects. 90 • ), are plotted in Fig. 12(b). If r z = 0, within the Derjaguin approximation, all pairs of interacting points have the same boundary conditions, and thus the value of the potential does not depend on the size of the patch. For θ p > 0, in this configuration, there is a singularity of type I, i.e., the potential has a 'V' shaped cusp around r z = 0; and the force is discontinuous, in the limit r z → 0 + it is negative (acting downwards). In this limit, the magnitude of the force depends on the size of the patch.
It is very small for θ p → 0 and for θ p → 90 • , and maximal for θ p ≈ 20 • (this behavior is not presented in the plot). Due to symmetry F c,z = 0 for r z = 0, which differs from the nonzero values in the limit r z → 0 + (see the bottom panel in Fig. 12(b)). This observation is in full agreement with the properties of the Derjaguin approximation: the jump of the force at r z = 0 is generated by the region Λ op which emerges when the second colloid is moved slightly.
For small shifts, this region has a ring-like shape around the circumference of the projections of the patches onto the projection plane. If θ p is small, the circumference of the patch is very small, while for θ p 90 • the region Λ op is located close to the boundary of Λ, where the surfaces of the spheres are almost perpendicular to the projection plane. In both cases the area of the region Λ op cannot be large, and thus, the jump of the force at r z = 0 is small. Upon increasing r z from 0, at first the magnitude of the z component of the force increases. This is the same effect as the one occurring for homogeneous spheres: the radial component dominates the force, and the increase of r z reduces the projection angle. Upon a further increase of r z , the influence of the patches becomes visible. The repulsion stemming from the region Λ op reduces the magnitude of the z component of the critical Casimir force. If the size θ p of the patch is sufficiently large, a further increase of r z changes the sign of F (2) c,z so that the force starts to act upwards. If r z /r y = tan θ p , i.e., if the projections of the edges of the patches on the two spheres are tangent, we observe a singularity of type II. (In Fig. 12(b) this occurs within the range of the plot only for θ p = 30 • .) For moderate values of θ p , the cusp with infinite slope in the plot of the force, which is characteristic for this type of singularity, coincides with the global maximum of the z component of the force; for larger values of θ p , the maximum develops on the left side of the cusp. Upon a further increase of r z , the force decays to zero as the distance between the patches grows.
In the third considered configuration both patches are in the uppermost position (β 1 = β 2 = 0). The corresponding critical Casimir potential and the z component of the force are plotted in Fig. 12(c). Like in the previous case, for r z = 0 the value of the potential does not depend on θ p and there is a singularity of type I -the potential has a 'V' shape around r z = 0, and the force jumps from positive values for r z < 0 to negative ones for r z > 0. In contrast to the previous case, here the absolute value of the force for r z → 0 + grows monotonically upon increasing the size of the patch. This observation can be understood from the fact that, if r z > 0 is very small, the typical distance between the points in the region Λ op decreases upon increasing the patch size (for θ p < 90 • ).
For small sizes of the patches the force is negative (i.e., it acts downward) for all values of r z . If θ p is sufficiently large, there appears a maximum where the force is positive (i.e., it acts upwards). Upon increasing θ p , at first, the maximum is very broad and it is located at very large values of r z . A further increase of θ p moves the maximum towards smaller values of r z and makes it more pronounced. In Fig. 12(c) the maximum is located within the range of the plot only for θ p = 90 • . This behavior can be understood on the basis of our formulae. For small values of r z , the particles attract each other because the area of the region Λ sm is much larger than the area of the region Λ op . Upon increasing r z , the area of the former region decreases, whereas the area of the latter region increases. If the patches are sufficiently large, this eventually produces the repulsion between the colloids, and the value of the z component of the force becomes positive.

D. Critical Casimir torque
In this subsection we discuss the critical Casimir torque acting on the colloidal particles. We calculate the scaling function for the torques T (1) and T (2) by using the numerical derivatives of U and by deriving the torque via Eq. (14b).
We first study the special case γ 1 = 0, in which the system is invariant under mirror reflection at the yz plane. This symmetry implies that the torque can only act in x direction. Moreover, the system possesses the additional so that it is sufficient to study the torque T c,x acting on the first particle.
In Fig. 13 the torque for θ p = 60 • , γ = 0, and β 2 = −90 • , with fixed distance D, radius R, and temperature T , is plotted as function of β 1 as a black curve. In this configuration, the patch on the second particle is located on the right hemisphere and (within the Derjaguin approximation) it is not participating in the interaction; the torque is the same as if the second particle had no patch. We note that in this case the torque is proportional (with a negative coefficient) to the derivative of the scaling function for the potential U with respect to β 1 (see Eq. (A8e)); this function U is plotted in Fig. 10(d).
For −90 • < β 1 < 60 • the torque is positive, which means that in Fig. 10(a) it acts as to rotate the first particle anticlockwise. The maximum is located slightly below β 1 = −30 • . These properties are not surprising because increasing β 1 (starting from −90 • ) adds points to the region Λ +− with large values of and removes those with smaller values of , which reduces the interaction free energy. The torque is largest close to the point, where the boundary conditions for the closest pair of points are changing; the neighborhood of these points is expected to dominate the interaction for small values of D/R. For 60 • < β 1 < 120 • the patch on the first particle is fully located on the left hemisphere, and thus in this interval there is no torque acting in the system.
The situation changes strongly for β 2 = 90 • , i.e., when the patch on the second particle is moved to the leftmost position. In this case, the x component of the critical Casimir torque is plotted in Fig. 13 as a red curve, and the scaling function for the underlying potential is presented in Fig. 10(c). If β 1 = −90 • , the patches face each other and there is no torque. Upon a slight increase of β 1 , the torque jumps from zero to a certain negative value (this jump is caused by a singularity of type I). Negative values of the torque mean that in Fig. 10(a) it acts as to rotate the first particle clockwise. Upon a further increase of β 1 , the torque decreases, reaches its minimum for β 1 slightly below 30 • , increases, changes sign, and has a positive maximum for β 1 = 30 • , where there is a singularity of type II with a characteristic cusp with infinite slope. For β 1 > 30 • (and β 1 < 150 • ), the torques for β 2 = 90 • and β 2 = −90 • (red and black curve in Fig. 13, respectively) are identical.
The observed behavior follows directly from the properties of the Derjaguin approximation. For β 1 = −90 • , the boundary conditions are the same everywhere and the potential has a minimum. Upon increasing β 1 the region Λ op becomes a nonempty set, and thus the free energy increases. The minimum of the torque is located close to the point where the boundary conditions for the closest points change. Upon a further increase of β 1 , a second effect appears according to which the patch on the first particle is moved to the left hemisphere and does not participate in the interaction. This increases the torque and eventually leads to the change of its sign. The maximum of the torque for β 1 = 30 • occurs at that position for which the region Λ ++ consists of a single point. For this configuration the area of Λ op is maximal. Finally, if β 1 > 30 • , there is no overlap between patches, the region Λ −+ does not change with β 1 , and thus the integral in Eq. (22) over this region does not depend on β 1 and does not contribute the derivative of U with respect to β 1 . Simultaneously, the integral over the region Λ +− is exactly the same for β 2 = −90 • and β 2 = 90 • . This explains why the x components of the torques are identical in these two cases (see β 1 > 30 • in Fig. 13).
We now discuss the case of γ 1 = 0. For simplicity, we fix β 1 = β 2 = 0 and ∆ = Θ = 0.5, and vary γ 1 from 0 • to 180 • . In this case all components of the critical Casimir torque can now be nonzero. The dependence of T (1) c,y and U c on γ 1 is presented in Fig. 14.
For γ 1 = 0, the patches on both particles face each other and, for any size θ p of the patch, one has Λ = Λ sm ; in this case and within the Derjaguin approximation, the critical Casimir potential attains its lowest possible value. Increasing the value of γ 1 renders the set Λ op nonempty, which increases U c and generates a y component of the torque acting against the rotation (T (1) c,y < 0). This situation changes if γ 1 exceeds the value of 2 θ p . In this case, there is no region Λ ++ anymore (i.e., this set is empty) and varying γ 1 shifts the region Λ +− , without changing its shape, in such a way that the integral in Eq. (22) remains the same. As a result, for 360 • − 2 θ p > γ 1 > 2 θ p , the potential is constant and there is no torque in y direction. In the special case of Janus particles (θ p = 90 • ), the edges of the patches form two diameters of the circle Λ on the projection plane, and γ 1 is the angle between them. Since depends only on the distance from the center of the circle Λ on the projection plane, the integral in Eq. (22) is linear in the area of region Λ op . As a result, the torque in y direction is constant and negative for 0 < γ 1 < 180 • , changes sign for γ 1 = 0 and γ 1 = 180 • , and is constant and positive for 180 • < γ 1 < 360 • .
We note that in the case of θ p < 90 • , for γ 1 = 2 θ p there is an unusual nonanalyticity. There is a jump in the second derivative of the potential with respect to γ 1 from a certain negative value for γ 1 < 2 θ p to 0 for γ 1 > 2 θ p . For γ 1 = 2 θ p , there is just one point in Λ ++ , but it is located at the boundary of Λ. This explains why the observed nonanalyticity cannot be classified into any of the three types described in Sec. IV D; we introduce a new type IV for such rare and untypical singularities (for further details see Appendix C 4).

E. Comparison with experimental results
One of the ways to test the validity of our calculation is to compare the results, obtained by using the Derjaguin approximation, with available experimental data. Up to our knowledge, most of the experimental papers focus on the problem of pattern formation by resorting to nonspherical colloids. Therefore, only a qualitative comparison is possible. In this subsection, we perform a simple check of our results with the pair potential measurements reported in Ref. [60].
In this experiment a binary liquid mixture of heavy water and 3-methylpyridine close to its lower critical demixing point is used as a solvent. The colloidal particles, immersed in the mixture, have an oblong shape with the surface preferring 3-methylpyridine on the two tips and heavy water in the middle of the particle (see Fig. 15(a)). For such a setup, the formation of superstructures by the colloids is investigated. Below the lower critical temperature (i.e., in the mixed region), depending on the concentration of the two components of the mixture, the colloids can align in two different ways: as needle-like chains, for which only the tips of two neighboring particles are in contact, and as fence-like chains in which the particles are parallel to each other (i.e., both the tips and the middle parts of neighboring particles are close to each other). Additionally, upon approaching the critical point a collapse of both types of chains is reported. Here, we are especially interested in the experimentally determined pair potential and bending stiffness of neighboring particles in needle-like chains.
Certainly, the actual particles as described above cannot be fully modeled within the present framework. Nonetheless, in order to proceed, we assume that the colloids can be approximated by a sphere of a radius R with two chemical inhomogeneities, i.e., two circular patches of size θ p in antipodal configuration (see, e.g., Fig. 15(b)). The optimal radius R and size θ p of the patches are yet to be determined.
In Ref. [60] the effective pair potential of the interaction between two neighboring colloids in needle-like chains is estimated and inferred from the observed spatial distribution of the colloids. The radial component of the interaction is dominated by the competition between critical Casimir attraction and electrostatic repulsion. This effect has already been studied in the context of homogeneous spheres [87]. Therefore we can focus on the angular dependence of the effective pair potential, which is generated by the chemical inhomogeneity of the surface of the colloids and by their nonspherical shape. The dependence of the pair potential on the so-called bending angle (which is an analogue of the angle λ in Fig. 15(c)) turns out to be parabolic, and the bending stiffness (which is proportional to the second derivative of the potential with respect to the bending angle) decreases upon approaching the critical point (see Fig. 2(c) in Ref. [60]).
In order to facilitate a comparison with our present results, it is necessary to determine the values of the relevant parameters in the experiment; most of them can be found in Ref. [60]: The critical temperature is T c ≈ 38.55 • C and the bulk correlation length is where t = (T c − T ) /T c is the reduced temperature, the critical exponent ν is given by Eq. (8), and ξ + 0 ≈ 1.6nm has been estimated in the course of the experiment. (For a molecular liquid the value of ξ + 0 is rather large and comparable to the one in nonionic micelle solutions [60,88].) Since the mean value of the surface-to-surface distance D between the colloids is not reported (most probably it was fluctuating during the measurement), we choose D = 0.2µm, for which the total radial potential (electrostatic plus critical Casimir) is found to exhibit a minimum (see Fig. 2(a) in Ref. [60]). Because in the experiment the particles are not spherical, we estimate the effective radius to be R = 1.15µm. This value is certainly smaller than the actual size of the particles; however, for this effective radius the curvatures of the surface at the closest points of the two particles -within the Derjaguin approximation this region provides the most important contribution to the interaction -are the same as in the experiment. The only parameter which cannot be reliably estimated is the effective size θ p of the patches; we have checked that in order to obtain the same order of magnitude of the potential as in Ref. [60], θ p should be roughly between 30 • and 40 • .
We assume that the bending potential is generated solely by the critical Casimir interaction, neglecting a possibly inhomogeneous distribution of the surface charge. Accordingly, fixing the surface-to-surface distance to its experimental equilibrium value allows us in the following to not consider the electrostatic repulsion.
In order to calculate the effective potential, we put two colloids in the configuration shown in Fig. 15(c). If λ < 90 • − θ p , one of the patches on each sphere is located on the distant hemisphere and thus, within the Derjaguin approximation, it does not participate in the mutual interaction. This implies that the potential can be calculated without introducing any modifications into our single patch framework (as described in Sec. II). It is sufficient to consider spherical particles with a single patch of size θ p in a configuration defined by β * 1 = −90 • − λ, β * 2 = 90 • , and γ * 1 = 0. We note that the measurements were done in the vicinity of the lower critical demixing point, for which the mixed phase is observed below the critical temperature. Therefore, here the amplitude ξ + 0 of the singular part of the correlation length is associated with T < T c .
In Fig. 15(d) the dependence of the critical Casimir potential on the bending angle λ is plotted for various temperatures below T c (i.e., in the disordered phase) for a patch size θ p = 30 • . This is our analogue of the experimental data presented in Fig. 2(c) in Ref. [60]. In both systems the potential becomes stronger upon increasing the bending angle λ; shifting the temperature away from T c reduces the growth of the bending potential as function of λ. This implies that λ = 0 is the stable configuration with respect to rotations, and that the bending stiffness grows upon lowering the temperature below the critical point (i.e., into the disordered phase). In contrast to the experimental results, the calculated potential is not parabolic around λ = 0. Instead, there is a singularity of type I and the function displays a 'V' shape. As discussed in Sec. V, the latter is an artifact of the Derjaguin approximation. Therefore, in our case the limit of the derivative lim λ→0 + ∂U c /∂λ is a suitable expression for the bending stiffness; we present a plot of this quantity in Fig. 15(e). As in the experiment (see the inset in Fig. 2(c) in Ref. [60]), the bending stiffness decreases upon shifting the temperature away from the critical point. We note that, because of the singularity at λ = 0, our definition of the bending stiffness differs significantly from the bending stiffness used in Ref. [60] and therefore a detailed quantitative comparison is not possible.
The analysis presented above allows us to conclude that for the system under consideration the calculated properties are in qualitative agreement with the experimental data. The observed differences can be traced back to artifacts of the Derjaguin approximation.

VII. CONCLUSIONS
We have studied the critical Casimir interaction between two identical patchy colloidal particles immersed in a binary liquid mixture close to demixing. The surface of the spherical colloids is chemically inhomogeneous in that it prefers one component of the binary solvent everywhere, except for a circular patch of angular size θ p , where there is an affinity to the second component of the binary mixture (see Fig. 2). When the mixture is close to its critical demixing point, the critical Casimir interaction between colloids emerges. It is anisotropic due to the inhomogeneous surface of the colloids and can be tuned by varying the thermodynamic parameters of the binary liquid mixture.
We have introduced the angles which, together with the surface-to-surface distance D and the radius R of the spheres, describe the relative configuration of the two particles (see Fig. 1). Using these parameters, we have formulated the scaling laws for the critical Casimir potential, force, and torque (see Eqs. (14) and (16)). We have derived the formulae which allow one to calculate the scaling functions for the force and the torque from that of the potential (see Eq. (A8)). This part of our study is very general; it can be applied to particles of arbitrary shapes and chemical surface patterns.
We have used the following method to calculate the aforementioned universal scaling functions: First, the scaling function for the critical Casimir interaction potential between two colloids has been calculated by applying the Derjaguin approximation. Within this approximation, the scaling function for a nonelemental geometry is expressed in terms of an integral over the relevant scaling functions for the system in the slab geometry (see Eq. (22)). Second, by numerical differentiation we have derived the universal scaling functions for all components of the forces and torques acting on both particles (see Sec. IV). The validity and accuracy of this approach have been discussed in Sec. V. Among the artifacts of the Derjaguin approximation we have identified various nonanalyticities of the scaling functions constructed upon this approximation (see Sec. IV D); their detailed analysis is a prerequisite for accurate numerical calculations.
In Sec. VI, we have applied the above method of calculation in order to study the critical Casimir potential, force, and torque acting between two colloids for various configurations and values of the relevant parameters. The analysis has revealed a complex interaction which can be tuned in a controlled way: The radial component of the force can be either repulsive or attractive and there can be a stable (i.e., with respect to radial shifts) position for an arbitrary value of D (see Fig. 11). If in the initial configuration one of the colloidal particles is shifted perpendicular to the radial direction, the magnitude of the force counteracting the shift can be tuned by changing the angular size θ p of the patches and the orientation of the particles. Moreover, for large distances the force can change sign (see Fig. 12). By varying the size of the patch a similar tuning has been observed for the torque (see Fig. 13).
Finally, we have performed a comparison of our results for the angular dependence of the critical Casimir interaction with available experimental data [60]. Even though the dependence of the free energy on the bending angle differs quantitatively between theory and experiment (due to the singularities introduced by the Derjaguin approximation), the temperature dependence of the bending stiffness shows qualitative agreement (see Sec. VI E).
In order to gain further insight into critical Casimir interactions between colloids with inhomogeneous surfaces, it is necessary to study analogous systems by using a variety of different techniques, such as Monte Carlo simulations or mean field theory. The comparison of such independent results with those presented here would allow one to determine the actual interaction more accurately. Moreover, the approach presented here allows one to extend our techniques to more complicated shapes of patches or to nonspherical colloids.
The results reported in this paper can also be used to perform molecular dynamics simulations in order to study clustering phenomena for patchy particles. Unfortunately, our numerical routines for calculating forces and torques are not fast enough to be directly usable for this purpose. Nevertheless, with additional efforts it is possible to overcome this problem: The scaling function for the potential U can be decomposed into its singular and its nonsingular part. Using the results from Appendix C, it is possible to propose exact formulae for the singular part, whereas the nonsingular part can be expanded into a series of appropriate orthogonal basis functions. The details of this procedure are beyond the scope of the present study. and cos β * 2, = cos β 2 .
For β 2 = ±90 • , the only solution of the above equations is ρ = (tan β 2 ) /r and By comparing various elements of the matrices in Eq. (A4), a straightforward calculation yields and where we have additionally assumed that β 1 = ±90 • . Inserting Eqs. (A3) and (A6) into Eq. (A1) leads to Finally, using Eqs. (14) and (16), after some algebra one obtains the relation between the scaling function for the potential and the scaling function for one of the components of force: The calculation above can be repeated for all components of the forces and torques. Here, we skip the details and provide only the final results: and A straightforward calculation shows that Eq. (A8) does satisfy the relations in Eq. (15). Additionally, the radial component of the force, given by Eq. (A8b), has been checked numerically to agree with the results obtained from the Derjaguin approximation for the force via Eq. (21).
The results for the force and the torque (Eq. (A8)) are not valid in the special cases β 1 = ±90 • or β 2 = ±90 • , for which certain coefficients in Eq. (A8) diverge. A careful investigation shows that in these cases the rotations around the z and y axes are not independent. Moreover, an infinitesimal rotation can lead to a non-infinitesimal change of Ω * . Since in most cases it is sufficient to consider β 1 and β 2 to be close but not equal to ±90 • , we refrain from reporting the formulae for the force and the torque in these special cases. In this appendix we discuss the formulae for the scaling functions ϑ sm and ϑ op for the critical Casimir force in the slab geometry with same and opposite boundary conditions, respectively. In spatial dimension d = 3 these functions are obtained by interpolating and extrapolating the data calculated numerically by using Monte Carlo simulations [73] in such a way, that all the known properties of the scaling functions are fulfilled.
In our analysis we are using the formulae obtained by A. Gambassi et al. [25,89]. They were successfully applied in several distinct studies [26,34,51] but, the fitted functions have not yet been documented in the literature.
The fit uses a scaling law for the critical Casimir force which differs slightly from Eq. (6): where P s is the scaling function and where the scaling variable x = L/ξ + 0 1/ν t is linear in the reduced temperature t. Thus the relation between the scaling functions P s and ϑ s is given by ϑ s (ω) = P s (ω 1/ν ) for ω 0, where A ξ is the ratio of the critical amplitudes of the bulk correlation length (see Eq. (8)). The fitted shapes of the scaling functions are as follows: where the numerical coefficients have been chosen to reproduce as good as possible the data marked as '(i)' in Figs. 9 and 10 in Ref. [73]. In the above formulae the critical index ν has been assumed to be 0.63. The accuracy of this fit is discussed in Sec. IV A. We note that the above forms of P sm and P op exhibit the correct behaviors [18] for x → ∞, x → −∞, and |x| 1; the functions are continuous but have a slight jump of the derivative at the gluing points.
The formula in Eq. (22) for the scaling function U can be transformed into U (∆, Θ, Ω * ) = Λ dθdφ J sm (θ, φ) + Λop dθdφ δJ (θ, φ) , (C2) where we have introduced δJ (θ, φ) = J op (θ, φ) − J sm (θ, φ) . (C3) The sets Λ and Λ op , as defined in Sec. III C, are parametrized by the spherical coordinates 0 θ π and 0 φ π on the first particle. In Eq. (C2), the first term on the right-hand side is equal to the scaling function for spheres without patches. Thus it depends neither on the relative configuration Ω * nor on the size θ p of the patches. Therefore, the nonanalyticity of the type I is produced by the second term in Eq. (C2).
If the patches on two spheres are exactly in a mirrorsymmetric configuration, Λ op is an empty set. If one of the spheres is slightly rotated, Λ op has a ring-like shape of variable thickness, following the circumferences of the projection of the patches onto the projection plane. For a more quantitative study, we assume that initially the patches on the spheres are in a mirror-symmetric configuration, opposite to each other, and that one of the spheres is rotated (around a certain axis) by a small angle σ. We denote as Λ σ op the region on the projection plane corresponding to opposing boundary conditions in the final configuration.
In the process of rotation, the circumference of the projection of the patch on the rotated sphere passes all points of Λ σ op . The circumference can be smoothly parametrized by θ =θ (u, v) and φ =φ (u, v), where u 0 u u 1 parametrizes the points of the circumference of the patch and 0 v σ measures the rotation angle. We note that before the rotation the circumferences of both patches are given byθ (u, 0) andφ (u, 0), while after the rotation the circumference of the patch on the rotating sphere is given byθ (u, σ) andφ (u, σ).
The second integral in Eq. (C2) can now be rewritten as where C 1 and C 2 are parameters which do not depend on σ. This explains the observed characteristic 'V' shape of the potential, and proves that the derivative of U with respect to σ does not exist for σ = 0 and its left and right limits have the same absolute value but different signs.

Singularity of type II
The singularity of type II appears if on the projection plane the circumferences of the patches are tangent in a single point. In order to describe these situations, we have introduced the overlap angle ζ II overlap (see Eq. (29)). In order to study the scaling function for the potential in more detail, we decompose the integral in Eq. (22): where we have introduced the sets Λ 1 = Λ ++ ∪ Λ +− and Λ 2 = Λ −+ ∪ Λ ++ which are the images on the projection plane of the patch on the first and the second particle, respectively. In the above equation, the integrand J sm is given by Eq. (C1) and δJ by Eq. (C3). The variables 0 θ π and 0 φ π are the spherical coordinates on the first particle and parametrize the set Λ on the projection plane.
The first term on right-hand side of Eq. (C6) describes the interaction of two spheres without patches, the second and the third term depend only on the position of the patch on one of the spheres. Since the nonanalyticity of type II depends on the position of both patches (via the overlap angle ζ II overlap ), it can only appear in the last term on the right-hand side of Eq. (C6): If ζ II overlap < 0, the patches do not overlap, Λ ++ is the empty set, and U II = 0. For 0 < ζ II overlap 1, the area of overlap is small and the integral in Eq. (C7) can be approximated by where θ 0 and φ 0 are the coordinates on the projection plane of the point at which the projections of both patches are tangent for ζ II overlap = 0. After some calculation, one obtains U II = − 8 3 sin θ 0 tan θ p δJ (θ 0 , φ 0 ) ζ II which is in full agreement with Eq. (31). In the above equation, Θ H (x) is the Heaviside step function. C 3 , C 4 , and C 5 are parameters which do not depend on ζ II overlap . The values of C 3 and C 4 are determined by first three terms in Eq. (C6), and is produced solely by the last term in Eq. (C6).

Singularity of type III
The singularity of type III occurs when the projection of the edge of one of the patches is tangent to the border of the circle Λ. This can happen in two different instances: (i) The patch is fully located on the hemisphere that is not participating in the interaction and, upon rotation, a part of it is moving towards the interacting hemisphere. (ii) The patch is fully located on the interacting hemisphere and the rotation moves a part of it to the other hemisphere.
Like in the previous case of the nonanalyticity of type II, the area of the relevant part of the patch is proportional to ζ III overlap,i 3/2 , where the overlap angle ζ III overlap,i is defined in Eq. (32). It is possible to derive the ex-pansion of the singular term in the interaction potential around the singularity. The result is very similar to Eq. (C9). Since for any of the limits θ, φ → 0, π we have δJ (θ, φ) / sin θ → 0, unlike the previous case, the leading order term in the expansion vanishes and, therefore, the next term, proportional to ζ III overlap,i 5/2 , dominates. Since this nonanalyticity is very weak, for our purposes a more detailed analysis is not necessary.

Singularities of type IV
All the other possible nonanalyticities require a special configuration according to which all three curves in the projection plane (i.e., the projection of the edges of the two patches and the circumference of the circle Λ) must meet in one point. In order to reach such a configuration, one has to tune all three angles, which makes these nonanalyticities rare. Moreover, the type of nonanalyticity detectable in the plots depends on how the parameters are changed. We note that, unlike for singularities of other types, the curves do not have to be tangent.
One of the possible situations, in which singularities of type IV become manifest, is presented in Fig. 14. In this case, the projections of the edges of the two patches are tangent at a point which is on the circumference of the circle Λ. As discussed in Sec. VI D, in this case the scaling function for the potential as a function of γ 1 is continuous, has a continuous first derivative, but exhibits a jump in the second derivative.
A detailed analysis of all possible cases is quite complicated and, because all these nonanalyticities are artifacts of the Derjaguin approximation, goes beyond the scope of the present study.