Thermodynamic Casimir Forces between a Sphere and a Plate: Monte Carlo Simulation of a Spin Model

Thermodynamic Casimir Forces between a Sphere and a Plate: Monte Carlo Simulation of a Spin Model

Martin Hasenbusch Martin.Hasenbusch@physik.hu-berlin.de Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany
August 17, 2019
Abstract

We study the thermodynamic Casimir force between a spherical object and a plate. We consider the bulk universality class of the three-dimensional Ising model, which is relevant for experiments on binary mixtures. To this end, we simulate the improved Blume-Capel model. Following Hucht, we compute the force by integrating energy differences over the inverse temperature. We demonstrate that these energy differences can be computed efficiently by using a particular cluster algorithm. Our numerical results for strongly symmetry breaking boundary conditions are compared with the Derjaguin approximation for small distances and the small sphere expansion for large distances between the sphere and the plate.

pacs:
05.50.+q, 05.70.Jk, 05.10.Ln, 68.35.Rh

I Introduction

In 1978 Fisher and de Gennes FiGe78 realized that when thermal fluctuations are restricted by the finite extension of the system, a force acts on the boundaries. Since this effect is similar to the Casimir effect, where the restriction of quantum fluctuations induces a force, it is called thermodynamic Casimir effect. Recently this force could be detected for various experimental systems and quantitative predictions could be obtained from Monte Carlo simulations of spin models Ga09 . The thermodynamic Casimir force is experimentally and maybe technologically interesting since it depends, in contrast to other forces, strongly on the temperature. Hence, since the temperature can be easily changed in experiments, the thermodynamic Casimir force can be easily manipulated. For a discussion of the thermodynamic Casimir force in the general context of nanoscale science see ref. Frenchetal .

Experiments on thin films of He and mixtures of He and He near the -transition GaCh99 ; GaCh02 ; GaScGaCh06 and on thin films of binary mixtures near the mixing-demixing transition FuYaPe05 ; RaBoJa07 have been performed. Theoretically these experiments are well described by a plate-plate geometry. For this geometry the thermodynamic Casimir force per area follows the finite size scaling form Krech

(1)

where is the thickness of the film, is the reduced temperature, is the dimension of the bulk system and is a universal function that depends on the bulk and boundary universality classes that characterize the film. The correlation length in the high temperature phase behaves as , where is the amplitude of the correlation length in the high temperature phase and is the critical exponent of the correlation length.

Field theoretic methods do not allow to compute for three-dimensional systems in the full range of the scaling variable Diehl86 ; Diehl97 . Therefore the fact that recently has been computed quite accurately by using Monte Carlo simulations of spin models constitutes valuable progress. In refs. Hucht ; VaGaMaDi08 ; Ha09 ; Ha10 the three-dimensional XY bulk universality class and a vanishing field at the boundary have been studied, which is relevant for the experiments on He. A quite satisfactory agreement between the experimental results and the theory was found. In refs. DaKr04 ; VaGaMaDi07 ; VaGaMaDi08 ; mybreaking ; PaToDi10 ; mycrossover ; VaMaDi11 ; HuGrSc11 ; mycorrection the Ising bulk universality class and various types of boundary conditions were studied. Note that the mixing-demixing transition of binary mixtures belongs to the Ising bulk universality class.

Experimentally it is easier to access the thermodynamic Casimir force for a sphere-sphere or a sphere-plate than for the plate-plate geometry. Several experiments for these geometries have been performed. In particular, the thermodynamic Casimir force between a colloidal particle immersed in a binary mixture of fluids and the substrate has been directly measured Nature08 ; Gaetal09 . In ref. NeHeBe09 it was demonstrated that the thermodynamic Casimir forces in a colloidal system can be continuously tuned by the choice of the boundary conditions. In addition to homogeneous substrates also chemically patterned ones have been studied So08 ; Tr11 ; Zvyagolskaya . In other experiments, the thermodynamic Casimir force is the driving force for colloidal aggregation, see e.g. refs. Bonnetal09 ; ZvArBe11 .

