- Email: [email protected]
Critical Casimir forces for a particle between two planar walls
Journal of Non-Crystalline Solids 407 (2015) 376–383
Contents lists available at ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Critical Casimir forces for a particle between two planar walls Oleg Vasilyev a,b, Anna Maciołek a,b,c,⁎ a b c
Max-Planck-Institut für Intelligente Systeme, Heisenbergstr. 3, D-70569 Stuttgart, Germany IV Institut für Theoretische Physik, Universität Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, PL-01-224 Warsaw, Poland
a r t i c l e
i n f o
Article history: Received 27 June 2014 Received in revised form 1 September 2014 Accepted 2 September 2014 Available online 19 September 2014 Keywords: Colloids; Phase transitions and critical phenomena; Classical spin models; Monte Carlo methods
a b s t r a c t We consider a cube placed between two planar surfaces bounding a three-dimensional Ising slab. This system can be viewed as the magnetic analog of a colloidal particle immersed in a fluid, which is confined by two parallel walls. Near the bulk critical point of the Ising ferromagnet the cube gets exposed to the potential of the fluctuation that induces the critical Casimir force. Using Monte Carlo simulations we study the dependence of this potential on the relative boundary conditions at the two surfaces and at the cube within a wide range of temperatures above and below the bulk critical temperature. In order to calculate the critical Casimir force and its potential, we adopt an approach based on an integration scheme of free energy differences. The scaling functions of the corresponding critical Casimir forces are also determined. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Two objects immersed in a fluid which is brought to its bulk critical point Tc experience the so-called critical Casimir interaction [1]. This effective interaction arises due to modifications of the structure of the fluid and due to the restrictions of the fluctuation spectrum of its order parameter by the surfaces of the objects. Accordingly, not only the temperature T but also the adsorption properties of both objects and their geometry determine the sign, the strength and the range of the critical Casimir forces (CCFs). In asymptotic regimes (T → Tc and all distances large compared with microscopic lengths) finite-size scaling [2] holds and the CCFs are described by universal scaling functions. Universality permits to use the simplest possible representative of the same universality class in order to calculate these scaling functions. In this sense the Ising model represents systems such as simple one-component fluids or binary liquid mixtures. The simple geometry of two parallel planes can be realized experimentally by growing a wetting film; exactly the wetting films were used to provide first reliable evidence for the existence of the CCFs [3–6]. On the other hand, even for this simple geometry determination of the scaling functions from corresponding model systems is a challenging task. Theoretical approaches which incorporate critical fluctuations beyond the Gaussian approximation and can take into account a dimensional crossover occurring in films are intractable for most of the relevant experimentally boundary conditions (BCs). Monte Carlo (MC) simulations offer a very useful alternative approach, ⁎ Corresponding author at: Max-Planck-Institut für Intelligente Systeme, Heisenbergstr. 3, D-70569 Stuttgart, Germany. E-mail address: [email protected] (A. Maciołek).
http://dx.doi.org/10.1016/j.jnoncrysol.2014.09.005 0022-3093/© 2014 Elsevier B.V. All rights reserved.
which is however not free from difficulties. Within this approach the range of sizes of the system is strongly limited by the steeply increasing computational costs. For the thicknesses and the aspect ratios of slabs accessible in simulations, numerical data do not collapse and corrections to scaling are necessary. Renormalization-group analyses reveal that there is a whole variety of sources for corrections to scaling which arise from bulk, surface, and finite-size effects [2]. In addition to the leading bulk corrections to scaling (with the correction-to-scaling exponent ω ≃ − 0.8) one expects corrections due to the finite aspect ratio, corrections due to the boundary conditions (which are especially large for the symmetry breaking boundary conditions) and, moreover, nextto-leading corrections might occur for narrow films. In spite of these difficulties the Casimir scaling functions in three spatial dimensions (3d) have been obtained rather accurately for experimentally relevant universality classes with a variety of BCs and surface fields [7–14]. In the present paper we study numerically the CCFs acting on a cube located between two planar walls as a function of temperature and separation for different BCs. The cube mimics a colloidal particle positioned near the wall. We use MC simulations of the Ising model on a cubic lattice as a representative of the Ising universality class of critical phenomena. The aim of this work is twofold. We would like to check whether the approach that we used in our previous studies for a slab geometry [8–11] can be successfully adopted for the present geometry. This approach is based on the so-called coupling parameter method where the difference of free energies is expressed as an integral over the mean energy difference (see, e.g., Refs. [8] and [15]). Moreover, we address the question of the mechanical stability of a cube particle for the cases when (a) both walls and the particle have the same adsorption preference (+| + |+), (b) the walls have opposite adsorption preferences (+| + |−), and (c) the cube particle is
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 1. The geometry under consideration (the cross-section in the xz plane): (a) the spherical particle of the diameter a at the distance D from the wall; (b) an appropriate geometry on the lattice for a cubic particle with the edge length a at the separation D from the bottom wall. The total system size is Lx × Ly × (Lz + 2), and the separation of the particle from the top wall is Lz − D − a, where Lz is the number of layers of free spins.
377
suppressed. The focus of that study was on symmetry breaking boundary conditions at the surfaces of the sphere and the plate. The use of the geometric cluster algorithm and the introduction of an effective radius of the sphere and the effective distance from the wall have allowed the author to obtain a good data collapse. Still another MC simulation method to compute the critical Casimir force acting between the disk and the wall in two-dimensional Ising systems has been proposed in Ref. [21]. This method is analogous to the one used in experiments reported in Refs. [19,20] for a colloidal particle where the Casimir potential is determined from the distribution of the position of the particle above the wall. The interaction between two particles in the presence of the bulk ordering field recently has been studied by using a local field integration method [22]. Our presentation is organized as follows. In Section 2 we briefly present the relevant theoretical background, introduce our model, and describe the numerical method employed in order to compute the CCFs and its scaling functions from the MC simulation data. Section 3 contains our results. We provide a summary and conclusions in Section 4. 2. The model and the method
neutral and the walls have opposite adsorption preferences (+|0|−). These situations can be mimicked in our Ising system by choosing an appropriate BC on the lattice sites forming boundaries of the three objects. Because of unavoidable limitations of the size of the system, we cannot expect to reach the asymptotic region of the universal behavior in our simulations. Nevertheless, after rescaling procedure for some BCs we have obtained a reasonable data collapse. Recently, a different MC simulation method for spin models which mimic the sphere–planar plane geometry has been developed [16]. The geometry of a sphere or a cube near a single planar surface is relevant for colloidal systems or in experimental setups of high-precision force measuring devices and methods, such as atomic force microscope (AFM) [17], surface force apparatus (SFA) [18] or total internal reflection microscopy (TIRM) [19,20]. For this geometry theoretical predictions for the CCFs are much more limited than for films. The MC simulation method employed in Ref. [16] is similar to the one used in [7] for the plate–plate geometry, where differences of the free energy are computed by integrating differences of the energy over the inverse temperature. Simulations have been performed for the improved Blume–Capel model on the simple cubic lattice, which shares the universality class of the three-dimensional Ising model but has the advantage over it that the amplitudes of the leading correction to scaling are substantially
(a)
2.1. Theoretical background Let us consider a spherical colloidal particle of a diameter a immersed in critical binary mixture near a wall, see Fig. 1(a). As already mentioned in the Introduction, at the critical point of the demixing transition of the binary solvent, the CCFs, which are the analogs of the electromagnetic Casimir force [23], emerge between the wall and the colloidal particle as a result of critical fluctuations of the order parameter of the solvent and the critical adsorption of a preferable component of the mixture on the surfaces of wall and the colloid. For a binary mixture the order parameter is the deviation of the concentration of the one component of the mixture from its value at the critical point. The CCFs between the wall and the particle are obtained from the free energy UC(β, D, a) of a system consisting of the colloidal particle of the characteristic linear size a immersed in the fluid at the inverse temperature β = 1/kBT at the separation D from the planar wall: f C ðβ; D; aÞ ¼ −
∂U C ðβ; D; aÞ : ∂D
ð1Þ
Note that the above definition of the CCFs differs from the one used to define the CCFs between two planar parallel plates immersed in the
(b)
(c)
Fig. 2. Bond arrangement for the computation of the free energy difference in Eq. (9) (see main text). The crossover Hamiltonian Hcr (a) describes the system which interpolates between those described by the Hamiltonians H0 (b) and H1 (c).
378
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 3. Results for (+[+]+) boundary conditions as a function of the distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The cube has the edge length of a = 4 (in units of lattice spacing) whereas the slab has the size of 42 × 42 × 45.
fluid. If D is the separation between the two plates then for large areas A of plates the total free energy F(β, D, A) of the system can be written as F ðβ; D; AÞ ≡ ADf ðβ; DÞ h i bulk −1 ex ¼ A Df ðβÞ þ β f ðβ; DÞ ;
ð2Þ
where f bulk(β) is the bulk free energy density at a given temperature. The excess free energy f ex per area contains two D-independent surface contributions in addition to the finite-size contribution f ex(β, D) − f ex(β, ∞) which vanishes for D → ∞. The D-dependence of the latter gives rise to the critical Casimir force fC per unit area A: ex
f C ðβ; DÞ=kB T ≡ −∂f ðβ; DÞ=∂D:
below (−) the bulk critical temperature Tc with ξ− 0 = 0.243(1) and ξ+ 0 = 0.501(2) [25] for the Ising model on the simple cubic lattice. The value of ξ+ 0 quoted here refers to the amplitude of the second moment correlation length ξ2nd ; however, ξ=ξ2nd ≃1 for β b βc for the Ising model [24]. The scaling function K also depends on the BCs imposed by the particle and the wall on the order parameter, which reflects the adsorption preferences of the surfaces of these objects for one of the two species of the mixture. 2.2. The model We consider the Ising model defined on a three-dimensional simple cubic lattice with the lattice spacing ℓ = 1 via the Hamiltonian
ð3Þ H ¼ −J
From the general theory of finite-size scaling [2] it follows that the CCFs between the wall and the particle can be cast in the following scaling form f C ðβ; D; aÞ≈
kB T D D K ; a ξ a
ð4Þ
where the universal scaling function K depends on the ratio D/ξ and also on the ratio D/a. In the critical region and vanishing ordering bulk field −ν h = 0, the bulk correlation length is ξ = ξ± , where the reduced 0 |τ| temperature is defined as τ = (T − Tc)/Tc = (βc − β)/β. The critical exponent takes the value of ν ≈ 0.63 for the Ising bulk universality class in 3d [12,24], ξ± 0 are nonuniversal amplitudes above (+) and
(a)
X
si s j ;
ð5Þ
hi; ji
where si = ± 1 denotes the spin variable. J, which we set equal to 1, is the spin–spin coupling constant and the sum 〈i, j〉 is taken over all nearest neighbor pairs of sites i and j on the lattice. Temperatures and energies are measured in units of J. The inverse critical temperature is βc = 0.2216544(3) [25]. We consider a slab geometry, i.e., a lattice of fluctuating spins of the size Lx × Ly × Lz with periodic BCs in the x and y directions (in which the system has linear extensions Lx and Ly). There are two additional top and bottom layers as shown in Fig. 1(b). All spins in the top (st, z = Lz + 1) layer are fixed at the value + 1 whereas all spins in the bottom (sb, z = 0) layer are fixed at the value + 1 or −1 for different types of BCs. We have modified the coupling constants
(b)
Fig. 4. Results for (−[+]+) boundary conditions as a function of the distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The system size and the edge length of a cube are the same as in Fig. 3.
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
379
(b)
Fig. 5. Results for (−[0]+) boundary conditions as a function of distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The system size and the edge length of a cube is the same as in Fig. 3.
between fixed spins at the top and bottom boundaries and those in the adjacent layers. In our present simulations we consider the infinite strength of these modified bonds, which correspond to the infinite surface fields acting on the next-to-surface layers of the system. We also fix a3 spins sp of the cubic particle of linear size a. The spins sp can take values sp = ± 1(for [+] and [−]) and 0 (for [0]). The 0 value corresponds to the neutral particle that does not interact with the components of the binary mixture. As for the walls, we have modified the bonds between the spins at the surface of a cube and their nearest neighbors and assume that their strength is infinite. We denote the type of BCs as (sb[sp]st) (as it was done in Ref. [26] for computation of torques in 2d systems), for example, (−[0]+) BCs correspond to negative spins sb = −1 at the bottom layer, zero spins sp = 0 for the particle and positive spins st = +1 for the top layer. In the present paper we consider three cases: (i) a particle between two attractive walls (+[+]+), (ii) a particle between the repulsive and attractive walls (−[+]+), and (iii) a neutral particle between walls of opposite signs (−[0]+).
2.3. Computation of free energy differences In general, the Monte Carlo methods are efficient for the computation of quantities which can be expressed as ensemble averages, such as energy, magnetization, specific heat or magnetic susceptibility of a system. The free energy differences can be nevertheless computed
(a)
accurately by using the so-called “coupling parameter approach” [15]. Recently, this method was used for the computation of the Casimir force scaling function for the slab geometry [8–11]. For the lattice with spacing ℓ = 1 we rewrite Eq. (1) as the finite difference: f C ðβ; D−1=2; aÞ ¼ −ΔU C ðβ; D; aÞ ¼ −½U C ðβ; D; aÞ−U C ðβ; D−1; aÞ ¼ −½ F 0 −F 1 :
ð6Þ The free energy UC(β, D, a) = F0 corresponds to the system where the cubic particle of the linear size a is at the separation D from the bottom wall (the wall is assumed to be located between the layer at z = 0 and the layer at z = 1) whereas the free energy UC(β, D − 1, a) = F1 corresponds to the similar system but with the separation D − 1. See Fig. 1 where distances are defined. In order to compute the free energy difference in Eq. (6) we introduce the crossover Hamiltonian H cr ðλÞ ¼ ð1−λÞH 0 þ λH 1 ;
ð7Þ
where λ ∈ [0, 1] is the crossover parameter. This Hamiltonian describes the system consisting of Lz + 1 layers of free spins in which the coupling constant associated with certain bonds is equal to λ and with some other bonds is equal to 1 − λ, as depicted in Fig. 2(a). These bonds of different strength are arranged around two selected horizontal layers of spins. The bottom selected layer is located at zbl = D − 2 (at distance 2 below the
(b)
Fig. 6. Results for (+[+]+) boundary conditions for different distances D: (a) the Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
380
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 7. Results for (−[+]+) boundary conditions for different distances D: (a) the critical Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
particle) and the top selected layer is located at ztl = D + a + 3 (at distance 2 above the particle), as shown in Fig. 2(a). In this figure, fixed spins of walls and of the particle are denoted by filled circles, while free (fluctuating) spins are denoted by open circles. For λ = 0, the bottom selected layer is incorporated into the bulk of the system and the top selected layer is excluded from interactions as shown in Fig. 2(b). In this configuration, which is described by the Hamiltonian Hcr(λ = 0) = H0 corresponding to the system with free energy F0 plus the free energy of isolated 2D layer, the bottom face of the cubic particle is at the distance D (D layers of free spins) from the bottom wall, and the upper face is at the distance Lz − D − a from the top wall. For λ = 1 the bottom selected layer is excluded from interactions and the top selected layer is incorporated into the bulk as shown in Fig. 2(c). The total system height is Lz (as for the case λ = 0). In this configuration, which is described by the Hamiltonian Hcr(λ = 1) = H1 corresponding to the system with free energy F1 plus the free energy of isolated 2D layer, the bottom face of the cubic particle is at the distance D − 1 from the bottom wall and the top face is at the distance Lz − D − a + 1 from the top wall. Locations of the selected layers are rather arbitrary; the removed layer should be somewhere in between the particle and the bottom wall and the inserted layer should be somewhere in between the particle and the top wall. As a function of the coupling parameter λ, Hcr(λ) interpolates between H0 and H1 as λ increases from 0 to 1, but the configurational space (amount and arrangement of free spins) remains the same. Only the interaction strength of some bond is changed. The free energy F cr ðλÞ ¼ −β1 ln ∑C expð−βHcr ðλÞÞ, where the sum is taken over all configurations fCg of free spins and interpolates between F0 and F1. The free
(a)
energy difference F1 − F0 can be expressed as F1 − F0 = ∫ 10Fcr′(λ)dλ where Fcr′ is the derivative of Fcr(λ) with respect to the coupling parameter λ: dF cr ðλÞ ¼ dλ
X C
ðH 1 −H 0 Þe−βHcr X e−βHcr ðλÞ C
ðλÞ
¼ hΔH ðλÞifCg ;
ð8Þ
which takes the form of the canonical ensemble average h…ifCg of ΔH(λ) ≡ H1 − H0 over all configurations fCg for the fixed value of the coupling parameter λ and therefore it can be computed via MC simulations. Actually, ΔH(λ) equals to the thermal average of the interaction energy of pair of spins connected by appearing bonds minus the interaction energy of spins connected by vanishing bonds in Fig. 2(a) for a given value of λ. As a result one can express the difference in free energies as an integral over canonical averages (see, e.g., Ref. [15]): Z ΔF ¼ F 1 − F 0 ¼
1 0
hΔH ðλÞifCg dλ:
ð9Þ
In that formula contributions of 2D layers for F0 and F1 are mutually canceled. Finally we obtain the expression for the Casimir force for separation D − 1/2 as Z f C ðβ; D−1=2; aÞ ¼ −
1 0
hΔH ðλÞifCg dλ:
ð10Þ
Here, the force as the finite difference between separations D and D − 1 is assigned to separation D − 1/2. Let us note that due to finite size
(b)
Fig. 8. Results for (−[0]+) boundary conditions for different distances D: (a) the critical Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
381
(b)
0.518 Fig. 9. The critical Casimir force at the critical point Tc multiplied by: (a) D21/2/(πa) and (b) D1.518 for different a and D for (+[+]+) boundary conditions. 1/2 /a
of our system the particle interacts also with the upper wall, which is included in Eq. (10). From the knowledge of the critical Casimir force f C ðβ; D−1=2; aÞ ¼ −∂U ¼ U C ðβ; D−1; aÞ−U C ðβ; D; aÞ , we can compute the potential of ∂D the CCFs UC as a function of a distance from the wall D with respect to some given separation D0. We select as a reference plane from which the potential is measured the midplane of the system D0 = Lz/2. It means that the potential of the CCFs is zero at this plane. For a lattice model (a discrete difference instead of derivative), the potential is computed in accordance to the formula C
U C ðβ; D; aÞ ¼
8 D X > > > > f C ðβ; z−1=2; aÞ; − > < > > > > > :
z¼D0 þ1 D−1 X
f C ðβ; z−1=2; aÞ;
D≥D0 :
ð11Þ
DbD0
z¼D0
3. Numerical results 3.1. Force and its potential as functions of the distance from the wall We compute the critical Casimir force fC (β, D − 1/2, a) and its potential UC(β, D, a) for a cubic particle of edge length a = 4 placed in the system of 42 × 42 × 45 free spins. Results as a function of the separation D − 1/2 for different values of the inverse temperature β are plotted in Fig. 3(a) for (+[+]+), in Fig. 4(a) for (−[+]+), and in Fig. 5(a) for (−[0]+). A potential UC(β, D, a) of the CCFs for these cases is presented in Fig. 3(b), Fig. 4(b), and Fig. 5(b), respectively. In all presented plots we have adopted the notation fC(β, D1/2, a) for fC(β, D − 1/2, a) and D1/2 for D − 1/2. (Note, that because of our choice of the selected layers for the computation of the free energy difference, the distance at which we can calculate the CCFs varies between 4.5 and 35.5). As can be expected, for (+[+]+) boundary conditions the cube is attracted to both walls equally strong so that the force is vanishing in the middle of the slab and its potential has a maximum there. Because the force is measured as a function of a distance from the bottom wall, it becomes positive (repulsive) for D1/2 N D0 = Lz/2 − 1/2 as a result of dominating attraction to the top wall. Interestingly, for temperatures sufficiently far below or above the critical temperature (the blue and the red curve in Fig. 3(a)), the maximum of the potential is very flat. This indicates that the particle can be positioned in the middle part of a slab and stay there for a relatively long time. The potential wells at the walls are the deepest at some temperature slightly above Tc. These features are consistent with the ones of the CCFs between two planar (+,+) surfaces [8], which are attractive, possess a minimum above Tc,
and decay exponentially fast with the deviation from the critical temperature τ. For (−[+]+) boundary conditions the symmetry is broken; the cube is attracted to the (+) wall and repelled from the (−) wall. Below Tc, the repulsion from the (−) wall is much stronger then the attraction to the (+) wall. This is due to the presence of the interface between predominantly positive and negative magnetized phases located between the (−) wall and the cube. The resulting potential decays rapidly and monotonically from the large positive values at the (−) wall towards its minimum located at the (+) wall. The depth of the minimum is comparable with that for the (+[+]+) case. Above Tc, the interface vanishes and the repulsion from the (−) wall, as we know it from the case of two planar walls with (−,+) boundary conditions, becomes very weak. In contrast, the attraction from the (+) walls attains its extreme value above Tc. The net effect for the cube is that the force acting on it is negligibly small around the midpoint of a slab and its potential is very flat. The CCFs and its potential are asymmetric. In the case of the “neutral” cube confined between the opposite walls, (−[0]+) boundary conditions have been numerically the most challenging one — notice the big error bars in Fig. 5. This is because the order parameter at the cube is very small and the resulting CCFs are relatively weak as compared to the previous cases. Recall that for the planar walls with (0,+) (or equivalently (0,−) boundary conditions [10]), the CCFs are repulsive with a maximum below Tc — just as for the (−,+) case, however, the amplitude at its maximum is about 10 times smaller than for the (+,−) boundary conditions. The neutral cube is repelled equally strong from both walls. Similar to the (−[+]+) case, below Tc the boundary conditions at the walls create an
Fig. 10. Results for (+[+]+) boundary conditions for different distances D: the rescaled Casimir force D1.518 1/2 fC as a function of D1/2/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
382
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 11. Results for a sphere of radius R = 5.6 placed near the wall with the same boundary conditions (surface spins fixed at the value of + 1) for different distances D: (a) the Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(2πR))fC as a function of D1/2/ξ. The system size is 50 × 50 × 105.
interface between oppositely magnetized phases which is distorted by a presence of the cube. However, for the neutral cube this distortion is symmetric. As a result, the critical Casimir force varies approximately linearly with the distance across the slab, except for the immediate vicinity of the walls, i.e., one or two lattice sites. The slope of this linear variation decreases as the temperature approaches the critical temperature. At T = Tc, the linear behavior is narrowed to the middle part of the slab. Accordingly, the critical Casimir potential is quadratic and strongly localizes the particle in the middle of the slab. Above Tc, the potential is very flat so that the cube can take stable positions within relatively large region near the center of the system.
e ðxÞ ¼ 2π∫∞ d‘ 1=‘3 θð‘xÞ, i.e., it does not where the scaling function K 1 depend on a. Here, θ(x) is the Casimir force scaling function for the planar plane (film) geometry. At the critical point β = βc one has f C ðβc ; D; aÞ πa ≈ 2 θð0Þ: kB T c D
In the opposite limit of D ≫ a/2, the form of the scaling for the CCFs follows from the small sphere expansion [28] f C ðβ; D; aÞ≈
3.2. The force as a function of temperature In order to study the temperature dependence of the CCFs acting on a cube, we compute fC(β, D1/2, a) for various values of inverse temperature β at several separations D1/2 = 4.5, 5.5, …, 19.5. Results are plotted in Fig. 6(a) for (+[+]+), Fig. 7(a) for (−[+]+) and Fig. 8(a) for (−[0]+). For the (+[+]+) case, the features of the CCFs are similar to the features of CCFs for the planar plane geometry with fixed (+,+) boundary conditions, i.e., the force is attractive and has its minimum above Tc. The same holds for (−[+]+) and (−[0]+) systems if compared to the planar plane geometry with (−,+) and (−, 0) boundary conditions, respectively. In these cases the force is attractive with a maximum above T c. Data for (−[0]+) systems suffer from large computation errors in the region of low temperatures. 3.3. Scaling function of the critical Casimir force Because the size of the systems for which we were able to perform our simulations is rather limited, we cannot expect to determine the Casimir force scaling function K Dξ ; Da . The additional difficulty is the presence of the second scaling variable y = D/a. We could plot the data as a function of the first variable while keeping the second one fixed. This, however, would require to consider a cube with different size a for each distance D, which is computationally very costly. What we try to do in the present paper is to investigate the leading dependence of the CCFs on D and on a. As a guide we use results for a spherical particle of radius a/2 located at a distance D of closest approach from the flat surface of a substrate. The CCFs acting on the particle in the limit of D ≪ a/2 is given by the Derjaguin approximation as [20,27–29]: f C ðβ; D; aÞ a=2 e D ≈ 2 K kB T ξ D
ð12Þ
ð13Þ
kB T D −xϕ −1 −2x −1 P ðD=aÞ þ O ðD=aÞ ϕ a a
ð14Þ
where xϕ ≈ 0.518 in three spatial dimensions and P(x) is the scaling function. Exactly at the critical point we then have f C ðβc ; D; aÞ a0:518 ≈ 1:518 P ð0Þ: kB T c D
ð15Þ
Although our system is not at the Derjaguin nor at the small sphere limit, we check whether the leading D dependence of the CCFs follows any of the two limiting cases. In Fig. 9(a) we plot the CCFs at the critical point βc multiplied by 0.518 D21/2/(πa) whereas in Fig. 9(b) multiplied by D1.518 for (+[+]+) 1/2 /a boundary conditions and for four different values of a. We observe that for the largest values of a, i.e., for a = 3 and 4, the CCFs multiplied by D21/2/(πa) do not vary much with the distance D for D ≳ 5.5. This observation has encouraged us to apply the “Derjaguinlimit” like rescaling of data obtained for a = 4 and presented in Figs. 6(a), 7(a), and 8(a): 2 D1=2 =ðπaÞ f C β; D1=2 ; a ≡ ϑ D1=2 =ξ :
ð16Þ
Results of this rescaling are plotted in Fig. 6(b), Fig. 7(b), and Fig. 8(b). For (+[+]+) systems we find an acceptable data collapse for larger values of D. 0.518 For the CCFs at the critical point multiplied by D1.518 (see 1/2 /a Fig. 9(b)), we observe that for a = 3 and 4 the CCFs do not vary much with the distance D for D ≳ 9, which is larger than in the previous case. The “small-sphere” like rescaling of data for (+[+]+) systems is shown in Fig. 10. We find a worse data collapse for ``small sphere'' limit than for the “Derjaguin-limit” like rescaling. We note that in the limit of small distance D ≪ a/2 to the wall, the two confining surfaces, i.e., the surface of the particle and the opposing wall, form a film geometry: therefore, for D ≪ a/2 one should simply recover the CCFs for a film geometry. This, however, is impossible in
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
our system because the distances D to which we are limited due to our computational method do not satisfy the condition D ≪ a/2. 4. Conclusion and outlook In this paper we have presented the results of the numerical evaluation of the critical Casimir force acting on the particle placed between two walls for various combinations of the boundary conditions. The coupling parameter approach for extracting the critical Casimir force from MC simulation data has proved to be satisfactory. We have found that the stable mechanical suspension, i.e., the particle not collapsing on either wall, is possible only for the case when neutral particle is placed between the opposite walls ((−[0]+) boundary conditions). However, for (+[+]+) and (−[+]+) boundary conditions a metastable mechanical suspension is possible for temperatures not too close to the bulk critical temperature, where the potential of the critical Casimir force has a flat extremum. Following Derjaguin approximation, we proposed the scaling form for the critical Casimir force. For the (+[+]+) case we have obtained an acceptable data collapse even though the size of the systems which we have been able to study is relatively small. We have studied in detail the cubic particle, but our numerical procedure can be extended to experimentally relevant sphere-plate geometry. Preliminary results for a sphere of radius R = 5.6 in the system size of 50 × 50 × 105 for (+[+]+) boundary conditions as a function of the inverse temperature β are shown in Fig. 11(a) for various distances D 1/2 from the bottom wall. In the lattice system we have modeled the sphere of a radius R as follows. We have selected for the center of the sphere the coordinates of a single spin. The coordinates of spins are integer numbers. Then we have assigned to the sphere all spins which are located closer to the central spin than the sphere radius R = 5.6. The radius is chosen non-integer to provide the right “volume ratio”, i.e., such that the number of spins in the sphere is close to 4/3πR3. We have also plotted our data rescaled according to the formula Eq. (16) which we have used for a cube. As can be seen from Fig. 11(b), one will have to take into account the second scaling variable R/D in the
383
full scaling form (see Eq. (4)) in order to obtain a data collapse. Also applying corrections to scaling will be necessary as well as studying larger systems in order to achieve the accuracy comparable to the one obtained in the MC simulations of the improved Blume–Capel model [16].
References [1] M.E. Fisher, P.G. de Gennes, C. R. Acad. Sci. Paris Ser. B 287 (1978) 207. [2] M.N. Barber, in: C. Domb, J.L. Lebowitz (Eds.), Phase Transitions and Critical Phenomena, vol. 8, Academic, New York, 1983, p. 149; V. Privman, in: V. Privman (Ed.), Finite Size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990, p. 1. [3] R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 1187 (1999); A. Ganshin, S. Scheidemantel, R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 97 (2006) 075301. [4] R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 88 (2002) 086101. [5] M. Fukuto, Y.F. Yano, P.S. Pershan, Phys. Rev. Lett. 94 (2005) 135702. [6] S. Rafaï, D. Bonn, J. Meunier, Physica A 386 (2007) 31. [7] A. Hucht, Phys. Rev. Lett. 99 (2007) 185301. [8] O. Vasilyev, A. Gambassi, A. Maciołek, S. Dietrich, EPL 80 (2007) 60009. [9] O. Vasilyev, A. Gambassi, A. Maciołek, S. Dietrich, Phys. Rev. E. 79 (2009) 041142. [10] O. Vasilyev, A. Maciołek, S. Dietrich, Phys. Rev. E. 84 (2011) 041605. [11] O. Vasilyev, S. Dietrich, EPL 104 (2013) 60002. [12] M. Hasenbusch, Phys. Rev. B 82 (2010) 104425. [13] M. Hasenbusch, Phys. Rev. E. 84 (2011) 041605. [14] M. Hasenbusch, Phys. Rev. B 85 (2012) 174421. [15] K.K. Mon, Phys. Rev. B 39 (1989) 467; K.K. Mon, K. Binder, Phys. Rev. B 42 (1990) 675. [16] M. Hasenbusch, Phys. Rev. E. 87 (2013) 022130. [17] B. Cappella Butt, M. Kappl, Surf. Sci. Rep. 59 (2005) 1–152. [18] J.N. Israelachvili, Intermolecular and Surface Forces, Academic, London, 1992. [19] C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, C. Bechinger, Nature 451 (2008) 172. [20] A. Gambassi, A. Maciolek, C. Hertlein, U. Nellen, L. Helden, C. Bechinger, S. Dietrich, Phys. Rev. E. 80 (2009) 061143. [21] H. Hobrecht, A. Hucht, EPL 106 (2014) 56005. [22] O. Vasilyev, Phys. Rev. E. 90 (2014) 012138. [23] H.B. Casimir, Proc. K. Ned. Akad. Wet. 793 (1948). [24] A. Pelissetto, E. Vicari, Phys. Rep. 368 (2002) 549. [25] C. Ruge, P. Zhu, F. Wagner, Physica A 209 (1994) 431. [26] O.A. Vasilyev, E. Eisenriegler, S. Dietrich, Phys. Rev. E. 88 (2013) 012137. [27] B. Derjaguin, Kolloid Z. 69 (1934) 155. [28] A. Hanke, F. Schlesener, E. Eisenriegler, S. Dietrich, Phys. Rev. Lett. 81 (1998) 1885. [29] F. Schlesener, A. Hanke, S. Dietrich, J. Stat. Phys. 110 (2003) 981.
Contents lists available at ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol
Critical Casimir forces for a particle between two planar walls Oleg Vasilyev a,b, Anna Maciołek a,b,c,⁎ a b c
Max-Planck-Institut für Intelligente Systeme, Heisenbergstr. 3, D-70569 Stuttgart, Germany IV Institut für Theoretische Physik, Universität Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, PL-01-224 Warsaw, Poland
a r t i c l e
i n f o
Article history: Received 27 June 2014 Received in revised form 1 September 2014 Accepted 2 September 2014 Available online 19 September 2014 Keywords: Colloids; Phase transitions and critical phenomena; Classical spin models; Monte Carlo methods
a b s t r a c t We consider a cube placed between two planar surfaces bounding a three-dimensional Ising slab. This system can be viewed as the magnetic analog of a colloidal particle immersed in a fluid, which is confined by two parallel walls. Near the bulk critical point of the Ising ferromagnet the cube gets exposed to the potential of the fluctuation that induces the critical Casimir force. Using Monte Carlo simulations we study the dependence of this potential on the relative boundary conditions at the two surfaces and at the cube within a wide range of temperatures above and below the bulk critical temperature. In order to calculate the critical Casimir force and its potential, we adopt an approach based on an integration scheme of free energy differences. The scaling functions of the corresponding critical Casimir forces are also determined. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Two objects immersed in a fluid which is brought to its bulk critical point Tc experience the so-called critical Casimir interaction [1]. This effective interaction arises due to modifications of the structure of the fluid and due to the restrictions of the fluctuation spectrum of its order parameter by the surfaces of the objects. Accordingly, not only the temperature T but also the adsorption properties of both objects and their geometry determine the sign, the strength and the range of the critical Casimir forces (CCFs). In asymptotic regimes (T → Tc and all distances large compared with microscopic lengths) finite-size scaling [2] holds and the CCFs are described by universal scaling functions. Universality permits to use the simplest possible representative of the same universality class in order to calculate these scaling functions. In this sense the Ising model represents systems such as simple one-component fluids or binary liquid mixtures. The simple geometry of two parallel planes can be realized experimentally by growing a wetting film; exactly the wetting films were used to provide first reliable evidence for the existence of the CCFs [3–6]. On the other hand, even for this simple geometry determination of the scaling functions from corresponding model systems is a challenging task. Theoretical approaches which incorporate critical fluctuations beyond the Gaussian approximation and can take into account a dimensional crossover occurring in films are intractable for most of the relevant experimentally boundary conditions (BCs). Monte Carlo (MC) simulations offer a very useful alternative approach, ⁎ Corresponding author at: Max-Planck-Institut für Intelligente Systeme, Heisenbergstr. 3, D-70569 Stuttgart, Germany. E-mail address: [email protected] (A. Maciołek).
http://dx.doi.org/10.1016/j.jnoncrysol.2014.09.005 0022-3093/© 2014 Elsevier B.V. All rights reserved.
which is however not free from difficulties. Within this approach the range of sizes of the system is strongly limited by the steeply increasing computational costs. For the thicknesses and the aspect ratios of slabs accessible in simulations, numerical data do not collapse and corrections to scaling are necessary. Renormalization-group analyses reveal that there is a whole variety of sources for corrections to scaling which arise from bulk, surface, and finite-size effects [2]. In addition to the leading bulk corrections to scaling (with the correction-to-scaling exponent ω ≃ − 0.8) one expects corrections due to the finite aspect ratio, corrections due to the boundary conditions (which are especially large for the symmetry breaking boundary conditions) and, moreover, nextto-leading corrections might occur for narrow films. In spite of these difficulties the Casimir scaling functions in three spatial dimensions (3d) have been obtained rather accurately for experimentally relevant universality classes with a variety of BCs and surface fields [7–14]. In the present paper we study numerically the CCFs acting on a cube located between two planar walls as a function of temperature and separation for different BCs. The cube mimics a colloidal particle positioned near the wall. We use MC simulations of the Ising model on a cubic lattice as a representative of the Ising universality class of critical phenomena. The aim of this work is twofold. We would like to check whether the approach that we used in our previous studies for a slab geometry [8–11] can be successfully adopted for the present geometry. This approach is based on the so-called coupling parameter method where the difference of free energies is expressed as an integral over the mean energy difference (see, e.g., Refs. [8] and [15]). Moreover, we address the question of the mechanical stability of a cube particle for the cases when (a) both walls and the particle have the same adsorption preference (+| + |+), (b) the walls have opposite adsorption preferences (+| + |−), and (c) the cube particle is
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 1. The geometry under consideration (the cross-section in the xz plane): (a) the spherical particle of the diameter a at the distance D from the wall; (b) an appropriate geometry on the lattice for a cubic particle with the edge length a at the separation D from the bottom wall. The total system size is Lx × Ly × (Lz + 2), and the separation of the particle from the top wall is Lz − D − a, where Lz is the number of layers of free spins.
377
suppressed. The focus of that study was on symmetry breaking boundary conditions at the surfaces of the sphere and the plate. The use of the geometric cluster algorithm and the introduction of an effective radius of the sphere and the effective distance from the wall have allowed the author to obtain a good data collapse. Still another MC simulation method to compute the critical Casimir force acting between the disk and the wall in two-dimensional Ising systems has been proposed in Ref. [21]. This method is analogous to the one used in experiments reported in Refs. [19,20] for a colloidal particle where the Casimir potential is determined from the distribution of the position of the particle above the wall. The interaction between two particles in the presence of the bulk ordering field recently has been studied by using a local field integration method [22]. Our presentation is organized as follows. In Section 2 we briefly present the relevant theoretical background, introduce our model, and describe the numerical method employed in order to compute the CCFs and its scaling functions from the MC simulation data. Section 3 contains our results. We provide a summary and conclusions in Section 4. 2. The model and the method
neutral and the walls have opposite adsorption preferences (+|0|−). These situations can be mimicked in our Ising system by choosing an appropriate BC on the lattice sites forming boundaries of the three objects. Because of unavoidable limitations of the size of the system, we cannot expect to reach the asymptotic region of the universal behavior in our simulations. Nevertheless, after rescaling procedure for some BCs we have obtained a reasonable data collapse. Recently, a different MC simulation method for spin models which mimic the sphere–planar plane geometry has been developed [16]. The geometry of a sphere or a cube near a single planar surface is relevant for colloidal systems or in experimental setups of high-precision force measuring devices and methods, such as atomic force microscope (AFM) [17], surface force apparatus (SFA) [18] or total internal reflection microscopy (TIRM) [19,20]. For this geometry theoretical predictions for the CCFs are much more limited than for films. The MC simulation method employed in Ref. [16] is similar to the one used in [7] for the plate–plate geometry, where differences of the free energy are computed by integrating differences of the energy over the inverse temperature. Simulations have been performed for the improved Blume–Capel model on the simple cubic lattice, which shares the universality class of the three-dimensional Ising model but has the advantage over it that the amplitudes of the leading correction to scaling are substantially
(a)
2.1. Theoretical background Let us consider a spherical colloidal particle of a diameter a immersed in critical binary mixture near a wall, see Fig. 1(a). As already mentioned in the Introduction, at the critical point of the demixing transition of the binary solvent, the CCFs, which are the analogs of the electromagnetic Casimir force [23], emerge between the wall and the colloidal particle as a result of critical fluctuations of the order parameter of the solvent and the critical adsorption of a preferable component of the mixture on the surfaces of wall and the colloid. For a binary mixture the order parameter is the deviation of the concentration of the one component of the mixture from its value at the critical point. The CCFs between the wall and the particle are obtained from the free energy UC(β, D, a) of a system consisting of the colloidal particle of the characteristic linear size a immersed in the fluid at the inverse temperature β = 1/kBT at the separation D from the planar wall: f C ðβ; D; aÞ ¼ −
∂U C ðβ; D; aÞ : ∂D
ð1Þ
Note that the above definition of the CCFs differs from the one used to define the CCFs between two planar parallel plates immersed in the
(b)
(c)
Fig. 2. Bond arrangement for the computation of the free energy difference in Eq. (9) (see main text). The crossover Hamiltonian Hcr (a) describes the system which interpolates between those described by the Hamiltonians H0 (b) and H1 (c).
378
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 3. Results for (+[+]+) boundary conditions as a function of the distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The cube has the edge length of a = 4 (in units of lattice spacing) whereas the slab has the size of 42 × 42 × 45.
fluid. If D is the separation between the two plates then for large areas A of plates the total free energy F(β, D, A) of the system can be written as F ðβ; D; AÞ ≡ ADf ðβ; DÞ h i bulk −1 ex ¼ A Df ðβÞ þ β f ðβ; DÞ ;
ð2Þ
where f bulk(β) is the bulk free energy density at a given temperature. The excess free energy f ex per area contains two D-independent surface contributions in addition to the finite-size contribution f ex(β, D) − f ex(β, ∞) which vanishes for D → ∞. The D-dependence of the latter gives rise to the critical Casimir force fC per unit area A: ex
f C ðβ; DÞ=kB T ≡ −∂f ðβ; DÞ=∂D:
below (−) the bulk critical temperature Tc with ξ− 0 = 0.243(1) and ξ+ 0 = 0.501(2) [25] for the Ising model on the simple cubic lattice. The value of ξ+ 0 quoted here refers to the amplitude of the second moment correlation length ξ2nd ; however, ξ=ξ2nd ≃1 for β b βc for the Ising model [24]. The scaling function K also depends on the BCs imposed by the particle and the wall on the order parameter, which reflects the adsorption preferences of the surfaces of these objects for one of the two species of the mixture. 2.2. The model We consider the Ising model defined on a three-dimensional simple cubic lattice with the lattice spacing ℓ = 1 via the Hamiltonian
ð3Þ H ¼ −J
From the general theory of finite-size scaling [2] it follows that the CCFs between the wall and the particle can be cast in the following scaling form f C ðβ; D; aÞ≈
kB T D D K ; a ξ a
ð4Þ
where the universal scaling function K depends on the ratio D/ξ and also on the ratio D/a. In the critical region and vanishing ordering bulk field −ν h = 0, the bulk correlation length is ξ = ξ± , where the reduced 0 |τ| temperature is defined as τ = (T − Tc)/Tc = (βc − β)/β. The critical exponent takes the value of ν ≈ 0.63 for the Ising bulk universality class in 3d [12,24], ξ± 0 are nonuniversal amplitudes above (+) and
(a)
X
si s j ;
ð5Þ
hi; ji
where si = ± 1 denotes the spin variable. J, which we set equal to 1, is the spin–spin coupling constant and the sum 〈i, j〉 is taken over all nearest neighbor pairs of sites i and j on the lattice. Temperatures and energies are measured in units of J. The inverse critical temperature is βc = 0.2216544(3) [25]. We consider a slab geometry, i.e., a lattice of fluctuating spins of the size Lx × Ly × Lz with periodic BCs in the x and y directions (in which the system has linear extensions Lx and Ly). There are two additional top and bottom layers as shown in Fig. 1(b). All spins in the top (st, z = Lz + 1) layer are fixed at the value + 1 whereas all spins in the bottom (sb, z = 0) layer are fixed at the value + 1 or −1 for different types of BCs. We have modified the coupling constants
(b)
Fig. 4. Results for (−[+]+) boundary conditions as a function of the distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The system size and the edge length of a cube are the same as in Fig. 3.
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
379
(b)
Fig. 5. Results for (−[0]+) boundary conditions as a function of distance D for different values of β: (a) the critical Casimir force fC and (b) the potential of the Casimir force FC. The system size and the edge length of a cube is the same as in Fig. 3.
between fixed spins at the top and bottom boundaries and those in the adjacent layers. In our present simulations we consider the infinite strength of these modified bonds, which correspond to the infinite surface fields acting on the next-to-surface layers of the system. We also fix a3 spins sp of the cubic particle of linear size a. The spins sp can take values sp = ± 1(for [+] and [−]) and 0 (for [0]). The 0 value corresponds to the neutral particle that does not interact with the components of the binary mixture. As for the walls, we have modified the bonds between the spins at the surface of a cube and their nearest neighbors and assume that their strength is infinite. We denote the type of BCs as (sb[sp]st) (as it was done in Ref. [26] for computation of torques in 2d systems), for example, (−[0]+) BCs correspond to negative spins sb = −1 at the bottom layer, zero spins sp = 0 for the particle and positive spins st = +1 for the top layer. In the present paper we consider three cases: (i) a particle between two attractive walls (+[+]+), (ii) a particle between the repulsive and attractive walls (−[+]+), and (iii) a neutral particle between walls of opposite signs (−[0]+).
2.3. Computation of free energy differences In general, the Monte Carlo methods are efficient for the computation of quantities which can be expressed as ensemble averages, such as energy, magnetization, specific heat or magnetic susceptibility of a system. The free energy differences can be nevertheless computed
(a)
accurately by using the so-called “coupling parameter approach” [15]. Recently, this method was used for the computation of the Casimir force scaling function for the slab geometry [8–11]. For the lattice with spacing ℓ = 1 we rewrite Eq. (1) as the finite difference: f C ðβ; D−1=2; aÞ ¼ −ΔU C ðβ; D; aÞ ¼ −½U C ðβ; D; aÞ−U C ðβ; D−1; aÞ ¼ −½ F 0 −F 1 :
ð6Þ The free energy UC(β, D, a) = F0 corresponds to the system where the cubic particle of the linear size a is at the separation D from the bottom wall (the wall is assumed to be located between the layer at z = 0 and the layer at z = 1) whereas the free energy UC(β, D − 1, a) = F1 corresponds to the similar system but with the separation D − 1. See Fig. 1 where distances are defined. In order to compute the free energy difference in Eq. (6) we introduce the crossover Hamiltonian H cr ðλÞ ¼ ð1−λÞH 0 þ λH 1 ;
ð7Þ
where λ ∈ [0, 1] is the crossover parameter. This Hamiltonian describes the system consisting of Lz + 1 layers of free spins in which the coupling constant associated with certain bonds is equal to λ and with some other bonds is equal to 1 − λ, as depicted in Fig. 2(a). These bonds of different strength are arranged around two selected horizontal layers of spins. The bottom selected layer is located at zbl = D − 2 (at distance 2 below the
(b)
Fig. 6. Results for (+[+]+) boundary conditions for different distances D: (a) the Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
380
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 7. Results for (−[+]+) boundary conditions for different distances D: (a) the critical Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
particle) and the top selected layer is located at ztl = D + a + 3 (at distance 2 above the particle), as shown in Fig. 2(a). In this figure, fixed spins of walls and of the particle are denoted by filled circles, while free (fluctuating) spins are denoted by open circles. For λ = 0, the bottom selected layer is incorporated into the bulk of the system and the top selected layer is excluded from interactions as shown in Fig. 2(b). In this configuration, which is described by the Hamiltonian Hcr(λ = 0) = H0 corresponding to the system with free energy F0 plus the free energy of isolated 2D layer, the bottom face of the cubic particle is at the distance D (D layers of free spins) from the bottom wall, and the upper face is at the distance Lz − D − a from the top wall. For λ = 1 the bottom selected layer is excluded from interactions and the top selected layer is incorporated into the bulk as shown in Fig. 2(c). The total system height is Lz (as for the case λ = 0). In this configuration, which is described by the Hamiltonian Hcr(λ = 1) = H1 corresponding to the system with free energy F1 plus the free energy of isolated 2D layer, the bottom face of the cubic particle is at the distance D − 1 from the bottom wall and the top face is at the distance Lz − D − a + 1 from the top wall. Locations of the selected layers are rather arbitrary; the removed layer should be somewhere in between the particle and the bottom wall and the inserted layer should be somewhere in between the particle and the top wall. As a function of the coupling parameter λ, Hcr(λ) interpolates between H0 and H1 as λ increases from 0 to 1, but the configurational space (amount and arrangement of free spins) remains the same. Only the interaction strength of some bond is changed. The free energy F cr ðλÞ ¼ −β1 ln ∑C expð−βHcr ðλÞÞ, where the sum is taken over all configurations fCg of free spins and interpolates between F0 and F1. The free
(a)
energy difference F1 − F0 can be expressed as F1 − F0 = ∫ 10Fcr′(λ)dλ where Fcr′ is the derivative of Fcr(λ) with respect to the coupling parameter λ: dF cr ðλÞ ¼ dλ
X C
ðH 1 −H 0 Þe−βHcr X e−βHcr ðλÞ C
ðλÞ
¼ hΔH ðλÞifCg ;
ð8Þ
which takes the form of the canonical ensemble average h…ifCg of ΔH(λ) ≡ H1 − H0 over all configurations fCg for the fixed value of the coupling parameter λ and therefore it can be computed via MC simulations. Actually, ΔH(λ) equals to the thermal average of the interaction energy of pair of spins connected by appearing bonds minus the interaction energy of spins connected by vanishing bonds in Fig. 2(a) for a given value of λ. As a result one can express the difference in free energies as an integral over canonical averages (see, e.g., Ref. [15]): Z ΔF ¼ F 1 − F 0 ¼
1 0
hΔH ðλÞifCg dλ:
ð9Þ
In that formula contributions of 2D layers for F0 and F1 are mutually canceled. Finally we obtain the expression for the Casimir force for separation D − 1/2 as Z f C ðβ; D−1=2; aÞ ¼ −
1 0
hΔH ðλÞifCg dλ:
ð10Þ
Here, the force as the finite difference between separations D and D − 1 is assigned to separation D − 1/2. Let us note that due to finite size
(b)
Fig. 8. Results for (−[0]+) boundary conditions for different distances D: (a) the critical Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(πa))fC as a function of (D − 1/2)/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
381
(b)
0.518 Fig. 9. The critical Casimir force at the critical point Tc multiplied by: (a) D21/2/(πa) and (b) D1.518 for different a and D for (+[+]+) boundary conditions. 1/2 /a
of our system the particle interacts also with the upper wall, which is included in Eq. (10). From the knowledge of the critical Casimir force f C ðβ; D−1=2; aÞ ¼ −∂U ¼ U C ðβ; D−1; aÞ−U C ðβ; D; aÞ , we can compute the potential of ∂D the CCFs UC as a function of a distance from the wall D with respect to some given separation D0. We select as a reference plane from which the potential is measured the midplane of the system D0 = Lz/2. It means that the potential of the CCFs is zero at this plane. For a lattice model (a discrete difference instead of derivative), the potential is computed in accordance to the formula C
U C ðβ; D; aÞ ¼
8 D X > > > > f C ðβ; z−1=2; aÞ; − > < > > > > > :
z¼D0 þ1 D−1 X
f C ðβ; z−1=2; aÞ;
D≥D0 :
ð11Þ
DbD0
z¼D0
3. Numerical results 3.1. Force and its potential as functions of the distance from the wall We compute the critical Casimir force fC (β, D − 1/2, a) and its potential UC(β, D, a) for a cubic particle of edge length a = 4 placed in the system of 42 × 42 × 45 free spins. Results as a function of the separation D − 1/2 for different values of the inverse temperature β are plotted in Fig. 3(a) for (+[+]+), in Fig. 4(a) for (−[+]+), and in Fig. 5(a) for (−[0]+). A potential UC(β, D, a) of the CCFs for these cases is presented in Fig. 3(b), Fig. 4(b), and Fig. 5(b), respectively. In all presented plots we have adopted the notation fC(β, D1/2, a) for fC(β, D − 1/2, a) and D1/2 for D − 1/2. (Note, that because of our choice of the selected layers for the computation of the free energy difference, the distance at which we can calculate the CCFs varies between 4.5 and 35.5). As can be expected, for (+[+]+) boundary conditions the cube is attracted to both walls equally strong so that the force is vanishing in the middle of the slab and its potential has a maximum there. Because the force is measured as a function of a distance from the bottom wall, it becomes positive (repulsive) for D1/2 N D0 = Lz/2 − 1/2 as a result of dominating attraction to the top wall. Interestingly, for temperatures sufficiently far below or above the critical temperature (the blue and the red curve in Fig. 3(a)), the maximum of the potential is very flat. This indicates that the particle can be positioned in the middle part of a slab and stay there for a relatively long time. The potential wells at the walls are the deepest at some temperature slightly above Tc. These features are consistent with the ones of the CCFs between two planar (+,+) surfaces [8], which are attractive, possess a minimum above Tc,
and decay exponentially fast with the deviation from the critical temperature τ. For (−[+]+) boundary conditions the symmetry is broken; the cube is attracted to the (+) wall and repelled from the (−) wall. Below Tc, the repulsion from the (−) wall is much stronger then the attraction to the (+) wall. This is due to the presence of the interface between predominantly positive and negative magnetized phases located between the (−) wall and the cube. The resulting potential decays rapidly and monotonically from the large positive values at the (−) wall towards its minimum located at the (+) wall. The depth of the minimum is comparable with that for the (+[+]+) case. Above Tc, the interface vanishes and the repulsion from the (−) wall, as we know it from the case of two planar walls with (−,+) boundary conditions, becomes very weak. In contrast, the attraction from the (+) walls attains its extreme value above Tc. The net effect for the cube is that the force acting on it is negligibly small around the midpoint of a slab and its potential is very flat. The CCFs and its potential are asymmetric. In the case of the “neutral” cube confined between the opposite walls, (−[0]+) boundary conditions have been numerically the most challenging one — notice the big error bars in Fig. 5. This is because the order parameter at the cube is very small and the resulting CCFs are relatively weak as compared to the previous cases. Recall that for the planar walls with (0,+) (or equivalently (0,−) boundary conditions [10]), the CCFs are repulsive with a maximum below Tc — just as for the (−,+) case, however, the amplitude at its maximum is about 10 times smaller than for the (+,−) boundary conditions. The neutral cube is repelled equally strong from both walls. Similar to the (−[+]+) case, below Tc the boundary conditions at the walls create an
Fig. 10. Results for (+[+]+) boundary conditions for different distances D: the rescaled Casimir force D1.518 1/2 fC as a function of D1/2/ξ. The system size and the edge length of a cube are the same as in Fig. 3.
382
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
(a)
(b)
Fig. 11. Results for a sphere of radius R = 5.6 placed near the wall with the same boundary conditions (surface spins fixed at the value of + 1) for different distances D: (a) the Casimir force fC as a function of β and (b) the rescaled force ϑ ≡ (D21/2/(2πR))fC as a function of D1/2/ξ. The system size is 50 × 50 × 105.
interface between oppositely magnetized phases which is distorted by a presence of the cube. However, for the neutral cube this distortion is symmetric. As a result, the critical Casimir force varies approximately linearly with the distance across the slab, except for the immediate vicinity of the walls, i.e., one or two lattice sites. The slope of this linear variation decreases as the temperature approaches the critical temperature. At T = Tc, the linear behavior is narrowed to the middle part of the slab. Accordingly, the critical Casimir potential is quadratic and strongly localizes the particle in the middle of the slab. Above Tc, the potential is very flat so that the cube can take stable positions within relatively large region near the center of the system.
e ðxÞ ¼ 2π∫∞ d‘ 1=‘3 θð‘xÞ, i.e., it does not where the scaling function K 1 depend on a. Here, θ(x) is the Casimir force scaling function for the planar plane (film) geometry. At the critical point β = βc one has f C ðβc ; D; aÞ πa ≈ 2 θð0Þ: kB T c D
In the opposite limit of D ≫ a/2, the form of the scaling for the CCFs follows from the small sphere expansion [28] f C ðβ; D; aÞ≈
3.2. The force as a function of temperature In order to study the temperature dependence of the CCFs acting on a cube, we compute fC(β, D1/2, a) for various values of inverse temperature β at several separations D1/2 = 4.5, 5.5, …, 19.5. Results are plotted in Fig. 6(a) for (+[+]+), Fig. 7(a) for (−[+]+) and Fig. 8(a) for (−[0]+). For the (+[+]+) case, the features of the CCFs are similar to the features of CCFs for the planar plane geometry with fixed (+,+) boundary conditions, i.e., the force is attractive and has its minimum above Tc. The same holds for (−[+]+) and (−[0]+) systems if compared to the planar plane geometry with (−,+) and (−, 0) boundary conditions, respectively. In these cases the force is attractive with a maximum above T c. Data for (−[0]+) systems suffer from large computation errors in the region of low temperatures. 3.3. Scaling function of the critical Casimir force Because the size of the systems for which we were able to perform our simulations is rather limited, we cannot expect to determine the Casimir force scaling function K Dξ ; Da . The additional difficulty is the presence of the second scaling variable y = D/a. We could plot the data as a function of the first variable while keeping the second one fixed. This, however, would require to consider a cube with different size a for each distance D, which is computationally very costly. What we try to do in the present paper is to investigate the leading dependence of the CCFs on D and on a. As a guide we use results for a spherical particle of radius a/2 located at a distance D of closest approach from the flat surface of a substrate. The CCFs acting on the particle in the limit of D ≪ a/2 is given by the Derjaguin approximation as [20,27–29]: f C ðβ; D; aÞ a=2 e D ≈ 2 K kB T ξ D
ð12Þ
ð13Þ
kB T D −xϕ −1 −2x −1 P ðD=aÞ þ O ðD=aÞ ϕ a a
ð14Þ
where xϕ ≈ 0.518 in three spatial dimensions and P(x) is the scaling function. Exactly at the critical point we then have f C ðβc ; D; aÞ a0:518 ≈ 1:518 P ð0Þ: kB T c D
ð15Þ
Although our system is not at the Derjaguin nor at the small sphere limit, we check whether the leading D dependence of the CCFs follows any of the two limiting cases. In Fig. 9(a) we plot the CCFs at the critical point βc multiplied by 0.518 D21/2/(πa) whereas in Fig. 9(b) multiplied by D1.518 for (+[+]+) 1/2 /a boundary conditions and for four different values of a. We observe that for the largest values of a, i.e., for a = 3 and 4, the CCFs multiplied by D21/2/(πa) do not vary much with the distance D for D ≳ 5.5. This observation has encouraged us to apply the “Derjaguinlimit” like rescaling of data obtained for a = 4 and presented in Figs. 6(a), 7(a), and 8(a): 2 D1=2 =ðπaÞ f C β; D1=2 ; a ≡ ϑ D1=2 =ξ :
ð16Þ
Results of this rescaling are plotted in Fig. 6(b), Fig. 7(b), and Fig. 8(b). For (+[+]+) systems we find an acceptable data collapse for larger values of D. 0.518 For the CCFs at the critical point multiplied by D1.518 (see 1/2 /a Fig. 9(b)), we observe that for a = 3 and 4 the CCFs do not vary much with the distance D for D ≳ 9, which is larger than in the previous case. The “small-sphere” like rescaling of data for (+[+]+) systems is shown in Fig. 10. We find a worse data collapse for ``small sphere'' limit than for the “Derjaguin-limit” like rescaling. We note that in the limit of small distance D ≪ a/2 to the wall, the two confining surfaces, i.e., the surface of the particle and the opposing wall, form a film geometry: therefore, for D ≪ a/2 one should simply recover the CCFs for a film geometry. This, however, is impossible in
O. Vasilyev, A. Maciołek / Journal of Non-Crystalline Solids 407 (2015) 376–383
our system because the distances D to which we are limited due to our computational method do not satisfy the condition D ≪ a/2. 4. Conclusion and outlook In this paper we have presented the results of the numerical evaluation of the critical Casimir force acting on the particle placed between two walls for various combinations of the boundary conditions. The coupling parameter approach for extracting the critical Casimir force from MC simulation data has proved to be satisfactory. We have found that the stable mechanical suspension, i.e., the particle not collapsing on either wall, is possible only for the case when neutral particle is placed between the opposite walls ((−[0]+) boundary conditions). However, for (+[+]+) and (−[+]+) boundary conditions a metastable mechanical suspension is possible for temperatures not too close to the bulk critical temperature, where the potential of the critical Casimir force has a flat extremum. Following Derjaguin approximation, we proposed the scaling form for the critical Casimir force. For the (+[+]+) case we have obtained an acceptable data collapse even though the size of the systems which we have been able to study is relatively small. We have studied in detail the cubic particle, but our numerical procedure can be extended to experimentally relevant sphere-plate geometry. Preliminary results for a sphere of radius R = 5.6 in the system size of 50 × 50 × 105 for (+[+]+) boundary conditions as a function of the inverse temperature β are shown in Fig. 11(a) for various distances D 1/2 from the bottom wall. In the lattice system we have modeled the sphere of a radius R as follows. We have selected for the center of the sphere the coordinates of a single spin. The coordinates of spins are integer numbers. Then we have assigned to the sphere all spins which are located closer to the central spin than the sphere radius R = 5.6. The radius is chosen non-integer to provide the right “volume ratio”, i.e., such that the number of spins in the sphere is close to 4/3πR3. We have also plotted our data rescaled according to the formula Eq. (16) which we have used for a cube. As can be seen from Fig. 11(b), one will have to take into account the second scaling variable R/D in the
383
full scaling form (see Eq. (4)) in order to obtain a data collapse. Also applying corrections to scaling will be necessary as well as studying larger systems in order to achieve the accuracy comparable to the one obtained in the MC simulations of the improved Blume–Capel model [16].
References [1] M.E. Fisher, P.G. de Gennes, C. R. Acad. Sci. Paris Ser. B 287 (1978) 207. [2] M.N. Barber, in: C. Domb, J.L. Lebowitz (Eds.), Phase Transitions and Critical Phenomena, vol. 8, Academic, New York, 1983, p. 149; V. Privman, in: V. Privman (Ed.), Finite Size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990, p. 1. [3] R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 1187 (1999); A. Ganshin, S. Scheidemantel, R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 97 (2006) 075301. [4] R. Garcia, M.H.W. Chan, Phys. Rev. Lett. 88 (2002) 086101. [5] M. Fukuto, Y.F. Yano, P.S. Pershan, Phys. Rev. Lett. 94 (2005) 135702. [6] S. Rafaï, D. Bonn, J. Meunier, Physica A 386 (2007) 31. [7] A. Hucht, Phys. Rev. Lett. 99 (2007) 185301. [8] O. Vasilyev, A. Gambassi, A. Maciołek, S. Dietrich, EPL 80 (2007) 60009. [9] O. Vasilyev, A. Gambassi, A. Maciołek, S. Dietrich, Phys. Rev. E. 79 (2009) 041142. [10] O. Vasilyev, A. Maciołek, S. Dietrich, Phys. Rev. E. 84 (2011) 041605. [11] O. Vasilyev, S. Dietrich, EPL 104 (2013) 60002. [12] M. Hasenbusch, Phys. Rev. B 82 (2010) 104425. [13] M. Hasenbusch, Phys. Rev. E. 84 (2011) 041605. [14] M. Hasenbusch, Phys. Rev. B 85 (2012) 174421. [15] K.K. Mon, Phys. Rev. B 39 (1989) 467; K.K. Mon, K. Binder, Phys. Rev. B 42 (1990) 675. [16] M. Hasenbusch, Phys. Rev. E. 87 (2013) 022130. [17] B. Cappella Butt, M. Kappl, Surf. Sci. Rep. 59 (2005) 1–152. [18] J.N. Israelachvili, Intermolecular and Surface Forces, Academic, London, 1992. [19] C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, C. Bechinger, Nature 451 (2008) 172. [20] A. Gambassi, A. Maciolek, C. Hertlein, U. Nellen, L. Helden, C. Bechinger, S. Dietrich, Phys. Rev. E. 80 (2009) 061143. [21] H. Hobrecht, A. Hucht, EPL 106 (2014) 56005. [22] O. Vasilyev, Phys. Rev. E. 90 (2014) 012138. [23] H.B. Casimir, Proc. K. Ned. Akad. Wet. 793 (1948). [24] A. Pelissetto, E. Vicari, Phys. Rep. 368 (2002) 549. [25] C. Ruge, P. Zhu, F. Wagner, Physica A 209 (1994) 431. [26] O.A. Vasilyev, E. Eisenriegler, S. Dietrich, Phys. Rev. E. 88 (2013) 012137. [27] B. Derjaguin, Kolloid Z. 69 (1934) 155. [28] A. Hanke, F. Schlesener, E. Eisenriegler, S. Dietrich, Phys. Rev. Lett. 81 (1998) 1885. [29] F. Schlesener, A. Hanke, S. Dietrich, J. Stat. Phys. 110 (2003) 981.