For simplicity, we restrict the following discussion to the sphere-plate geometry. Furthermore we assume that both the surfaces of the sphere and the plate are homogeneous. In this case, the thermodynamic Casimir force follows the scaling form HaScEiDi98

(2)

where is the radius of the sphere, the distance between the plate and the sphere and . Note that in the literature often , where is the bulk correlation length in the high temperature phase, is used instead of as argument of the universal scaling function. The thermodynamic Casimir force is given by

(3)

where is the reduced free energy of the system.

In the limit , the function can be obtained from by using the Derjaguin approximation (DA) Derja . Instead, in the limit the small sphere expansion BuEi97 can be used. In ref. HaScEiDi98 the function is discussed for the case where both surfaces strongly adsorb the same component of the binary mixture. The behaviour of for three-dimensional systems for is inferred from the exactly known one in two dimensions and from mean-field calculations. In ref. troendle the thermodynamic Casimir force is discussed for colloids close to chemically patterned substrates. To this end the authors employ the Derjaguin approximation. They check the range of applicability of the Derjaguin approximation by using mean-field calculations. In section III of troendle , they discuss for a start also homogeneous substrates.

Here we discuss a method to determine the thermodynamic Casimir force for the sphere-plate geometry by using Monte Carlo simulations of spin models. The generalization to the sphere-sphere geometry is straightforward. The method is based on the geometric cluster algorithm of Heringa and Blöte HeBl98 . In the present work, we simulate the improved Blume-Capel model on the simple cubic lattice. We focus on strongly symmetry breaking boundary conditions at the surfaces of the sphere and the plate. We assume that the surface of the plate prefers positive magnetisation. We study the two cases and , where is the value of the spins at the surface of the sphere. Preliminary results for , indicate that also this case can be efficiently simulated. In our numerical tests we studied spheres of the radii and lattice spacings. For and we considered distances between the sphere and the plate up to . We introduce an effective radius of the sphere which is determined by a particular finite size scaling method. To this end we compute the ratio of partition functions by using a cluster algorithm which is closely related to the boundary flip algorithm MH93A ; MH93B . The index of gives the sign of . In this finite size scaling study we consider , which is large compared with the values used otherwise. Analysing our data we find that replacing by leads to a good scaling behaviour for the whole range of studied here. As result we obtain reliable estimates for the scaling functions and for , where the subscript of gives the sign of . The most striking observation is that in the low temperature phase deviates quite strongly from the Derjaguin approximation already for and likely smaller.

The present work should be seen as pilot study that allows for a number of improvements as discussed in section VI. The method discussed here should apply to a broader range of problems as we also shall discuss in section VI.

The paper is organised as follows. First we define the model. Then we discuss the geometry of the systems that we study. In section IV we discuss the algorithm that is used to compute the thermodynamic Casimir force. In section V we present our numerical results. First we discuss the performance of the algorithm. Then we give our results for the scaling functions . In section VI we conclude and give an outlook. In Appendix A we discuss the finite size scaling method that we used to determine the effective radius .

Ii The Model

We simulate the Blume-Capel model on the simple cubic lattice. It is defined by the reduced Hamiltonian

(4)

where and denotes a site on the simple cubic lattice, where takes integer values and is a pair of nearest neighbours on the lattice. The partition function is given by , where the sum runs over all spin configurations. The parameter controls the density of vacancies . In the limit vacancies are completely suppressed and hence the spin-1/2 Ising model is recovered. Here we consider a vanishing external field throughout. In dimensions the model undergoes a continuous phase transition for at a that depends on . For the model undergoes a first order phase transition, where DeBl04 .

Numerically, using Monte Carlo simulations it has been shown that there is a point on the line of second order phase transitions, where the amplitude of leading corrections to scaling vanishes. Our recent estimate is mycritical . In mycritical we simulated the model at close to on lattices of a linear size up to . From a standard finite size scaling analysis of phenomenological couplings like the Binder cumulant we find . Furthermore the amplitude of leading corrections to scaling is at least by a factor of smaller than for the spin-1/2 Ising model. In order to set the scale we use the second moment correlation length in the high temperature phase. Its amplitude is given by mycorrection

(5)

In the high temperature phase there is little difference between and the exponential correlation length which is defined by the asymptotic decay of the two-point correlation function. Following pisaseries :

(6)

for the thermodynamic limit of the three-dimensional system. Note that in the following always refers to .

Iii Geometry of the System

We have used fixed boundary conditions in -direction and periodic boundary conditions in and directions. In most of our simulations we have fixed for and . In order to check the effect of the second boundary, we have performed some simulations with for and for . In both cases there are layers with fluctuating spins. For we take ,,, and throughout.

To model the spherical object, we have fixed all spins to for

(7)

The distance between the sphere and the plate is given by

(8)

Analysing our data we shall introduce an effective radius of the sphere that we determine by using a finite size scaling analysis in Appendix A below. Taking also into account the extrapolation length at the boundary of the plate we arrive at

(9)

where for the model and boundary conditions discussed here mycrossover . Note that given in ref. mycrossover refers to as location of the boundary. Mostly we have simulated the two cases and . In our study, nearest neighbour pairs that contain a fixed and a fluctuating spin are coupled by the same as nearest neighbour pairs that contain two fluctuating spins. One could imagine to chose the couplings depending on the position on the surface of the sphere, in order to effectively improve the spherical shape of the object.

Iv The algorithm

In order to determine the thermodynamic Casimir force , eq. (3), between the plate and the sphere, we compute finite differences of the reduced free energy of the system

(10)

Note that we chose the difference of the distances equal to one due to the discrete translational invariance of the lattice. Similarly to the proposal of Hucht Hucht , we determine by integrating the corresponding difference of energies

(11)

where and

(12)

The value of is chosen such that , which is the case for . The integration is done numerically, using the trapezoidal rule. It turns out that similar to the study of the plate-plate geometry nodes are needed to compute the thermodynamic Casimir force in the whole range of temperatures that is of interest to us.

Figure 1: Two-dimensional sketch of two spheres at and . The -axis is drawn in vertical direction. The overlapp of the two spheres is given by black squares with white dots. Sites that only belong to the lower sphere are represented by squares with vertical lines and those that only belong to the upper sphere by squares with horizontal lines.

Routinely one would carry out this programme by simulating the systems for and separately using standard Monte Carlo algorithms. In order to avoid systematic errors, should be taken. However the variance of and the CPU-time required for the simulations are increasing with the system size. In fact, first numerical experiments have shown that even for moderate values of and , this approach fails to give sufficiently accurate results.

The expectation value only differs strongly between the two systems for nearest neighbour pairs close to the sphere and this difference decays at least power-like with increasing distance from the sphere. Therefore, computing one would like to concentrate on the part of the system, where the difference in is large and to avoid to pick up variance from the remainder. Routinely one could implement this idea by ad hoc restricting the range of summation in eq. (12). The problem is that this ad hoc restriction leads to a systematic error that needs to be controlled, e.g. by comparing the results of different restrictions. In the following we shall discuss an algorithm that generates the restriction of the summation to the neighbourhood of the sphere automatically and at the same time allows to compute without bias.

This algorithm is closely related to the geometric cluster algorithm HeBl98 . Similar to ref. SwWa86 , two systems are updated jointly. These two systems share the values of and . The distance of the sphere from the boundary is and for system 1 and 2, respectively. In the following the spins are denoted by , where the second index indicates whether the spin belongs to system 1 or 2. Configurations are denoted by . The combined partition function is given by

(13)

where is the reduced Hamiltonian of the system . The cluster update comprises the exchange of spins

(14)

for all sites in the cluster. Performing this exchange of spins for all sites within a cluster is called flipping the cluster in the following. A cluster is build in the usual way: A pair of nearest neighbours is deleted with the probability and frozen otherwise. In our case HeBl98

(15)

Two sites of the lattice belong to the same cluster if they are connected at least by one chain of frozen pairs of nearest neighbours. Various choices of the selection of the clusters to be updated are discussed in the literature. In the case of the Swendsen-Wang SwWa87 algorithm all clusters are constructed and then all spins within a given cluster are flipped with probability . In the single cluster algorithm Wolff one first selects randomly one site of the lattice. Then one constructs only the cluster that contains this site. All spins within this cluster are flipped. Further choices have been discussed in the literature Kerler ; wall .

In our case, the choice of the clusters to be flipped is dictated by the requirement that the fixed spins have to keep their value under the update. The exchange of spins should not be allowed for those sites , where is fixed while is not, or vice versa. In the following these sites are called frozen sites. In Fig. 1 these sites are represented by squares with vertical or horizontal lines on them. In the following a cluster that contains at least one frozen sites is called frozen cluster. In the update all clusters except for the frozen ones are flipped. In order to perform this operation only the frozen clusters have to be constructed.

Now let us discuss how , eq. (11), is computed. Let us assume that after equilibration we have generated a sequence of configurations that are labelled by . We get

(16)

where

(17)

is the standard estimator of . Combining the right hand side of eq. (16) and the particular properties of the update discussed above we arrive at the improved estimator

(18)

where means that at least one of the sites or belongs to a frozen cluster. Note that for our choice of the update

(19)

for all nearest neighbour pairs where neither nor belongs to a frozen cluster.

Whether we have achieved progress this way depends on the average volume taken by frozen clusters. Our numerical study discussed in detail below shows indeed that for the whole range of the total size of all frozen clusters is small compared with the volume of the lattice. Hence the variance of is drastically reduced compared with the standard estimator. Furthermore, since the sum runs only over the frozen clusters, the numerical effort to compute is much smaller than the one to compute the standard estimator.

In the following we shall discuss the details of our implementation. We simulate a large number of distances jointly. We perform the exchange cluster update and the measurement of for all pairs with . In our C-program we store the spins in the array int spins[I_D][L_0][L][L]; where I_D equals the number of distances that are simulated. Routinely one would swap the spins stored in spins[i_c][i_0][i_1][i_2] and spins[i_c+1][i_0][i_1][i_2], where the site represented by [i_0][i_1][i_2] is not a member of a frozen cluster. Instead, in order to save CPU time, we do that for the sites that belong to frozen clusters. This way, the systems with the distances and of the sphere from the plate interchange their position in the array spins. In order to keep track of where the systems are stored in the array spins, we introduce the array int posi[I_D]; where the index i_d equals and posi[i_d] indicates where the system with the distance of the sphere from the plate is stored in spins. Note that parallel tempering simulations are usually organized in a similar way. Instead of swapping configurations, the location of the systems with the temperatures and is interchanged.

In order to get an ergodic algorithm we supplement the exchange cluster algorithm discussed above with standard updates of the individual systems. To this end we perform Swendsen-Wang SwWa87 updates and heat-bath updates. In case of the Swendsen-Wang cluster-update we keep the clusters that contain fixed spins fixed. Those clusters that do not contain fixed spins are flipped with probability .

These updates are performed in a fixed sequence: First we sweep once for each systems over the whole lattice using the heat-bath algorithm. It follows one Swendsen-Wang cluster-update for each system. Then we perform -times the following two steps:

  • For up to perform in a sequence the exchange cluster update and the measurement of for the pairs of distances and .

  • Perform for each system one sweep with the heat-bath algorithm over a sub-lattice characterized by , and .

In the following we shall denote this sequence of updates as one cycle of the algorithm. Since the estimator takes mainly data from the neighbourhood of the sphere, it seems useful to perform additional updates in this region of the lattice. Therefore we introduced the heat-bath updates of the sub-lattice. Here we made no attempt to optimize the shape and size of this sub-lattice nor the frequency of these updates.

V Numerical Results

As a first test of our implementation we simulated systems with the radius and lattices of the size and . We considered the distances between the sphere and the plate. We fixed for and . We simulated both and at the two values and of the inverse temperature. We simulated with the exchange cluster algorithm and without. For these lattice sizes we still get differences of the energies that are clearly larger than the statistical error from simulations that take about one day of CPU-time on a single core of the CPU in the standard way. Our results with and without the exchange cluster algorithm are compatible within the statistical errors. But even for the small lattices sizes used here, the statistical error, using an equal amount of CPU-time, is reduced by more than a factor of 10 by using the exchange cluster algorithm and the improved estimator (18). After these preliminary studies were concluded successfully, we performed an extensive study aiming at physically relevant results. First we summarize the parameters of our simulations. Then we discuss the performance of the algorithm, focussing on the size of the frozen clusters minus the number of fixed spins. Finally we present our results for the thermodynamic Casimir force and compare them with the literature.

v.1 Our simulations

We studied spheres of the radii , and for both and and for only. For these radii, we simulated systems with , where , , , and for , , and , respectively. In the following we shall denote the number of values of by . We simulated these systems for about 100 values of . We adapted the step-size in . Close to , the step-size is the smallest. For all values of we simulated lattices of the size , , and , respectively. In the case of these simulations, we fixed for and . We fixed the parameter of the update cycle to . The sizes of the sub-lattice that is updated in excess times per cycle using the heat-bath algorithm is , , and for , , and , respectively. These parameters are chosen ad hoc. For lack of human time we did not try to optimize them.

In order to check for the effect of the finite values of and , we redid the simulations for in the neighbourhood of on lattices of the sizes and , where we considered both and for . For we simulated in addition lattices of the size with for and with for in the neighbourhood of . For these additional simulations we chose . The sizes of the sub-lattice that is updated in excess are the same as above. In order to keep the figures readable we shall give in sections V.3, V.4 and V.5 below results for at only. The results obtained for at confirm our final conclusions.

Focussing on larger distances we simulated for . For in the range and for in the range we simulated a lattice of the size . The size of the sub-lattice that is updated in excess times per cycle using the heat-bath algorithm is . For less close to we simulated smaller lattices. For these simulations we took . Finally we performed simulations for at the distances using a lattice of the size in the range also taking .

For each of these simulations we performed up to update cycles. For and , , and we simulated at only. In all our simulations we used the Mersenne twister algorithm twister as pseudo-random number generator. Our code is written in C and we used the Intel compiler. Simulating at , our Swenden-Wang update takes about s and the local heat-bath update about s per site on a single core of a Quad-Core AMD Opteron(tm) 2378 CPU. The exchange cluster, including the measurement (18) takes roughly s per site of the frozen clusters. In total our simulations took roughly the equivalent of 50 years of CPU time on one core of a Quad-Core AMD Opteron(tm) 2378 CPU.

v.2 Size of the frozen clusters

The CPU-time needed for one update cycle depends on the size of the frozen clusters. Also the variance of the estimator (18) compared with the standard one can be small only if the size of the frozen clusters is small compared with the volume of the lattice. Therefore we shall discuss the size of the frozen clusters and its dependence on the parameters , , and the lattice size in some detail below.

To give a first impression, we plot

Figure 2: We plot the average size of the frozen clusters minus the number of frozen sites for the radius and the distance as a function of the inverse temperature . The upper and the lower curve represent the results for and , respectively.

in Fig. 2 the average size of the frozen clusters minus the number of frozen sites for the radius and the distance , i.e. for the pair of systems with and , for and . For decreasing in the high temperature phase, the average size of the frozen clusters for and seems to become identical. As the critical point is approached, the size for becomes clearly larger than for . This remains true in the low temperature phase. In both cases, the maximum is located close to . We find and for and , respectively. Note that this is only a small fraction of the total volume of the lattice, .

Let us now discuss the behaviour of in more detail. Since for , the frozen clusters consist of the frozen sites only, i.e. . At low temperatures, the configurations become ordered. Therefore also in this limit the frozen clusters consist of the frozen sites only.

From our numerical data we see that for fixed and lattice size, is monotonically increasing with increasing distance . For , seems to converge to a finite value . In the high temperature phase, , this value is the same for and . We fitted our data in the high temperature phase with the Ansatz . For we get quite reasonable fits with the result . Fits for larger radii, where we have less data apparently give smaller values of . However, since should increase with increasing radius , this observation seems to be an artifact of our finite data set. A theoretical understanding of the behaviour of , relating with known critical exponents, analogous to ref. HeBl98 , would be desirable.

Next let us study the dependence of on the linear lattice sizes and in the

Figure 3: We plot as a function of for , and focussing on the neighbourhood of the critical point. We give results for the lattices sizes , and .

neighbourhood of the critical point. In Fig. 3 we plot as a function of for , and for the lattice sizes , and . For and the estimates of for the different lattice sizes are consistent within the statistical error. Instead for a dependence of on the lattice size can be seen. For most values of in this range is increasing with increasing lattice size. However for the value of for seems to be a bit larger than those for the lattices sizes and . For we get a qualitatively similar behaviour. The same holds for for both and .

Finally let us discuss the behaviour of at the critical point in more detail. In Fig. 4 we plot as a function of

Figure 4: We plot as a function of for , and . We give results for the lattice sizes , and .

for , and . The value of seems to increase somewhat faster than linearly with the distance . Differences between the results for the lattice sizes , and are clearly visible. Similar observations can be made for and also for other radii.

In table 1 we give at for the radii , and for

3.5 1096(3) 1495(5) 269(1)
4.5 2464(6) 04201(19) 565(2)
7.5 3862(8) 11203(67) 970(2)
Table 1: for , and at the critical point for . The results are taken from the lattice sizes , , for , and , respectively.

the distance , i.e. , , and , respectively. As one might expect, is increasing with increasing radius . The ratio of for and is roughly the same for all three radii. In contrast, we find that is increasing more rapidly with increasing for than for . For the results for and suggest that is increasing less than linearly with the radius , while the opposite seems to be true for .

In section IV we argued that the estimator (18) reduces the variance if . For we have verified that this is indeed the case for the full range of and all radii and lattice sizes studied here. Our preliminary results for obtained at indicate that in this case the exchange cluster algorithm is less efficient than for . Still it seems plausible that also for meaningful results for the thermodynamic Casimir force can be obtained.

v.3 The thermodynamic Casimir force: Effects of the finite lattice

Figure 5: We plot as a function of for and . The lattice sizes are , , and . In the upper and the lower part of the figure we give results for and , respectively.

For corrections to the limit decay exponentially with increasing lattice size. Therefore corrections should be negligible for . Instead, at the critical point, we expect that corrections decay power-like. For we checked explicitly how depends on the lattice size. To this end we simulated in the neighbourhood of in addition to the lattice of the size ones of the sizes and . In Fig. 5 we plot for our results for the difference of energies as a function of . We give our results for and in the upper and the low part of Fig. 5, respectively. We see differences between the results obtained for the and the larger lattices that are statistically significant for and for and , respectively. Note that , using eq. (5).

In Fig. 6 we give the numerical results of

Figure 6: We plot as a function of for , , and . In the upper part of the figure we give the result obtained by using a lattice for the full range of that we consider. In the lower part of the figure we focus on the neighbourhood of , plotting the results for the lattice sizes , , and .

for , and . In the upper part of the figure we give the result of the numerical integration of obtained by using the data for the lattice. In the lower part of the figure we focus on the neighbourhood of , where effects due to finite and become important. As check, we computed in an alternative way: We replaced the estimates of from the lattice, where available, by those from the or the lattice. We find that the behaviour of is changed only little. The difference between these results is of a similar size as the statistical error. Note that the statistical error of is the accumulated error of over the range of the integration.

We conclude that in a neighbourhood of the critical point the effect of the finite lattice size on is clearly visible. However, since the range of where this is the case, is quite small, is affected less. The amplitude of this effect can be estimated by comparing the results obtained for different lattice sizes.

v.4 The scaling function for

In this section we discuss our results for the scaling function for . In a first step we check how well our results obtained for different radii collapse on a single curve. To this end we make use of the effective radius determined in Appendix A, see table 5, and the effective distance (9). In particular in the following analysis of our data. In the plots given below we show the statistical error of the data only. The uncertainty of and still comes on top of this. Essentially it is a two percent uncertainty of the scales.

Figure 7: We plot as a function of for and the radii , and . In the upper part we give results for , which means , and for , and , respectively. In the lower part we give results for , which means , and for , and , respectively. Note that is used.

In figure 7 we plot as a function of for and . We multiplied by to reduce the dependence on . Note that for small values of , and large ones HaScEiDi98 ; BuEi97 , where and mycritical . We took the data obtained by using , and lattices for the radii , and , respectively. For this choice is roughly the same for the three radii considered.

In particular for we see an excellent collapse of the data obtained for the three different radii. Note that things look far worse when and are used instead of and . For the collapse of the data is slightly worse than for . In particular for we see some deviation of the result for from that for the two other radii. This is not too surprising, since means for and also the values of that correspond to are far from .

We conclude that the data scale well and our numerical results should describe the scaling limit. At the level of our accuracy, one should interpret data obtained for with caution.

In Fig. 8 we plot as a function of for various values of in the range . For each value of , we took the largest radius, where data are available.

Figure 8: We plot as a function of for various values of . Note that the index of indicates . We compare the Derjaguin approximation in the limit with our Monte Carlo data for various finite values of . In particular we have used our data for and giving , and , , , and giving , , , and , respectively, and and , and giving , , and , respectively.

For comparison we plotted the Derjaguin approximation in the limit . For a brief discussion of the Derjaguin approximation see e.g. ref. HaScEiDi98 . As input we used the results for and obtained in ref. mycorrection .

We find that the shape of changes only little with varying . In particular the location of the minimum changes very little with varying . The absolute value of is increasing with increasing , which is consistent with the fact that for large , BuEi97 . These observations are fully consistent with those of ref. troendle . In Fig. 2 a of troendle is plotted as a function of . Note that the index corresponds to our and to our . For , the mean-field result for is fully consistent with the Derjaguin approximation for four dimensions, while for a small deviation can be seen.

In the following we study in more detail how the minimum of as function of depends on . Fortunately, is large enough, such that finite and effects are well under control. We located the minimum of by finding the zero of . To this end, we interpolated linearly in . Using these results, we computed . In Fig. 9 we plot as a function of using our data obtained for the radii , , and . We see that for the results obtained for different radii agree quite nicely. At the estimate seems plausible. With further increasing the value of is slightly decreasing. For we get . For clear deviations between the results obtained by using different radii can be observed. Therefore we abstain to give a conclusion for the scaling limit in this range. Note that the Derjaguin approximation gives for the limit . The Derjaguin approximation predicts that is increasing with increasing distance . If that prediction is correct, should show a maximum in the interval .

Figure 9: We plot our estimates of as a function of . Monte Carlo results obtained by using , , and are given. These are compared with the Derjaguin approximation in the limit . Note that .

Our results can be compared with the upper part of FIG. 2 of ref. HaScEiDi98 . Note that the authors of ref. HaScEiDi98 consider as variable. They obtain for using the Derjaguin approximation. This discrepancy should be due to the different numerical estimates of the scaling function of the plate-plate geometry. E.g. following FIG 3 of ref. HaScEiDi98 , while our present estimate is mycorrection . Then, according to ref. HaScEiDi98 increases and reaches a maximum at assuming the value . Then, with increasing it is decreasing and reaches in the limit . In particular for the value assumed at we see a clear discrepancy with our results. It is beyond the scope of the present work to discuss the reliability of the first order -expansion that has been used in ref. HaScEiDi98 to get the small-sphere expansion.

Figure 10: We plot our estimates of as a function of . Monte Carlo results obtained by using , , and are given. The Derjaguin approximation in the limit is given by a black circle. The Derjaguin approximation for finite values of is represented by a dashed black line. Note that .

In Fig. 10 we plot as a function of . We find a nice collapse of the data obtained from different radii of the sphere. In contrast to , this is also the case for . Taking into account our Monte Carlo data only, seems to be a monotonically decreasing function of . However we note that for our numerical estimates of are overshooting obtained from the Derjaguin approximation. Also note that the Derjaguin approximation predicts that increases with increasing , which is indicated by a dashed black line in Fig. 10. Therefore should show a maximum in the interval . In the lower part of Fig. 2. of ref. HaScEiDi98 is shown as a function of . For it shows a maximum (minimum with the sign of ref. HaScEiDi98 ) taking a value slightly larger than . Taking our Monte Carlo data, is still slightly increasing at and takes the value .

Finally we discuss the scaling function at the critical point, i.e. . In figure 11 we plot our numerical results for as a function of . At the critical point we have to consider finite and effects. Therefore, for we plotted our results obtained by using lattices of the sizes , and . In particular for large we see a clear deviation of the results obtained for the lattice from the other two. On the other hand the estimates obtained for the and lattices are consistent within the statistical error. Therefore deviations from the limit should be of similar size as the statistical error. In the case of we plot the data for the largest lattice size available. We find a quite good collapse of the data obtained for different radii. For we get . We see that is almost linearly decreasing with increasing . For we read off . The last three data points for seem to suggest that is increasing again for . However this is not statistically significant and is likely an artifact.

Using the Derjaguin approximation, we get and as for , is increasing with increasing . Therefore it seems plausible that has a maximum in the range .

Figure 11: We plot as a function of . Monte Carlo data obtained by using , , and are given. In the case of we compare results obtained by using the lattice sizes , , and . Note that .

v.5 The scaling function for

Next let us turn to the scaling function for . In figure 12 we plot as a function of for and . The limit is obtained by using the Derjaguin approximation. Unfortunately we do not have data for the scaling function for the plate-plate geometry for and the amplitude of in this range is still significant. Therefore we have extrapolated in two different ways: (1) for and (2) is linear in the interval and , , for we take . The true should be enclosed by these two extrapolations. In figure 12, the results based on these two extrapolations are denoted by Derjaguin 1 and Derjaguin 2, respectively. For comparison we plot our Monte Carlo result obtained for and which corresponds to . We observe that the two curves are very similar in the high temperature phase. Since has its maximum in the low temperature phase, the same holds for computed by using the Derjaguin approximation. In contrast, for we find that the maximum is located in the high temperature phase, even though very close to the critical point. Furthermore, at is much smaller than in the low temperature phase. Also the decay of as seems to be much faster for than for .

Figure 12: We plot as a function of for obtained by using the Derjaguin approximation and using our Monte Carlo data for and . Note that the index of indicates . In the case of the Derjaguin approximation we have used two different extrapolations of the scaling function for as discussed in the text.
Figure 13: We plot as a function of for various values of . Note that the index of indicates . We use Monte Carlo data for