text:
# Cut-and-permute algorithm for self-avoiding walks in the presence of surfaces

## 1 Introduction

In this paper we introduce and discuss the properties of a hybrid Monte Carlo algorithm which can be used to study the equilibrium properties of a polymer molecule grafted to a surface.

The system can be experimentally obtained either by chemically grafting one polymer end in an irreversible way, or, as for surfactants, via physical adsorption of an endgroup or one of the two blocks in a diblock copolymer. In this last case, the process is reversible, since the attached end can desorb both when the temperature increases and when the solvent quality changes.

Systems of this kind include polymers grafted at colloidal particles or surfaces in solution which can help stabilizing against flocculation  or polymers grafted at flexible lipid membranes exerting on the membrane a bending force which is proportional to the monomer concentration at the membrane .

In what follows we will consider a polymer in a good solvent which is grafted at a surface $`S`$ and interacts repulsively with $`S`$. We will focus in particular on the case in which $`S`$ is a flat surface. The polymer will be modelled by an $`N`$-step lattice self-avoiding walk (SAW), which provides a good description of the critical behavior of polymer molecules in the bulk as well as in confined geometries .

Efficient simulations of SAWs can be obtained by using nonlocal algorithms. For instance, the pivot algorithm (see Refs.  and references therein) is optimal, up to a constant factor, for sampling global observables in the fixed-$`N`$ free-endpoint ensemble in the bulk since the autocorrelation time in CPU units is simply proportional to $`N`$. In the presence of surfaces, this algorithm is not as efficient and in some cases it is not even ergodic (this is the case of a two-dimensional strip, see Ref. ). In the presence of a single impenetrable plane, the algorithm is ergodic, but still one expects it to be inefficient since the initial part of the walk will be rarely updated. For instance, we will show that, in order to update the direction of one of the first links pointing along the normal vector to the surface, the walk should have an extremely unlikely geometrical shape. This means that, in order to explore the relevant phase space, the algorithm has to go through highly improbable configurations. Therefore, even if ergodic, the algorithm has an autocorrelation time that grows rapidly with $`N`$.

To overcome these problems, we propose here a hybrid algorithm based on the pivot move and on a cut-and-permute move. This last move has been introduced in Ref.  in a hybrid algorithm working in the fixed-endpoint, variable-$`N`$ ensemble. It consists in cutting the walk into two parts and in rebuilding it by concatenating the two parts in the opposite order. As we will discuss, such a move is quite efficient in changing the configuration of the walk near the grafted endpoint. Thus, if we combine the pivot move and the cut-and-permute move, we obtain an algorithm which does not have the quasi-ergodicity problems of the pivot algorithm. Also, by a careful implementation, it is possible to obtain the same optimal behavior as the pivot algorithm for polymers in the bulk: the autocorrelation times of global observables in CPU units increase only as $`N`$. It must be noted that, even though we consider here only the case of a flat surface, the results should apply to any convex surface, for instance to polymers grafted outside a sphere: in this respect, the plane should be the worst case. Note also that we do not consider here any interaction between the surface and the walk other than the excluded-volume interaction. Nonetheless, the results should also hold in the presence of attractive interactions as long as they are sufficiently weak and the walk is not absorbed. In this case, it is probably important to add local moves (and, perhaps, bilocal moves as defined in Ref. ) to speed up the dynamics near the surface. Finally, we want to notice that this algorithm is also needed if we want to apply the join-and-cut algorithm of Ref.  to walks in the presence of a surface. Indeed, the ideas that are presented here apply directly to that algorithm, so that by using the cut-and-permute together with the pivot move, one should be able to have a version of the join-and-cut algorithm which works reasonably well also in the presence of a surface.

The paper is organized as follows. In Sec. 2 and 3 we consider the pivot and the cut-and-permute move and discuss in detail the acceptance fraction, paying particular attention to those moves in which the pivot or the cutting point is near the grafted endpoint of the walk. In Sec. 4 we discuss the implementation of the two moves and compute the scaling behavior of the average CPU time for each move. In Sec. 5 we discuss the full algorithm, showing its optimal behavior.

## 2 The Pivot move in the presence of a surface

In this paper we will consider $`N`$-step self-avoiding walks on a $`d`$-dimensional lattice $`^d`$ in the presence of an excluded surface $`S`$ of equation $`S(\stackrel{}{x})=0`$, $`\stackrel{}{x}^d`$. A SAW $`\omega `$ is a sequence of lattice points $`\{\omega _0,\mathrm{},\omega _N\}`$, such that $`\omega _i`$, $`\omega _{i+1}`$ are lattice neighbours and $`\omega _i\omega _j`$ for $`ij`$. The walk is confined in the outward half-space $`S(\stackrel{}{x})0`$ with its first vertex grafted at $`S`$ and fixed at position $`\omega _0`$. The state space is therefore

$$\mathrm{\Omega }_N^\mathrm{S}(\omega _0)=\{\omega =\{\omega _0,\mathrm{},\omega _N\}|S\left(\omega _i\right)0,i=1,\mathrm{},N\mathrm{and}\omega _i\omega _j,ij\}.$$
(1)
As probability measure on the space $`\mathrm{\Omega }_N^\mathrm{S}(\omega _0)`$, we will consider the uniform one, which gives equal weight to every allowed walk, i.e. $`\pi \left(\omega \right)=\frac{1}{c_N^\mathrm{S}}`$, where $`c_N^\mathrm{S}`$ is the cardinality of $`\mathrm{\Omega }_N^\mathrm{S}(\omega _0)`$. In the following we will assume $`S`$ to be a $`(d1)`$ hyperplane of equation $`z=0`$.

In $`\mathrm{\Omega }_N^S(\omega _0)`$ we will first consider the pivot algorithm . The elementary move of the algorithm in the bulk consists of the following steps:

1. Given a walk configuration $`\omega =\{\omega _0,\mathrm{},\omega _N\}`$, an integer $`k\{0,\mathrm{},N1\}`$ is chosen at random and the corresponding monomer $`\omega _k`$ is taken as pivot point.
2. An element $`g`$ of the symmetry group of the lattice is chosen at random with probability $`p_g`$ and a new walk is built by applying it to the part of the chain which follows the pivot point, using $`\omega _k`$ as fixed point of the transformation. A walk $`\omega ^{}=\{\omega _0^{},\mathrm{},\omega _N^{}\}`$, where the new monomer coordinates are given by

$$\omega _{}^{}{}_{i}{}^{}=\{\begin{array}{cc}\omega _i\hfill & \text{for }1ik\text{,}\hfill \\ \omega _k+g(\omega _i\omega _k)\hfill & \text{for }k+1iN\text{ ,}\hfill \end{array}$$
(2)
is obtained.
3. If the walk $`\omega ^{}`$ is not self-avoiding the move is rejected and the original walk $`\omega `$ is kept. Otherwise, $`\omega ^{}`$ is taken as new current walk.

We require $`p_g=p_{g^1}`$ in order to satisfy detailed balance. In the bulk the algorithm is ergodic if all axes reflections and either all $`\pi /2`$ rotations or all diagonal reflections are given nonzero probability . Note that in the bulk, in step 1, we can restrict $`k`$ to belong to $`\{1,\mathrm{},N1\}`$, since pivot moves with pivot $`\omega _0`$ are symmetry transformations. However, this is not the case in the presence of a surface, and it is thus important to include $`k=0`$.

A further requirement has to be satisfied in order to take into account the excluded region:

4. If $`\omega ^{}\mathrm{\Omega }_N^S(\omega _0)`$, i.e. if $`S\left(\omega _i^{}\right)<0`$ for some $`i=1,\mathrm{},N`$, the move is rejected and the old walk is counted again.

Because of the presence of the excluded region, ergodicity is not always satisfied. If $`S`$ is a $`(d1)`$-dimensional hyperplane, one can prove that axes reflections and $`\pi /2`$ rotations are enough to ensure ergodicity. However, even if ergodicity is satisfied, one may be worried by the fact that in order to explore the relevant phase space, the algorithm has to go through highly improbable configurations. For instance, if we denote by $`\widehat{z}`$ the direction perpendicular to the hyperplane $`S`$, we expect the algorithm to be quite slow in updating the links at the beginning of the walk that are directed in the $`\widehat{z}`$ direction. Indeed, transformations that have the pivot near the surface and that modify the $`\widehat{z}`$ direction—moves involving the inversion of the $`z`$-axis, reflections with respect to the diagonals, $`\pi /2`$ rotations in the $`x=0`$ or $`y=0`$ planes—are unlikely to be accepted because of the excluded region.

Therefore, it is important to study the dynamic behaviour of the algorithm as a function of the pivot point in order to determine whether the algorithm is still efficient. For a generic move based on a lattice symmetry transformation $`g`$ applied in the pivot point $`\omega _k`$ we can define a local acceptance fraction $`f_{g,piv}(N,k)`$. In the scaling limit $`N\mathrm{}`$, $`k\mathrm{}`$ with $`k/N`$ finite, the local acceptance fraction is expected to satisfy the scaling Ansatz

$$f_{g,piv}(N,k)N^{p_{g,piv}}h_{g,piv}\left(k/N\right),$$
(3)
where $`p_{g,piv}`$ is the exponent which asymptotically governs the decay of the global acceptance fraction for $`N\mathrm{}`$ and $`h_{g,piv}\left(k/N\right)`$ is a scaling function. For small values of $`\alpha k/N`$ the behaviour of the scaling function $`h_{g,piv}(\alpha )`$ depends drastically on whether or not the transformation $`g`$ preserves the $`z`$ coordinate of the walk.

In order of analyse quantitatively the dynamic behaviour of the algorithm, we group the pivot moves into equivalence classes. Two pivot moves which are based on the symmetries $`g`$ and $`g^{}`$ are said to be equivalent if there exists a global symmetry $`\sigma `$—i.e. a transformation of the whole walk—which preserves the geometry of the system (i.e. a lattice symmetry which does not modify the $`z`$ coordinates) such that, for every walk $`\omega \mathrm{\Omega }^S(\omega _0)`$ and for every pivot point $`\omega _k`$, the walk which is obtained by applying a pivot move based on $`g^{}`$ to the globally $`\sigma `$-transformed walk $`\sigma (\omega )`$ in the pivot point $`\omega _k^{}=\sigma (\omega _k)`$ and then transforming back via $`\sigma ^1`$ coincides with the walk obtained from $`\omega `$ by means of a pivot move based on the symmetry $`g`$ with pivot point $`\omega _k`$. In formulae, if we denote by $`G_k`$ the operator associated with the pivot move based on $`g`$ applied on monomer $`k`$, then $`g`$ and $`g^{}`$ are equivalent, if, for every $`k`$, $`G_k=\sigma ^1G_k^{}\sigma `$. For instance, consider the transformations $`g:(x,y,z)(x,z,y)`$ and $`g^{}:(x,y,z)(x,z,y)`$. It is a simple matter to show that they are equivalent: It is enough to consider the $`\pi /2`$ rotation $`\sigma :(x,y,z)(x,z,y)`$.

In $`3`$ dimensions the original $`47`$ different moves are classified in $`15`$ different classes. Moves belonging to a given class have exactly the same critical behaviour. Therefore, in the following we will study Eq. (3) for different classes of moves.

In Fig. 1 we report the acceptance fraction for the 15 classes as a function of $`\alpha `$, for $`N=100`$. For small $`\alpha `$, that is for pivot points that are near the surface, there are essentially three types of behaviour. Moves that involve the inversion of the $`z`$-axis have a very low acceptance rate: for $`\alpha =0.1`$ the acceptance of, say, reflections with respect to the $`z`$-axis is $`0.04`$ and it drops further to $`0.0045`$ for $`\alpha =0.05`$. Transformations that involve $`\pi /2`$ rotations and diagonal reflections have a better behaviour although the acceptance still drops as $`\alpha 0`$. Transformations that do not modify the $`z`$-direction are instead unaffected by the presence of the surface.

Although the behaviour of the acceptance fraction as a function of $`\alpha `$ is radically different between SAW in the bulk and SAW in a half-space, the critical behaviour of the global acceptance fraction averaged over $`k`$ is very similar.

Denoting with $`f_{i,piv}(N,k)`$ the acceptance fraction of the class of moves $`i`$ applied in the pivot point $`\omega _k`$, and with $`f_N^{i,piv}`$ its average over $`k`$, we expect $`f_N^{i,piv}`$ to vanish as the length of the walk increases as

$$f_N^{i,piv}=\frac{1}{N}\underset{k=0}{\overset{N1}{}}f_{i,piv}(N,k)N^{p_{i,piv}}\left(A_i+B_i/N^{\mathrm{\Delta }_i}+\mathrm{}\right),$$
(4)
where the index $`i`$ indicates a symmetry class, $`\mathrm{\Delta }_i`$ is the corresponding leading correction-to-scaling exponent, and $`A_i`$, $`B_i`$ are non-universal amplitudes.

Before giving a quantitative estimate of the exponents $`p_{i,piv}`$, we try to give heuristically a rough estimate of the relative order of magnitude of the acceptance exponents in the bulk and in the presence of the surface. The argument which we use was already introduced in Ref.  for SAW’s in the bulk and we simply extend it to the case of SAW’s in a half-space.

If the two parts of the walk $`\omega [0,k]=\{\omega _0,\mathrm{},\omega _k\}`$ and $`\omega [k,N]=\{\omega _k,\mathrm{},\omega _N\}`$ are considered as independent, the local acceptance fraction for any applied lattice symmetry in the bulk would be given by

$$f_{i,piv}(N,k)\frac{c_N}{c_kc_{Nk}},$$
(5)
where $`c_k`$ is the number of walks of length $`k`$. For $`k\mathrm{}`$,

$$c_k\mu ^kk^{\gamma 1},$$
(6)
where $`\mu `$ is the critical fugacity and $`\gamma 1.1575(5)`$ , so that we have

$$f_N^{i,piv}N^{1\gamma },$$
(7)
and $`p_{i,piv}\gamma 10.16`$.

In the presence of the surface the same argument gives

$$f_{i,piv}(N,k)\frac{\stackrel{~}{c}_N(0)}{\stackrel{~}{c}_k(0)\stackrel{~}{c}_{Nk}(\omega _k)},$$
(8)
where $`\stackrel{~}{c}_k(\stackrel{}{r})`$ is the number of walks of length $`N`$ starting from $`\stackrel{}{r}`$. For $`k\mathrm{}`$, $`\stackrel{~}{c}_k(0)\mu ^kk^{\stackrel{~}{\gamma }1}`$, where $`\mu `$ is a fugacity which is identical to that appearing in Eq. (6) and $`\stackrel{~}{\gamma }`$ a new critical exponent. In three dimensions, $`\stackrel{~}{\gamma }0.68`$ . Assuming $`\omega _k`$ to be at a (macroscopically) finite distance from the surface we can take $`\stackrel{~}{c}_{Nk}(\omega _k)c_{Nk}`$, where $`c_{Nk}`$ scales according to Eq. (6). Then, Eq. (8) gives again $`f^{i,piv}N^{1\gamma }`$, i.e. $`p_{i,piv}\gamma 10.16`$. Thus, heuristically, we expect the acceptance exponents $`p_{i,piv}`$ to be of the same order in the presence of the surface and in the bulk. This is confirmed by the numerical estimates.

In order to compute the acceptance exponents $`p_i`$ we considered walks of length $`N=100`$, $`200`$, $`400`$, $`800`$, $`4000`$, $`8000`$, $`16000`$, $`32000`$ and performed fits of the global acceptance fraction for each class of moves to a simple power law. In order to study the effect of the subleading terms appearing in (4) we performed different fits using only data for $`NN_{\mathrm{min}}`$ for increasing values of $`N_{\mathrm{min}}=100,\mathrm{},16000`$. We obtained as a result the effective exponents which are reported in table 1 and which approach $`p_i`$ as the lower cutoff $`N_{\mathrm{min}}`$ increases.

Except for a few cases, which in the following we indicate with a star, the results are stable for walks of length $`N4000`$. Our best estimates for the different classes of moves are:

1. (a) $`z`$-axis inversion: $`p_{piv}=0.0919\left(47\right)^{}`$,
(b) $`x`$ or $`y`$-axis inversion: $`p_{piv}=0.0925\left(11\right)`$;
2. (a) $`\pm \pi /2`$ rotation in $`yz`$ \[or $`zx`$\] planes: $`p_{piv}=0.10202\left(81\right)`$,
(b) $`\pm \pi /2`$ rotation in the $`xy`$ plane: $`p_{piv}=0.10002\left(60\right)`$;
3. (a) $`\pi `$ rotation in $`yz`$ \[or $`zx`$\] planes: $`p_{piv}=0.1338\left(12\right)`$,
(b) $`\pi `$ rotation in the $`xy`$ plane: $`p_{piv}=0.1329\left(17\right)`$;
4. (a) diagonal reflection in the $`yz`$ \[or $`zx`$\] planes: $`p_{piv}=0.08875\left(94\right)^{}`$,
(b) diagonal reflection in the $`xy`$ plane: $`p_{piv}=0.0887\left(14\right)^{}`$;
5. (a) diagonal reflection in the $`yz`$ \[or $`zx`$\] plane and $`x`$ \[resp. $`y`$\] axis refl.: $`p_{piv}=0.13024\left(83\right)`$,
(b) diagonal reflection in the $`xy`$ plane and $`z`$-axis inversion : $`p_{piv}=0.1327\left(23\right)`$;
6. (a) $`\pm \pi /2`$ rotation in the $`yz`$ \[or $`zx`$\] planes and $`x`$ \[resp. $`y`$\] axis refl.: $`p_{piv}=0.13556\left(76\right)`$,
(b) $`\pm \pi /2`$ rotation in the $`xy`$ plane and $`z`$-axis inversion: $`p_{piv}=0.13371\left(69\right)`$;
7. 3-axes reflection: $`p_{piv}=0.1637\left(17\right)`$;
8. diagonal reflection in the $`yz`$ or $`zx`$ planes and $`\pm \pi /2`$ rotation in the $`xy`$ plane: $`p_{piv}=0.11811\left(77\right)`$;
9. two diagonal reflections, one in the $`yz`$ or in the $`zx`$ plane and the other in the $`xy`$ plane: $`p_{piv}=0.11813\left(94\right)`$.

Here we have grouped the equivalence classes in the presence of the surface in 9 groups. These groups correspond to the 9 equivalence classes for the algorithm in the bulk, where one can consider global transformations $`\sigma `$ that do not preserve the $`z`$ coordinate. From the results we observe that the exponent $`p_{piv}`$ depends only on the equivalence class in the bulk. In agreement with the above-reported heuristic argument, the presence of the surface has no influence on the acceptance exponents of the different classes of moves, but only on the shape of the scaling functions $`h_{i,piv}(\alpha )`$ introduced in Eq. (3). On the other hand, one should observe that the acceptance exponents which correspond to different bulk equivalence classes are different, with the only exception of classes 8. and 9., which can be hardly distinguished even in the presence of the surface. One should also notice that the exponent is larger for $`3`$-axes reflections than for diagonal reflections. As it has been already noticed in Ref.  this can be understood intuitively. Indeed, two subwalks on the opposite sides of a pivot point tend to be directed and occupy on average opposite regions of space. We can imagine, for instance, that the two subwalks occupy two opposite octants which touch at the pivot point. The subwalk $`\omega [k,N]`$ is not moved in the opposite octant by pivot moves based, for instance, on a one-axis reflection, a diagonal reflection or a $`\pi /2`$ rotation, while a 3-axes reflection will move the whole subwalk in the opposite octant. Therefore, it should have a higher probability of rejection.

The stability of most of the fits with $`N_{\mathrm{min}}`$ seems to rule out the possibility that all equivalence classes have the same exponent and that the discrepancies are due to residual corrections to scaling. Different symmetries have apparently different acceptance exponents. If this is correct, the global acceptance exponent $`p_{piv}`$ averaged over all transformations would coincide with the exponent of the pivot move which, in the limit $`N\mathrm{}`$, has the highest probability of being accepted, that is $`p_{piv}0.089`$. This exponent is lower than that reported in Ref.  which was obtained by averaging over all transformations. We mention that if we also perform the group average we obtain $`p_{piv}0.116(1)`$. The group average exponent is only slightly higher than what can be found using the data reported in Ref.  for the pivot algorithm in the bulk, for which $`p0.113`$111 We have been informed that the acceptance fractions $`f`$ reported in Ref.  for $`d=3`$ (but not for $`d=2`$) are in error because the program performed with probability $`1/48`$ an identity move (which is always accepted). Therefore, data there should be corrected by the map $`f(48/47)(f1/48)`$. This mistake was guessed by Tom Kennedy on the basis of his own pivot simulations (see Ref. ) and was confirmed by Madras and Sokal’s examination of their program. Aware of this problem, we recomputed the acceptance fraction exponent and found $`p0.113`$..

Now we can use our best estimates of the acceptance exponents for testing our scaling Ansatz (3) and determining the behaviour of the scaling functions $`h_{i,piv}(\alpha )`$ for $`\alpha k/N0`$. We will show that for $`\alpha 0`$

$$h_{i,piv}(\alpha )\alpha ^{q_i},$$
(9)
with $`q_i<0`$ for the “good” moves which do not change the $`z`$ coordinate of the walk, and $`q_i>0`$ otherwise.

In Fig. 2 we report $`f_{i,piv}(N,k)(N/100)^{p_{i,piv}}=\widehat{h}_{i,piv}(\alpha )`$ for 4 different classes of moves: 1.(b), 3.(b), 9. and 3.(a). They are representative of the different types of behaviour for $`\alpha 0`$. One observes very good scaling, the data for different $`N`$ falling one on top of each other, except at small values of $`\alpha `$ for those moves for which $`\widehat{h}_{i,piv}`$ is nonzero at small $`\alpha `$ in the range of $`N`$ values we have considered. Indeed, for moves 1.(b) and 3.(b) the estimates increase as $`N`$ increases: for $`N=100`$ we would estimate $`\widehat{h}_{1.(\mathrm{b}),piv}(0)0.92`$ and $`\widehat{h}_{3.(\mathrm{b}),piv}(0)0.84`$, while for $`N=32000`$ we have $`\widehat{h}_{1.(\mathrm{b}),piv}(0)1.12`$ and $`\widehat{h}_{3.(\mathrm{b}),piv}(0)1.10`$. For the moves of class 9., instead, the curve for $`N=100`$ is slightly higher than that for $`N=4000`$, $`32000`$: for $`N=100`$ we would estimate $`\widehat{h}_{9.,piv}(0)0.15`$, while for $`N=4000`$, $`32000`$ we would obtain $`\widehat{h}_{9.,piv}(0)0.07`$, although the data for $`N=32000`$ seem to be even lower. Clearly, for small $`\alpha `$, there are significant corrections to scaling and indeed we will now show that $`h_{1.(\mathrm{b}),piv}(\alpha )`$ and $`h_{3.(\mathrm{b}),piv}(\alpha )`$ diverge for $`\alpha 0`$, while $`h_{9.,piv}(0)=0`$.

Consider first moves in classes 1.(b) and 3.(b) and in general all the moves preserving the $`z`$-coordinate of the walks. They are only marginally affected by the presence of the surface. The only effect is that, since the surface induces a monomer depletion near the surface, a move applied in a pivot point $`\omega _k`$ with small $`k`$ has higher probability of being accepted than the same move applied in $`\omega _{Nk}`$, explaining the slight asymmetry of the scaling curves $`h_{i,piv}(\alpha )`$. Moreover, we expect the local acceptance fraction $`f_{i,piv}(N,k)`$ to remain finite as $`N`$ increases with $`k`$ fixed and small. By comparing with the scaling Ansatz (3), it follows that $`h_{i,piv}(\alpha )\alpha ^{p_{i,piv}}`$, i.e. $`q_i=p_{i,piv}`$ and the scaling function diverges for $`\alpha 0`$. This is confirmed by the numerical results: indeed, $`\widehat{h}_{1.(\mathrm{b}),piv}(\alpha )`$ and $`\widehat{h}_{3.(\mathrm{b}),piv}(\alpha )`$ increase as $`N\mathrm{}`$.

Let us now consider the classes of moves which change, but do not invert, the $`z`$-coordinate, for instance class 9. In order to understand the behaviour of $`h_{i,piv}(\alpha )`$ for $`\alpha 0`$ we must compute the local acceptance probability for $`k`$ fixed and small. For this purpose, we have considered the acceptance fraction averaged over the pivot points $`kk_{\mathrm{max}}`$, i.e.

$$\stackrel{~}{f}_{9.,\mathrm{piv}}(N,k_{\mathrm{max}})=\frac{1}{k_{\mathrm{max}}+1}\underset{k=0}{\overset{k_{\mathrm{max}}}{}}f_{9.,piv}(N,k).$$
(10)
For fixed $`k_{\mathrm{max}}`$, it decays faster than the global acceptance fraction. For instance, for $`k_{\mathrm{max}}=20`$, for the moves that involve $`\pi /2`$ rotations or diagonal reflections in the $`(y,z)`$, $`(z,x)`$ planes (classes 8. and 9.), we have numerically found $`\stackrel{~}{f}_{8.,9.,,piv}(N,20)N^{0.49}`$. It will be shown in the following that this exponent is close to the one which characterizes the probability that a generic bulk $`N`$-step SAW grafted at the surface belongs to $`\mathrm{\Omega }_N^S(\omega _0)`$. Since we expect that $`f_{8.,9.,piv}(N,k_{\mathrm{max}})\stackrel{~}{f}_{8.,9.,piv}(N,k_{\mathrm{max}})`$ for small $`k`$, the scaling function $`h_{8.,9.,piv}(\alpha )`$ is expected to vanish with $`q_{8.,9.}0.49p_{8.,9.,piv}0.37`$. Similar behaviour is expected for all classes of moves ( 2.(a), 4.(a), 5.(a), 6.(a), 8. and 9.) which do not invert the $`\widehat{z}`$ axis.

The behavior of moves like 1.(a) that invert the $`z`$-axis is much worse. For instance, we performed a long run for $`N=800`$, in which such moves were attempted on the first 20 monomers approximately 21000 times, and none of the attempts was successful. Fitting directly $`\widehat{h}_{1.(\mathrm{a})}(\alpha )`$ for $`\alpha =k/N0.1`$ for $`N=100,800,4000`$, one finds $`q_{1.(\mathrm{a})}=3.8(3)`$.

It can be useful to get heuristically a rough idea of the order of magnitude of the exponent $`q_{1.(\mathrm{a})}`$. Let us denote with $`m_N(\omega ,k)`$ the maximum elongation in the $`z`$ direction of the monomers following the pivot point

$$m_N(\omega ,k)\underset{i=k+1,\mathrm{},N}{\mathrm{max}}\omega _i^z,$$
(11)
where $`\omega _i^z`$ is the $`z`$-coordinate of the monomer $`\omega _i`$. If we denote with $`\omega _k^z`$ the $`z`$-coordinate of the pivot point, the proposed move does not give rise to collisions with the surface if the inequality

$$\omega _k^z>\frac{m_N(\omega ,k)}{2}$$
(12)
is satisfied. Thus, if the pivot point is one of the first vertices following the grafted end, Eq. (12) states that the move may be successful only for those walks that have nearly all their monomers within few lattice planes from the surface and therefore we expect that $`h(\alpha )`$ vanishes for $`\alpha 0`$. It is known (see Ref. ) that for $`N\mathrm{}`$, $`z\mathrm{}`$ and $`\zeta z/N^\nu `$ fixed and small, the probability distribution $`P_N(z)`$ of the fraction of monomers which lay at a distance $`z`$ from the surface scales as

$$P_N(z)\frac{1}{N^\nu }\zeta ^{\frac{1}{\nu }1}.$$
(13)
Therefore, the fraction $`W(z)`$ of monomers that lay within a distance $`z`$ from the surface scales as

$$W(z)=_0^zP_N(z^{})𝑑z^{}\zeta ^{1/\nu }.$$
(14)
If $`k`$ is small but already in the scaling regime, we expect that $`\omega _k^zk^\nu `$ and therefore $`W(2\omega _k^z)k/N=\alpha `$. Therefore, since the average fraction of monomers laying in the strip $`z<2\omega _k^z`$ is $`\alpha `$, the average fraction of walks such that all their monomers lay in the strip is less than $`\alpha `$. It follows that $`q_{1.(\mathrm{a})}1`$.

Clearly, the pivot move is inefficient when the pivot is near the surface. Therefore, we expect the dynamics of observables that strongly depend on the behavior of the walk near the surface—for instance, the number of monomers on the surface—to be much slower than that of global observables. Also, note that $`\tau _{exp}`$ becomes rapidly large, making it very difficult to thermalize the system using only the pivot algorithm. In the following Section we discuss a second nonlocal move that solves the problems we have discussed.

## 3 The Cut-and-permute move

In this Section we want to consider a different nonlocal move that is able to modify the walk near the surface. By adding it to the pivot algorithm we will obtain an algorithm without the quasi-ergodicity problems we have discussed above. Of course, we do not want to destroy the optimal dynamic behaviour of the pivot algorithm, and thus we want to introduce a move for which the mean CPU time per succesful move scales simply with $`N`$. The move we introduce here consists in cutting the walk in two parts, and in rebuilding the walk in reverse order.

The elementary cut-and-permute move works as follows:

1. Given an $`N`$-step SAW $`\omega \mathrm{\Omega }_N^S(\omega _0)`$ , choose with probability $`p_c`$ a cut-point $`\omega _c`$ with $`c\{1,2,\mathrm{},N1\}`$. The point $`\omega _c`$ divides the walk in two subwalks $`\omega ^1=\{\omega _1,\mathrm{},\omega _c\}`$ and $`\omega ^2=\{\omega _c,\mathrm{},\omega _N\}`$.
2. Cut the walk in the cut-point and rebuild it arranging the two subwalks $`\omega ^1`$ and $`\omega ^2`$ in reverse order. The resulting walk $`\omega ^{}`$ has vertices at positions

$$\omega _{}^{}{}_{i}{}^{}=\{\begin{array}{cc}\omega _0+\omega _{c+i}\omega _c\hfill & \text{for }0iNc\text{,}\hfill \\ \omega _N\omega _c+\omega _{iN+c}\hfill & \text{for }Nc+1iN\text{.}\hfill \end{array}$$
(15)
3. The proposed move is accepted if the resulting walk $`\omega ^{}`$ is self-avoiding and does not intersect the excluded region.

It is easy to see that the move satisfies detailed balance as long as $`p_c=p_{Nc}`$.

It is also possible to incorporate in the cut-and-permute move the “good” pivot transformations, i.e. those that do not change the $`z`$-coordinate. We thus define an “improved” cut-and-permute move as follows (the move is illustrated in Fig. 3):

1. Given an $`N`$-step SAW $`\omega \mathrm{\Omega }_N^S`$ , choose with probability $`p_c`$ a cut-point $`\omega _c`$ with $`c\{1,2,\mathrm{},N1\}`$. The point $`\omega _c`$ divides the walk in two subwalks $`\omega ^1=\{\omega _1,\mathrm{},\omega _c\}`$ and $`\omega ^2=\{\omega _{c+1},\mathrm{},\omega _N\}`$.
2. Choose with probability $`P_g`$ an element $`g`$ belonging to the symmetry group of the plane which is parallel to the surface. In three dimensions the symmetry group is the dihedral group $`D_4`$ in the $`z=0`$ plane, whose 8 elements are $`\pm \pi /2`$ rotations, $`\pi `$ rotation, axis inversions, diagonal reflections, and the identity.
3. Apply the chosen symmetry to the subwalk $`\omega ^2`$ to obtain the subwalk $`\omega _{}^{2}{}_{}{}^{}=\{\omega _c,\omega _c+g(\omega _{c+1}\omega _c),\mathrm{},\omega _c+g(\omega _N\omega _c)\}`$.
4. Cut the walk in the cut-point and rebuild it arranging the two subwalks $`\omega ^1`$ and $`\omega _{}^{2}{}_{}{}^{}`$ in reverse order. The resulting walk $`\omega ^{}`$ has vertices at positions

$$\omega _{}^{}{}_{i}{}^{}=\{\begin{array}{cc}\omega _0+g(\omega _{c+i}\omega _c)\hfill & \text{for }0iNc\text{,}\hfill \\ g(\omega _N\omega _c)+\omega _{iN+c}\hfill & \text{for }Nc+1iN\text{.}\hfill \end{array}$$
(16)
5. The proposed move is accepted if the resulting walk $`\omega ^{}`$ is self-avoiding and does not intersect the excluded region.

The improved cut-and-permute move does not satisfy detailed balance. Indeed, the walk $`\omega ^{}`$ obtained from $`\omega `$ with a cut-and-permute move in the cut-point $`\omega _c`$ cannot be transformed back to the original walk $`\omega `$ by any move applied in $`\omega _{Nc}`$, since the subwalk of $`\omega `$ following the cut-point has changed its orientation in space and its orientation cannot be modified in the second cut-and-permute move. However, the move still leaves the probability measure invariant as long as $`p_c=p_{Nc}`$. Indeed, if we perform a move based on $`g`$ in the point $`\omega _c`$ and then a move based again on $`g`$ in the point $`\omega _{Nc}`$, we obtain the walk $`\omega ^{\prime \prime }=g\left(\omega \right)`$ corresponding to the application of symmetry $`g`$ to the original walk. If we denote by $`p(\omega \omega ^{})`$ the probability that the move applied to $`\omega `$ gives as result $`\omega ^{}`$, from what we have described it follows that $`p(\omega \omega ^{})=p(\omega ^{}g(\omega ))`$. Thus,

$$\underset{\omega ^{}}{}\pi \left(\omega ^{}\right)p\left(\omega ^{}\omega \right)=\underset{\omega ^{}}{}\frac{1}{c_N^S}p\left(g^1(\omega )\omega ^{}\right)=\frac{1}{c_N^S}=\pi \left(\omega \right),$$
(17)
as required. In the following we will consider the improved cut-and-permute move, but, as we shall show explicitly, nothing would change by using the simpler version.

As it has been done for the pivot algorithm, we can define equivalence classes of cut-and-permute moves and study the acceptance fraction $`f_{i,cp}(N,c)`$ as a function of the position $`c`$ of the cut-point and its average $`f_N^{i,cp}`$ over all points $`c`$. For large $`N`$, we expect a scaling behavior of the form

$$f_N^{i,cp}N^{p_{cp}}.$$
(18)
We have not added an index $`i`$ to $`p_{cp}`$ because, as we shall see, this exponent does not depend on the equivalence class.

To give an estimate of the exponent $`p_{cp}`$, we can use again a heuristic argument in which we consider the two subwalks in which the original walk $`\omega `$ is divided at the cut-point $`\omega _c`$ as independent. We also assume that the probability that the attempted move is accepted is the product of the probability $`P_1(c)`$ that the concatenation of the two subwalks gives a walk which does not have self-collisions times the probability $`P_2(c)`$ that the subwalk $`\omega _2`$, transformed under a $`D_4`$ symmetry and translated with its first vertex on the surface does not intersect the forbidden region. Repeating the argument presented for the pivot case, we have

$$P_1(c)=\frac{\stackrel{~}{c}_N}{\stackrel{~}{c}_cc_{Nc}},$$
(19)
which, by averaging over $`c`$, gives $`P_1N^{(\gamma 1)}N^{0.1575(5)}`$ . The second probability is given by

$$P_2(c)=\frac{\stackrel{~}{c}_{Nc}}{c_{Nc}},$$
(20)
which, after averaging, gives $`P_2N^{(\gamma \stackrel{~}{\gamma })}N^{0.48}`$. The global acceptance fraction of the cut-and-permute move is expected to scale as $`f_{cp}=P_1P_2N^{p_{cp}}N^{2\gamma +1+\stackrel{~}{\gamma }}N^{0.63}`$. This argument is likely to give an overestimate of the exponent $`p_{cp}`$ for the following reasons. The original subwalk $`\omega _1`$ which is grafted at the surface is directed, so that $`P_1(c)`$ is underestimated. Moreover, the walk $`\omega _2`$ is expected to have a residual directionality in the $`+\widehat{z}`$ direction so that the probability that $`\omega _2`$ does not intersect the forbidden region is expected to be larger than $`P_2(c)`$. Table 2 contains numerical estimates of the acceptances $`f_N^{cp,i}`$ for different classes of moves in $`D_4`$ at different values of $`N`$, averaged over the cutting point $`c`$. We have fitted the data for $`NN_{\mathrm{min}}`$ to a power law considering increasing values of $`N_{\mathrm{min}}`$. Systematic errors due to corrections to scaling should become negligible for $`N_{\mathrm{min}}\mathrm{}`$. The $`\chi ^2`$ values of fits with $`N_{\mathrm{min}}4000`$ indicate that the fits are stable, with the exception of two cases. These two cases are denoted with a star in the following list containing our best estimates for the exponent $`p_{cp}`$:

1. identity $`p_{cp}=0.4969\left(35\right)`$;
2. diagonal reflections $`p_{cp}=0.4903\left(29\right)`$;
3. $`\pm \pi /2`$ rotations $`p_{cp}=0.4911\left(21\right)`$;
4. $`\pi `$ rotation $`p_{cp}=0.4873\left(40\right)^{}`$;
5. 1-axis reflection $`p_{cp}=0.4929\left(13\right)^{}`$.

At variance with the pivot case the results for the different equivalence classes are compatible within error bars: the transformation $`g`$ seems to play little role. If we average over all symmetries we obtain $`p_{cp}=0.4922(20)`$. Note that, as expected, this result is somewhat lower than the heuristic estimate reported above.

The exponent $`p_{cp}`$ is significantly higher than $`p_{piv}`$ so that cut-and-permute moves will be more rarely accepted. Numerically, we find that a cut-and-permute move is accepted every $`2`$ (resp. $`20`$) succesful pivot moves for $`N=100`$ (resp. $`N=32000`$). However, they represent the moves that most effectively change the conformation of the walk near the grafted end and thus they will play an important role in decorrelating the walk configurations.

As for the pivot case, it is interesting to study the local acceptance fraction as a function of $`c`$ and $`N`$. For $`N\mathrm{}`$, $`\zeta c/N`$ fixed we expect a scaling form

$$f_{cp}(N,c)N^{p_{cp}}h_{cp}(\zeta ),$$
(21)
where $`h_{cp}(\zeta )`$ is a scaling function which encodes the dependence on the different cut-points. In Figure 4 we report the function $`\widehat{h}_{cp}(\zeta )=f_{cp}(N,c)(N/100)^{p_{cp}}`$ for 5 values of $`N`$ using the exponent $`p_{cp}`$ given by our fits. There is no difference in the local behaviour of the 5 classes of equivalent moves, so the function we report refers to the total acceptance obtained by averaging over the equivalence classes.

Note the very good scaling: the points fall on top of each other for all values of $`N`$.

Let us analyze now the scaling behaviour of the acceptance at fixed $`c`$. We define the observable

$$\stackrel{~}{f}_{cp}(N,c_{\mathrm{max}})=\frac{1}{2c_{\mathrm{max}}}\left[\underset{c=1}{\overset{c_{\mathrm{max}}}{}}f_{cp}(N,c)+\underset{c=Nc_{\mathrm{max}}}{\overset{N1}{}}f_{cp}(N,c)\right],$$
(22)
and study its scaling with $`N`$ as $`N\mathrm{}`$. For $`c_{\mathrm{max}}=20`$ a power-law fit gives $`\stackrel{~}{f}_{cp}(N,20)N^{0.01}`$. This is just an effective exponent and indeed it decreases in modulus as the minimum value of $`N`$ considered in the fit is increased. The results are therefore compatible with $`\stackrel{~}{f}_{cp}(N,20)\mathrm{const}`$ for $`N_{\mathrm{min}}4000`$. Such a result is easy to understand. Indeed, suppose first that $`Nc_{\mathrm{max}}c<N`$. Then the move is accepted if $`\omega ^{}[0,c^{}]`$, $`c^{}=Nc`$, does not intersect the surface—which happens with a probability independent of $`N`$—and does not intersect $`\omega ^{}[c^{},N]`$, which, for large $`N`$, should be roughly independent of $`N`$. Thus, $`f_{cp}(N,c)`$ should be constant for $`N`$ large. If $`1cc_{\mathrm{max}}`$, note that $`f_{cp}(N,c)=f_{cp}(N,Nc)`$, to conclude the argument. In conclusion, the cut-and-permute move has a good scaling behavior exactly in the region in which the pivot moves behave badly. The addition of the cut-and-permute move in the hybrid dynamics is therefore crucial in speeding up the dynamics of the first steps of the walk.

As in the pivot case, since the local acceptance fraction converges to a constant for small or large fixed values of $`\zeta `$, the scaling function $`h_{cp}(\zeta )`$ increases as $`\zeta ^{p_{cp}}`$ for $`\zeta 0`$ (resp. $`(1\zeta )^{p_{cp}}`$ for $`\zeta 1`$.) This suggests that in step 1. one can take $`p_c`$ different from 0 only for $`\zeta `$ near 0 and 1, i.e. set $`p_c=0`$ if $`\alpha NcN(1\alpha )`$. For instance, by taking $`\alpha =1/10`$, we would approximately increase the acceptance by a factor $`1.7`$.

## 4 Computational complexity

In order to study the effective dynamic behaviour of the algorithm in CPU time, it is important to determine the scaling behaviour of the CPU time needed to generate pivot and cut-and-permute successful moves.

The CPU time depends on the data structure that is used. For the walk we consider a sequentially allocated linear list and in order to check for self-collisions and collisions with the surface we use a hash table as described in Ref. . This allows to insert a single monomer in an average CPU time of order one. In the following we will define a CPU time unit as the CPU time which is needed to upgrade a single monomer position.

In order to optimize the algorithm, it is important to mimimize the CPU time spent in a failed attempt. This requires a careful choice of the order in which the new walk $`\omega ^{}`$ is rebuilt.

For the pivot algorithm, Ref.  suggested to build the new walk and perform the self-avoidance check in the following order: $`\omega _k^{}`$, $`\omega _{k+1}^{}`$, $`\omega _{k1}^{}`$, $`\mathrm{}`$,$`\omega _{k+i}`$, $`\omega _{ki}`$, $`\mathrm{}`$ till all monomers have been checked. The motivation was that the probability of self-intersections is higher for those monomers which are close to the pivot point. However, in the presence of a surface a move that changes the $`z`$-coordinate can also fail because the new walk intersects the surface.

Intersections with the surface occur with higher probability for those monomers of the subwalk $`\omega [k,N]`$ which are more likely to undergo a large displacement, while monomers belonging to the grafted subwalk $`\omega [0,k]`$ are not modified and satisfy the geometrical constraint automatically. In order to keep into account this effect, for the moves that change the $`z`$-coordinates, we can use a different strategy: we insert the monomers of the proposed walk $`\omega ^{}`$ in the hash table in the following order: $`\omega _k^{}`$, $`\omega _N^{}`$, $`\omega _{k+1}^{}`$, $`\omega _{k1}^{}`$, $`\omega _{N1}^{}`$, $`\mathrm{}`$, $`\omega _{k+i}^{}`$, $`\omega _{ki}^{}`$, $`\omega _{Ni}`$, $`\mathrm{}`$ , till every monomer has been checked once. As a measure of the CPU time spent by the algorithm we can use the number of walk monomers $`I(\omega ^{})`$ that are checked in the move. If the move is successful, clearly $`I(\omega ^{})=N`$, while for a failed attempt

$`I(\omega ^{})`$ $`=`$ $`\mathrm{min}\{i:\omega ^{}[ki,k1]\omega ^{}[k,k+i]\mathrm{}`$ (23)
orω[Ni,N]

S
orω[ki,k+i]

S
},\displaystyle{\rm\,or\,}\omega^{\prime}\left[N-i,N\right]\,\cap\,\vbox{\hbox{\hbox{}\kern 2.5pt$\scriptstyle\circ$}\vskip-11.62494pt\hbox{$S$}}\neq\emptyset{\rm\,or\,}\omega^{\prime}\left[k-i,k+i\right]\,\cap\,\vbox{\hbox{\hbox{}\kern 2.5pt$\scriptstyle\circ$}\vskip-11.62494pt\hbox{$S$}}\neq\emptyset\},
where, following Ref. , we denote the subwalk $`\{\omega _{\mathrm{max}(i,0)}^{},\omega _{\mathrm{max}(i,0)+1}^{},\mathrm{}\omega _{\mathrm{min}(j,N)}^{}\}`$, with the symbol $`\omega ^{}[i,j]`$ and the half-space $`z<0`$ with  $``$ $`S`$ . We want now to evaluate the average time spent in checking a failed move that is expected to scale as

$$T_N^{piv,failed}I(\omega ^{})|_{\mathrm{failed}}N^{y_{piv}}.$$
(24)
where the average over all failed moves is taken. An estimate of $`y_{piv}`$ can be obtained by a heuristic argument. Inserting monomers in the hash table in the way we described before, the probability that $`I(\omega ^{})`$ exceeds the value $`i`$ can be estimated from the scaling of the acceptance fraction. Indeed, we expect

$`\mathrm{Prob}\{I(\omega ^{})>i\}`$ $`=`$ $`\mathrm{Prob}\{\omega ^{}[ki,k+i]\mathrm{is}\mathrm{SAW}\mathrm{and}\mathrm{does}\mathrm{not}\mathrm{intersect}\mathrm{the}\mathrm{surface},`$ (25)
$`\mathrm{and}\omega ^{}[Ni,N]\mathrm{does}\mathrm{not}\mathrm{intersect}\mathrm{the}\mathrm{surface}\}i^{p_{piv}}.`$
Therefore, the average time spent in checking a failed move is expected to scale like

$$T_N^{piv,failed}=\underset{i=0}{\overset{N}{}}\mathrm{Prob}\{I>i\}N^{1p_{piv}},$$
(26)
so that $`y_{piv}=1p_{piv}`$. Estimates of $`T_N^{piv,failed}`$ are reported in Table 3, for the different classes of pivot moves and different values of $`N`$. We performed fits of the form (24) in order to determine $`y_{piv}`$, including each time only data with $`NN_{\mathrm{min}}`$ (in Table 3, $`N_{\mathrm{min}}`$ is the value $`N`$ reported in the first column). We observed that in all cases $`y_{piv}`$ increases with $`N_{\mathrm{min}}`$ and thus our results are probably a lower estimate of $`y_{piv}`$. Here, we report the results with $`N_{\mathrm{min}}=4000`$ and indicate with a star those cases in which the $`\chi ^2`$ is still too large:

1. $`z`$-axis inversion: $`y_{piv}=0.8750\left(16\right)`$,
$`x`$ or $`y`$-axis inversion: $`y_{piv}=0.8234\left(17\right)^{}`$;
2. $`\pm \pi /2`$ rotation in $`yz`$ \[or $`zx`$\] planes: $`y_{piv}=0.8554\left(13\right)^{}`$,
$`\pm \pi /2`$ rotation in the $`xy`$ plane: $`y_{piv}=0.8205\left(23\right)`$;
3. $`\pi `$ rotation in $`yz`$ \[or $`zx`$\] planes: $`y_{piv}=0.82756\left(97\right)^{}`$,
$`\pi `$ rotation in the $`xy`$ plane: $`y_{piv}=0.7989\left(46\right)`$;
4. diagonal reflection in the $`yz`$ \[or $`zx`$\] planes: $`y_{piv}=0.8635\left(22\right)`$,
diagonal reflection in the $`xy`$ plane: $`y_{piv}=0.8470\left(25\right)^{}`$;
5. diagonal reflection in the $`yz`$ \[or $`zx`$\] plane and $`x`$ \[resp. $`y`$\] axis refl.: $`y_{piv}=0.8381\left(13\right)^{}`$,
diagonal reflection in the $`xy`$ plane and $`z`$-axis inversion : $`y_{piv}=0.8342\left(16\right)`$;
6. $`\pm \pi /2`$ rotation in the $`yz`$ \[or $`zx`$\] planes and $`x`$ \[resp. $`y`$\] axis refl.: $`y_{piv}=0.8250\left(18\right)`$,
$`\pm \pi /2`$ rotation in the $`xy`$ plane and $`z`$-axis inversion: $`y_{piv}=0.8339\left(23\right)`$;
7. 3-axis reflection: $`y_{piv}=0.8138\left(35\right)`$;
8. diagonal reflection in the $`yz`$ or $`zx`$ planes and $`\pm \pi /2`$ rotation in the $`xy`$ plane: $`y_{piv}=0.8476\left(20\right)^{}`$;
9. two diagonal reflections, one in the $`yz`$ or in the $`zx`$ plane and the other in the $`xy`$ plane: $`y_{piv}=0.85121\left(67\right)^{}`$.

If we average over all equivalence classes we obtain $`y_{piv}=0.8428(13)`$. These results should not be trusted too much, and in many cases the correct estimate is probably higher. However, the important thing that emerges is that $`y_{piv}\mathrm{}<1p_{piv}`$, a result that will be used below.

We also computed $`T_N^{piv,failed}`$ for the moves that change the $`z`$-coordinate when we insert the points as proposed in Ref. , i.e. without modifying the order in which we insert monomers into the hash table because of the presence of the surface. We found that, by using this method, the average number of points inserted in the hash table before detecting a failure is approximately $`1.4`$ times larger. The difference is more significant for pivot points $`\omega _k`$ with $`k<N/2`$—for $`k/N0`$ the difference is a factor of two—while for $`k\mathrm{}>0.6`$ the two methods are equivalent. However, this has no influence on the exponent $`y_{piv}`$, which is approximately the same in both cases.

At this point we can estimate the average CPU time spent in a pivot move. Clearly

$`T_{piv}`$ $`=`$ $`(1f_N^{piv})T_N^{piv,failed}+f_N^{piv}T_N^{piv,succ}`$ (27)
$``$ $`1N^{y_{piv}}+N^{p_{piv}}NN^{1p_{piv}},`$
where we have taken into account that $`y_{piv}\mathrm{}<1p_{piv}`$. Notice also that more time is spent in a successful move than in an unsuccesful one.

Let us now consider the cut-and-permute move. As in the previous case we tried to build the new walk in the most efficient way in order to reduce the CPU time necessary to detect a failure. Unlike the pivot case, for the cut-and-permute move the heuristic argument indicates that the dominant cause of failure is given by intersections with the surface. As it is illustrated in Fig. 3, the cut-point $`\omega _c`$ divides the walk in two subwalks, $`\omega [0,c]`$, which is grafted at the surface, and $`\omega [c,N]`$. The move transforms $`\omega [c,N]`$ via a lattice symmetry $`g`$ which preserves $`z`$ and translates the transformed subwalk grafting its first vertex at the surface. The subwalk $`\omega [0,c]`$ is instead translated in order to join the free end of the grafted subwalk. The only possibility of hitting the surface comes from the new grafted subwalk $`\omega ^{}[0,Nc]`$ and the monomers which are more likely to intersect the forbidden region are the ones following the new grafted vertex. Collision between monomers are instead more likely to occur at the joining of the two subwalks.

As for the pivot move, in order to check at the same time the two regions where failures are more likely to occur, we inserted the monomers of the proposed new walk $`\omega ^{}`$ in the hash table in the following order: $`\omega _{Nc}^{}`$, $`\omega _0^{}`$, $`\omega _{Nc1}^{}`$, $`\omega _{Nc+1}^{}`$, $`\omega _1^{}`$, $`\mathrm{}`$, $`\omega _{i1}^{}`$, $`\omega _{Nci}^{}`$, $`\omega _{Nc+i}^{}`$, $`\mathrm{}`$.

In this way the CPU time needed to check a proposed walk is

$`I(\omega ^{})`$ $`=`$ $`\mathrm{min}\{i:\omega ^{}[Nci,Nc1]\omega ^{}[Nc+1,Nc+i]\mathrm{}`$ (28)
orω[0,i]

S
orω[Nci,Nc+i]

S
},\displaystyle{\rm\,or\,}\omega^{\prime}\left[0,i\right]\,\cap\,\vbox{\hbox{\hbox{}\kern 2.5pt$\scriptstyle\circ$}\vskip-11.62494pt\hbox{$S$}}\neq\emptyset{\rm\,or\,}\omega^{\prime}\left[N-c-i,N-c+i\right]\,\cap\,\vbox{\hbox{\hbox{}\kern 2.5pt$\scriptstyle\circ$}\vskip-11.62494pt\hbox{$S$}}\neq\emptyset\},
if $`\omega ^{}`$ is rejected and $`I(\omega ^{})=N`$ if it is accepted. The average time spent in detecting a failed move is expected to scale as

$$T_N^{cp,failed}=I(\omega ^{})|_{\mathrm{failed}}N^{y_{cp}}.$$
(29)
We can estimate $`y_{cp}`$ heuristically as before. With the above procedure the probability of using CPU time which is longer than $`i`$ goes as

$`\mathrm{Prob}\{I(\omega ^{})>i\}=`$
$`\mathrm{Prob}\{\omega ^{}[Nki,Nk+i]\mathrm{is}\mathrm{SAW}\mathrm{and}\omega ^{}[0,i]\mathrm{does}\mathrm{not}\mathrm{intersect}\mathrm{the}\mathrm{surface}\}`$
$`i^{p_{cp}},`$ (30)
and the average total time spent in a failure scales as

$$T_N^{cp,failed}=\underset{i=0}{\overset{N}{}}Prob\{I>i\}N^{1p_{cp}},$$
(31)
so that $`y_{cp}=1p_{cp}`$. The data for $`T_N^{cp,failed}`$ in the cut-and-permute move and the estimates of the $`y_{cp}`$ exponent obtained from data with $`NN_{min}`$, with $`N_{min}`$ appearing in the first column, are reported in Table 4. From those data we have computed our final estimates as we did in the pivot case obtaining:

1. diagonal reflections $`y_{cp}=0.4870\left(19\right)`$;
2. identity $`y_{cp}=0.4788\left(70\right)`$;
3. $`\pm \pi /2`$ rotations $`y_{cp}=0.4921\left(13\right)`$;
4. $`\pi `$ rotation $`y_{cp}=0.4933\left(43\right)^{}`$;
5. 1 axis reflection $`y_{cp}=0.4869\left(13\right)^{}`$.

The group average $`y_{cp}=0.4890(15)`$ is slightly lower than $`1p_{cp}0.5078(20)`$ as it is also happened for the pivot move.

In order to confirm the role played by intersections with the surface, we also measured $`I(\omega ^{})`$ when monomers are inserted in the hash table in the order $`\omega _{Nc}^{}`$, $`\omega _{Nc+1}^{}`$, $`\omega _{Nc1}^{}`$, $`\mathrm{}`$, without inserting at the same time the monomers near the grafted end. In this case we find $`y_{cp}0.9`$ which is very similar to the exponent for the pivot move. The proposed procedure is thus very inefficient.

Finally, we can estimate the average CPU time spent in a cut-and-permute move. Clearly

$`T_{N,cp}`$ $`=`$ $`(1f_N^{cp})T_N^{cp,failed}+f_N^{cp}T_N^{cp,succ}`$ (32)
$``$ $`1N^{y_{cp}}+N^{p_{cp}}NN^{1p_{cp}},`$
where we have taken into account $`y_{cp}\mathrm{}<1p_{cp}`$. Note also that $`T_{N,cp}/T_{N,piv}N^{p_{cp}p_{piv}}N^{0.4}`$, so that the average CPU time for successful move scales identically.

## 5 The hybrid cut-and-permute algorithm: autocorrelation times

In this Section we want finally to define our hybrid algorithm that consists in performing both pivot and cut-and-permute moves. More precisely, the algorithm is specified by a number $`0<q<1`$. The hybrid algorithm works by performing a pivot move with probability $`q`$ and a cut-and-permute move with probability $`(1q)`$. Since the pivot alone is ergodic, the full algorithm is ergodic.

In order to understand the dynamic behavior we will distinguish three different classes of observables: (a) global observables that depend on the global shape of the SAW, for instance, the end-to-end distance or the radius of gyration; (b) global surface observables that depend on the interaction of the SAW with the surface, for example, the number of monomers lying on $`S`$; (c) local observables that depend on local properties of the SAW.

Let us analyse at first the behaviour of global observables. Both pivot and cut-and-permute are nonlocal moves and can decorrelate global observables in few accepted moves, but the cut-and-permute move cannot modify a class of observables, for instance the $`z`$-coordinate of the free endpoint. We denote with $`O_z`$ the class of observables which are updated only in the pivot dynamics and with $`O_{x,y}`$ all the others. For a polymer of length $`N`$ the autocorrelation time $`\tau _{int,O_z}`$ is therefore of order $`1/(qf_{N,piv})`$, while $`\tau _{int,O_{x,y}}`$ is of order $`1/f_{N,max}`$, where $`f_{N,max}=\mathrm{max}\{qf_{N,piv},(1q)f_{N,cp}\}`$. Since in the limit $`N\mathrm{}`$ the acceptance fraction of the cut-and-permute move decreases faster than the acceptance fraction of the pivot move, if $`q`$ is constant, the scaling behaviour of the autocorrelation time $`\tau _{int,glob}`$ of any global observable is determined by the pivot dynamics only, i.e.

$$\tau _{int,glob}\tau _{int,O_z}\tau _{int,O_{x,y}}q^1N^{p_{piv}},$$
(33)
so that $`z_{glob}p_{piv}0.10`$. Table 5 reports the autocorrelation times for different values of $`N`$ for the $`z`$ component of the free endpoint obtained setting $`q=1/2`$. The dynamic exponent $`z`$, which is obtained by fitting the data to a power law, is $`z0.15`$, in reasonable agreement with the prediction.

Global surface observables will be changed in a few successful cut-and-permute moves and thus we expect

$$\tau _{int,surf}(1q)N^{p_{cp}},$$
(34)
so that $`z_{surf}p_{cp}0.49`$. Unfortunately, we have not measured surface quantities in our simulation and therefore we are not able to test this prediction. Nonetheless, we believe that the correct value of $`z`$ is not too far from $`p_{cp}`$.

Finally, let us consider a local observable, for instance the average number of monomers at which the SAW makes a 90o turn. Such a quantity is of order $`N`$ with a variance of order $`\sqrt{N}`$. By using a standard random-walk argument we obtain

$$\tau _{int,loc}\mathrm{min}\{q^1N^{1+p_{piv}},(1q)^1N^{1+p_{cp}}\}N^{1+p_{piv}}N^{1.09},$$
(35)
since $`p_{cp}<p_{piv}`$. Since local moves should be the slowest ones we expect $`\tau _{exp}\mathrm{}>N^{1.1}`$.

Finally, we estimate the behaviour of the autocorrelation time for global observables in CPU units. Since the average CPU time spent in performing one move of the cut-and-permute algorithm is proportional to

$$T=qT_{piv}+(1q)T_{cp}qN^{1p_{piv}}+(1q)N^{1p_{cp}}qN^{1p_{piv}},$$
(36)
the autocorrelation time of global observables in CPU units scales as

$`\tau _{int,glob}^{CPU}`$ $``$ $`N^{z_{glob}+1p_{piv}}N^{1.06},`$ (37)
$`\tau _{int,surf}^{CPU}`$ $``$ $`N^{z_{surf}+1p_{piv}}N^{1.40},`$ (38)
where, in the absence of a numerical estimate, we have used our heuristic estimate, $`z_{surf}p_{cp}`$. The behavior of global observables is nearly optimal, while for surface observables the behavior is not as good. If one is interested in computing surface observables, it is possible to improve the algorithm by increasing the frequency of the cut-and-permute moves. Indeed, if we scale $`q^1N^{p_{cp}p_{piv}}`$—essentially we keep the relative frequency of successful pivot and cut-and-permute moves constant—then our previous arguments give

$$\tau _{int,glob}N^{p_{cp}},\tau _{int,surf}N^{p_{cp}}.$$
(39)
The time per move scales now as $`N^{1p_{cp}}`$, so that

$$\tau _{int,glob}^{CPU}\tau _{int,surf}^{CPU}N.$$
(40)
In this case, the dynamics is optimal for both types of observables.

Note that one cannot obtain the same result using the pure pivot algorithm. Indeed, the autocorrelation time for global surface observables would be proportional to the number of attempts occurring between two accepted pivot moves changing the $`z`$-coordinate of the walk and applied in a pivot point $`\omega _k`$ with $`k`$ small. Because of the monomer depletion near the surface, indeed, one does not expect that a move at a pivot point far from the grafted end would affect macroscopically the number of monomers near the surface. In Sec.2 we have shown that among the $`z`$-changing moves only the ones which do not invert $`z`$ have a small probability of being accepted for small $`k`$. This probability decays as $`N^x`$ with $`xp_{cp}0.49`$. Therefore, if we choose the pivot point with uniform probability, we get

$$\tau _{int,surf}^{\mathrm{pivot}}N^{z_{surf}^{\mathrm{pivot}}}N^{1+p_{cp}},$$
(41)
and the autocorrelation time in CPU units scales as

$$\tau _{int,surf}^{CPU,\mathrm{pivot}}N^{z_{surf}^{\mathrm{pivot}}+1p_{piv}}N^{2.4}.$$
(42)
This indicates clearely that the pivot algorithm alone is inefficient for simulating polymers grafted at an inpenetrable interface.

Acknowledgements

The author is glad to acknowledge stimulating discussions with Sergio Caracciolo, Peter Grassberger, Andrea Pelissetto, and Stu Whittington.


format: text

text:
# Spinodal curve of a three-component molecular system

## 1 Introduction

Two almost completely immiscible liquids may increase their mutual solubility in the presence of a solvent. By increasing the solvent concentration, the two liquid layers become more alike and, at some (temperature dependent) composition, a single homogeneous phase is formed. In order to describe this phenomenon, Wheeler and Widom  considered three types of diatomic molecules, $`AA`$, $`BB`$, and $`AB`$, covering the bonds of a regular lattice (one molecule per bond) and subject to the condition that only $`A`$ atoms or only $`B`$ atoms may meet at a given lattice site (infinite repulsion between $`A`$ and $`B`$ atoms of a common site, no interaction otherwise). The model reduces to the standard (spin-1/2, nearest-neighbor interaction) Ising model on the same lattice and its ferromagnetic transition corresponds to a phase-separation transition for the ternary solution. Although the model catches the main aspects of the phase-separation transition, it is not temperature dependent. The model can be extended to include temperature as a significant variable by allowing finite interactions between any two atoms of a common site. This *extended Wheeler-Widom model* is still equivalent to the standard Ising model but on a decorated lattice in which each site of the original lattice is replaced by a cluster of sites. The exact coexistence (binodal) curves of this extended Wheeler-Widom model on the honeycomb and Bethe lattices were determined in  and , respectively. Further extensions of the model have been considered in .

A solution suddenly quenched from an initial one-phase state in the two-phase region reaches its equilibrium through either a nucleation or a spinodal decomposition process . The states below the coexistence curve can be then classified as metastable (close to the binodal) and unstable (deep in the two-phase region), respectively. The *spinodal curve* is supposed to separate these two types of behavior. According to the classical (mean-field) theory of first-order phase transitions, the spinodal is defined as the locus for which an appropriate susceptibility diverges; for a multi-component solution, it follows that the determinant of the diffusion coefficients must be zero and this provides an experimental determination of the spinodal curve. Although a mean-field-like approach oversimplifies the problem, it is still a convenient way to distinguish between the metastable and unstable states.

The aim of the present paper is to determine the spinodal curve for the extended Wheeler-Widom model. The molecules are considered to occupy the bonds of a honeycomb lattice (see Fig. 1) and the lattice of the corresponding Ising model is called 3-12 (see Fig. 2). As shown in Section 2, by applying at first a star-triangle transformation and then a double-decoration, the Ising model on the 3-12 lattice becomes the Ising model on the honeycomb lattice, the original lattice of the molecular model. In Section 3 we shortly review two approximate treatments of the ferromagnetic Ising model: *mean-field* and *Bethe-lattice*. The resulting coexistence and spinodal curves for the three-component model, the effect of temperature, and the limit of a two-component system are all presented in Section 4. In the last section we discuss some possible developments of the present work.

## 2 Extended Wheeler-Widom model on the honeycomb lattice

The grand-canonical partition function of the Wheeler-Widom model with finite interactions on the honeycomb lattice ($`N`$ sites) is proportional to the partition function of the standard Ising model on the 3-12 lattice

$$Z_{\text{3-12}}=\underset{\{\sigma _i\}}{}\mathrm{exp}\left(R\underset{(i,j)C_3}{}\sigma _i\sigma _j+L\underset{(i,j)C_2}{}\sigma _i\sigma _j+h\underset{i}{}\sigma _i\right)$$
(2.1)
with $`\sigma _i,\sigma _j=\pm 1`$ ($`i,j=1,\mathrm{},N`$) depending on the spin orientation. The first sum in (2.1) is over all possible configurations, the second one runs over all the triangles ($`C_3`$ graphs) of the 3-12 lattice, the third one is over all the bonds ($`C_2`$ graphs) connecting the triangles, and the last one is over the whole lattice (see Fig. 2). The parameters $`R`$, $`L`$, and $`h`$ from (2.1) are related to the parameters of the original model as follows

$$R=\frac{1}{4k_BT}(2\epsilon _{AB}\epsilon _{AA}\epsilon _{BB}),$$
(2.2)
$$L=\frac{1}{4k_BT}(\mu _{AA}+\mu _{BB}2\mu _{AB}),$$
(2.3)
and

$$h=\frac{1}{4k_BT}(\mu _{AA}\mu _{BB}2\epsilon _{AA}+2\epsilon _{BB}).$$
(2.4)
Here $`k_B`$ is Boltzmann’s constant and $`T`$ is the temperature. We use $`\epsilon _{AA}`$, $`\epsilon _{AB}`$, and $`\epsilon _{BB}`$ to denote the interaction energies between the molecular ends of a common lattice site. The chemical potentials $`\mu _{AA}`$, $`\mu _{BB}`$, and $`\mu _{AB}`$ all tend to infinity (no vacant sites are allowed), but differences of any two of them are finite thermodynamic variables.

The partition function of the standard Ising model on the 3-12 lattice can be connected to the partition function of the standard Ising model on the honeycomb lattice

$$Z_{\mathrm{hon}}=\underset{\{\sigma _i\}}{}\mathrm{exp}\left(K\underset{<i,j>}{}\sigma _i\sigma _j+3h_{}\underset{i}{}\sigma _i\right)$$
(2.5)
via two transformations: star-triangle followed by double-decoration (see Fig. 3). The result is

$$Z_{\text{3-12}}(R,L,h)=A^N(L_1)B^{3N/2}(L_1,L,h)Z_{\mathrm{hon}}(K,3h_{}),$$
(2.6)
where

$$A^4(L_1)=16\mathrm{cosh}(3L_1)\mathrm{cosh}^3(L_1)$$
(2.7)
and

$$\mathrm{cosh}(2L_1)=(e^{4R}+1)/2,R>0,$$
(2.8)
are obtained from the star-triangle transformation, and

$`B^4`$ $`=`$ $`16e^{4L}\left[1+e^{2L}\mathrm{cosh}(2L_1+2h)\right]\left[1+e^{2L}\mathrm{cosh}(2L_12h)\right]`$ (2.9)
$`\times \left[\mathrm{cosh}(2L_1)+e^{2L}\mathrm{cosh}(2h)\right]^2,`$
$$e^{4K}=\frac{\left[1+e^{2L}\mathrm{cosh}(2L_1+2h)\right]\left[1+e^{2L}\mathrm{cosh}(2L_12h)\right]}{\left[\mathrm{cosh}(2L_1)+e^{2L}\mathrm{cosh}(2h)\right]^2},$$
(2.10)
and

$$e^{4h_{}}=\frac{1+e^{2L}\mathrm{cosh}(2L_1+2h)}{1+e^{2L}\mathrm{cosh}(2L_12h)}$$
(2.11)
follow from the double-decoration transformation. The parameter $`L_1`$ from equations (2.6) to (2.11) is an interaction constant of the Ising model on the intermediate lattice (see Fig. 3); its sign is undetermined and, without loss of generality, we may choose $`L_1>0`$.

Let us denote by $`X_{AA}`$, $`X_{BB}`$, and $`X_{AB}`$ the mole fractions respectively corresponding to the three types of the molecules. The mole fractions obey the conservation equation

$$X_{AA}+X_{BB}+X_{AB}=1$$
(2.12)
and they can be related to the quantities of the equivalent Ising model by

$$X_{AA}+X_{BB}X_{AB}=\sigma _{\text{3-12}}$$
(2.13)
and

$$X_{AA}X_{BB}=m_{\text{3-12}},$$
(2.14)
where $`\sigma _{\text{3-12}}=<\sigma _i\sigma _j>_{C_2}`$ describes the spin-spin correlation on the $`C_2`$ graphs and $`m_{\text{3-12}}=<\sigma _i>`$ is the magnetization of the 3-12 lattice. Using (2.6), $`\sigma _{\text{3-12}}`$ and $`m_{\text{3-12}}`$ can be expressed in terms of the similar quantities $`\sigma `$ and $`m`$ on the honeycomb lattice as

$$\sigma _{\text{3-12}}=\left(\frac{\mathrm{ln}B}{L}\right)_{L_1,h}+\left(\frac{K}{L}\right)_{L_1,h}\sigma +2\left(\frac{h_{}}{L}\right)_{L_1,h}m$$
(2.15)
and

$$m_{\text{3-12}}=\frac{1}{2}\left(\frac{\mathrm{ln}B}{h}\right)_{L_1,L}+\frac{1}{2}\left(\frac{K}{h}\right)_{L_1,L}\sigma +\left(\frac{h_{}}{h}\right)_{L_1,L}m.$$
(2.16)
The Ising model has two parameters, $`K`$ and $`h_{}`$, while the Wheeler-Widom model has three independent combinations of the original parameters: $`R`$, $`L`$, and $`h`$. One can choose $`1/R`$ as the *reduced temperature* of the molecular model ; we will use $`t`$ to denote it and keep it constant throughout the calculation. (When $`\epsilon _{AB}\mathrm{}`$, then $`t0`$ and the model reduces to the original model introduced by Wheeler and Widom ; see also Appendix.) Because we describe the molecular model by means of the Ising model, all the derivatives from the equations (2.15) and (2.16) must be expressed in terms of $`R`$, $`K`$, and $`h_{}`$ rather than $`L_1`$, $`L`$, and $`h`$. The interaction constant $`L_1`$ as a function of $`R`$ is given by (2.8) and from (2.10) and (2.11) we can find $`L`$ and $`h`$ as functions of $`K`$ and $`h_{}`$ at constant $`L_1`$ (or, equivalently, at constant $`t`$). Namely,

$$\mathrm{tanh}(2h)=\frac{\mathrm{sinh}(2L_1)\mathrm{sinh}(2h_{})}{\mathrm{cosh}(2L_1)\mathrm{cosh}(2h_{})e^{2K}}\chi $$
(2.17)
and

$$\frac{e^{2L}}{\sqrt{1\chi ^2}}=\frac{\mathrm{cosh}(2L_1)\mathrm{cosh}(2h_{})e^{2K}}{\mathrm{cosh}(2L_1)e^{2K}\mathrm{cosh}(2h_{})}\lambda .$$
(2.18)
The expressions (2.17) and (2.18) are well defined only if $`|\chi |<1`$ and $`\lambda >0`$, which amounts to

$$|h_{}|<L_1$$
(2.19)
and

$$\frac{1}{2}\mathrm{ln}[\mathrm{cosh}(2L_12|h_{}|)]<K<\frac{1}{2}\mathrm{ln}\left[\frac{\mathrm{cosh}(2L_1)}{\mathrm{cosh}(2h_{})}\right]K_{\mathrm{max}},$$
(2.20)
the latter inequalities determining the region in the $`(K,h_{})`$-space which can be mapped into $`(L,h;R)`$-space. This region becomes the whole space when $`R\mathrm{}`$. The derivatives from (2.15) and (2.16) as functions of $`K`$, $`h_{}`$, and $`R`$ are given in Appendix. Provided $`\sigma `$ and $`m`$ (both functions of $`K`$ and $`h_{}`$) are known, we have now the complete “dictionary” between our molecular system and the Ising model.

The ferromagnetic transition of the Ising model ($`K>K_c>0`$, $`h_{}=0`$) corresponds to a phase-separation transition for the Wheeler-Widom model ($`L>0`$, $`h=0`$). If $`\sigma _0(K)\sigma (K,0)`$ and the spontaneous magnetization $`m_0(K)m(K,0)`$ from (2.15) and (2.16) are specified, then one can determine the isothermal coexistence curves for our molecular system by eliminating $`K`$ between the equations (2.13) and (2.14). The coexistence curves based on the exact solution of the Ising model on the honeycomb lattice and on the three-coordinated Bethe lattice were determined in the  and , respectively. In finding the spinodal curve, we need $`\sigma _s(K)\sigma [K,h_{}^{(s)}(K)]`$ and $`m_s(K)m[K,h_{}^{(s)}(K)]`$ calculated along a certain path $`h_{}^{(s)}(K)`$ of the Ising model. We consider two approximate methods below: mean-field and Bethe-lattice.

## 3 Ising ferromagnet in the mean-field approximation and on the Bethe lattice

In the mean-field approximation, the magnetization of the Ising model with the partition function (2.5) is the solution of the equation

$$m=\mathrm{tanh}(qKm+3h_{}),$$
(3.1)
where $`q`$ is the coordination number of the lattice. This equation is invariant under the transformation $`(m,h_{})(m,h_{})`$, i.e., $`h_{}`$ is an odd function of $`m`$. The spin-spin correlation function $`\sigma `$ can be determined, for example, by using its relation to the internal energy per spin $`u`$,

$$\sigma =\frac{2}{qK}\left(\frac{u}{k_BT}+3h_{}m\right),$$
(3.2)
and substituting the expression for $`u`$ from  into (3.2). As a result,

$$\sigma =m^2.$$
(3.3)
In zero magnetic field and for $`K>K_c=1/q`$ the model exhibits a ferromagnetic transition. The spontaneous magnetization is given by (3.1) with $`h_{}=0`$: if $`m_0`$ is a solution, then $`m_0`$ is also a solution and both these values minimize the free energy as function of the magnetization. Once having determined the spontaneous magnetization, the relation (3.3) gives $`\sigma _0=m_0^2`$. At constant $`K`$, the “van der Waals loops” of $`h_{}`$ as function of $`m`$ (for $`|m|<|m_0|`$) are interpreted within mean-field-like theories as non-equilibrium states: metastable for $`|m|>|m_s|`$ and unstable for $`|m|<|m_s|`$, where $`\pm m_s`$ give the extremum points of $`h_{}`$ (which, at the same time, are the inflexion points of the free energy as a function of the magnetization). By changing $`K`$, these two extremum points describe the spinodal curve, and the mean-field value of $`m_s`$ is

$$m_s=\pm \sqrt{1\frac{K_c}{K}}\mathrm{with}K_c=\frac{1}{q}.$$
(3.4)
The spin-spin correlation function along the spinodal curve is then given by $`\sigma _s=m_s^2`$.

In an improved scheme, with the same results as the Bethe approximation , we consider the Ising model on the Bethe lattice (see Fig. 4). Defining

$$z=e^{2K},\mu =e^{6h_{}},\mathrm{and}\mu _1=\mu x^{q1},$$
(3.5)
the magnetization of the Ising model on the Bethe lattice with the coordination number $`q`$ is

$$m=\frac{1\mu _1^2}{1+2z\mu _1+\mu _1^2},$$
(3.6)
while the spin-spin correlation function is

$$\sigma =\frac{12z\mu _1+\mu _1^2}{1+2z\mu _1+\mu _1^2},$$
(3.7)
where $`\mu _1`$ obeys the equation

$$\frac{\mu _1}{\mu }=\left(\frac{z+\mu _1}{1+z\mu _1}\right)^{q1}.$$
(3.8)
From (3.6) and (3.7) we can get $`\sigma `$ as a function of $`z`$ and $`m`$,

$$\sigma =12z\frac{1m^2}{z+\sqrt{1(1z^2)m^2}}.$$
(3.9)
This is the analog of Eq. (3.3) of the mean-field approximation and, for $`z1`$ ($`K0`$, i.e., small couplings), they coincide. In zero magnetic field and for $`z<z_c=12/q`$ the model undergoes a ferromagnetic transition. Both $`m_0`$ and $`\sigma _0`$ can be determined in principle by solving at first (3.8) with $`\mu =1`$ and then substituting the solution $`\mu _1(z)`$ into (3.6) and (3.7). In the particular case of $`q=3`$ one gets

$$m_0=\pm \frac{1}{12z}\left(\frac{13z}{1+z}\right)^{1/2}$$
(3.10)
and

$$\sigma _0=1+2z\left[\frac{1}{1z}\frac{1}{3(1+z)}\frac{2}{3(12z)}\right].$$
(3.11)
In determining the spinodal curve, we should first derive $`\mu `$ as a function of $`m`$ using (3.6) and (3.8), and then one readily shows that $`h_{}`$ has extremum points at

$$m_s=\pm \sqrt{\frac{1z^2/z_c^2}{1z^2}}\mathrm{with}z_c=1\frac{2}{q}.$$
(3.12)
Combining (3.12) with (3.9) we easily find $`\sigma _s`$.

## 4 Phase diagrams

The ferromagnetic transition of the Ising model occurs in zero magnetic field and for a coupling constant $`K`$ greater than a critical one $`K_c`$. This transition has a correspondent in the molecular model only if $`K_{\mathrm{max}}`$ from (2.20), calculated at $`h_{}=0`$, is greater than $`K_c`$. The condition $`\mathrm{ln}[\mathrm{cosh}(2L_1)]=2K_c`$ then determines the *upper consolute temperature* $`t_u`$ of the model, i.e., the maximum temperature at which a phase separation transition can occur in the solution,

$$t_u=\frac{4}{\mathrm{ln}\left[2e^{2K_c}1\right]}.$$
(4.1)
In the mean-field approximation of the Ising model ($`q=3`$) we get $`t_u3.762`$; for the Bethe lattice, one has $`t_u=4/\mathrm{ln}52.485`$.

In drawing the isothermal coexistence and spinodal curves, we use the triangular diagram: a fixed composition is represented by a point inside an equilateral triangle of the unit height, the mole fraction of each component being given by the distance from that point to the corresponding side of the triangle. For a triangular diagram, the conservation equation (2.12) is automatically fulfilled, and the correspondence with the three components of the molecular solution is chosen as shown in Fig. 5.

Both the coexistence and the spinodal curves are symmetric under reflections around the ($`X_{AA}=X_{BB}`$)-axis, a consequence of the symmetry properties of the derivatives from (2.15) and (2.16) when $`h_{}h_{}`$ (see Appendix). The analog of the Curie point for a ferromagnet is the *plait point*; it is given by the limit $`KK_c`$ and it lies always on the ($`X_{AA}=X_{BB}`$)-axis. At the plait point, the spinodal curve is tangent to the coexistence curve. In the vicinity of the plait point, both the coexistence and the spinodal curves can be approximated by parabolas. In the opposite limit, at $`K=K_{\mathrm{max}}`$, it is shown in Appendix that $`X_{AB}=0`$ and, thus, both the coexistence and spinodal curves reach the base of the triangular diagram.

The phase diagrams at some particular values of the reduced temperatures are presented in Figs. 5 and 6, respectively corresponding to the mean-field and Bethe-lattice approximations of the equivalent Ising model. In the mean-field approximation (Fig. 5), at $`t=0`$ the spinodal and the coexistence curve coincide; they are described by the same parabola $`y=(13x^2)/2`$. This means that no metastable states are predicted at zero temperature, an expected result from more general considerations. Technically, this is a consequence of the relation (3.3) between $`\sigma `$ and $`m`$, which is independent of the coupling constant $`K`$: the distinct $`K`$-dependencies of $`m`$ for the spinodal and coexistence curves play only a scaling role (the same curve is covered in different ways). In the Bethe-lattice case (Fig. 6), the relation (3.9) between $`\sigma `$ and $`m`$ explicitly depends on the coupling constant, and, in general, it should be so. Consequently, the three-component system can be in a metastable state even at zero temperature, a fact that might reveal an inconsistency of the theory. (The implications of this result and a possible way to overcome the difficulty within mean-field-like theories of first-order phase transitions will be discussed in a forthcoming paper.) By increasing the temperature, the curves go down and shrink to a point when $`t=t_u`$. The limiting case $`X_{AB}=0`$ corresponds to a two-component system. In the mean-field approximation of the equivalent Ising model, the spinodal and the coexistence curves depend on temperature as in Fig. 7; the case of the Bethe lattice is illustrated in Fig. 8.

## 5 Discussion

The Wheeler-Widom model (with finite, two-body interactions) for three-component molecular systems is equivalent to the standard (spin-1/2, nearest-neighbor interaction) Ising model. Based on this equivalence, we derived the spinodal curve for the molecular model on a honeycomb lattice using approximate (mean-field and Bethe-lattice) results for the Ising model. Both the spinodal and the (previously determined) coexistence curves have been drawn at different temperatures.

The present work can be extended in several ways. One direction is to consider *three-body* interactions between the molecular ends associated with the same site of the honeycomb lattice. The coexistence curves then become *asymmetric*, with the plait point located off the ($`X_{AA}=X_{BB}`$)-axis of the triangular diagram ; a similar asymmetry should be observed for the corresponding spinodal curves. Another possible development is to determine the spinodals for systems with both upper and *lower* consolute temperatures. Within the present formalism, a lower bound on the transition temperature follows from the assumption that the molecular ends may form *bonds*; the coexistence surface for a three-coordinate Bethe lattice has been determined in . A more general situation, with two possible transitions (one with both upper and lower consolute temperatures, the second one only with an upper consolute temperature below the minimum transition temperature of the first one) has been analyzed in . It is also possible to generate phase diagrams with all the characteristics mentioned above (asymmetric, upper and lower consolute temperatures, etc.) within a single ternary solution model of the Wheeler-Widom type .

In the framework of the formalism presented above, the spinodal curve of a ternary solution is determined starting from a simple microscopic model reducible to two components. Together with analysis of experimental data, it could contribute to a better understanding of interesting phenomena such as diffusion and nucleation in multi-component systems. Another important aspect should also be mentioned: since a Wheeler-Widom-type model is equivalent to the standard Ising model, it perhaps provides the only way to get exact results for systems with three components.

### Acknowledgements

We thank Dale A. Huckaby for many valuable discussions on this subject. The support of the NASA Microgravity Biotechnology Program through Grant NAG8-1356 is gratefully acknowledged.

## Appendix

Let us denote the derivatives that appear in (2.15) and (2.16) by

$$\{\begin{array}{cc}c_{11}=\left(\frac{\mathrm{ln}B}{L}\right)_{L_1,h}\hfill & ,c_{21}=\frac{1}{2}\left(\frac{\mathrm{ln}B}{h}\right)_{L_1,L}\hfill \\ & \\ c_{12}=\left(\frac{K}{L}\right)_{L_1,h}\hfill & ,c_{22}=\frac{1}{2}\left(\frac{K}{h}\right)_{L_1,L}\hfill \\ & \\ c_{13}=2\left(\frac{h_{}}{L}\right)_{L_1,h}\hfill & ,c_{23}=\left(\frac{h_{}}{h}\right)_{L_1,L}.\hfill \end{array}$$
(A.1)
Using the notations

$$p=\left(e^{4R}+1\right)/2\mathrm{and}z=e^{2K},$$
(A.2)
the quantities $`\chi `$ and $`\lambda `$ respectively defined in (2.17) and (2.18) become

$$\chi =\frac{\sqrt{p^21}\mathrm{sinh}(2h_{})}{p\mathrm{cosh}(2h_{})z}\mathrm{and}\lambda =\frac{p\mathrm{cosh}(2h_{})z}{pz\mathrm{cosh}(2h_{})}.$$
(A.3)
Introducing

$$\omega =(1+p\lambda )^2(p^21)\chi ^2\lambda ^2,$$
(A.4)
the $`c`$-quantities from (A.1) are given in terms of $`K`$, $`h_{}`$, and $`R`$ as

$$\{\begin{array}{cc}c_{11}=\frac{\lambda }{p+\lambda }\frac{1+p\lambda }{\omega }\hfill & ,c_{21}=\frac{1+z^2}{2}\frac{\chi \lambda }{p+\lambda }\hfill \\ & \\ c_{12}=\frac{p}{p+\lambda }\frac{1+p\lambda }{\omega }\hfill & ,c_{22}=\frac{1z^2}{2}\frac{\chi \lambda }{p+\lambda }\hfill \\ & \\ c_{13}=2\sqrt{p^21}\frac{\chi \lambda }{\omega }\hfill & ,c_{23}=\sqrt{p^21}[1+p(1\chi ^2)\lambda ]\frac{\lambda }{\omega }.\hfill \end{array}$$
(A.5)
Taking the limit $`R\mathrm{}`$ in (A.5), we get $`c_{11}=c_{13}=c_{21}=c_{22}=0`$ and $`c_{12}=c_{23}=1`$. This limit, understood as $`\epsilon _{AB}\mathrm{}`$ \[see (2.2)\], corresponds to the original model introduced by Wheeler and Widom .

Let us remark that $`\chi `$ is an odd function of $`h_{}`$, while $`\lambda `$ and $`\omega `$ are even functions of $`h_{}`$ \[see (A.3) and (A.4)\]. As a consequence, the quantities $`c_{11}`$, $`c_{12}`$, and $`c_{23}`$ are symmetric under the transformation $`h_{}h_{}`$, whereas $`c_{13}`$, $`c_{21}`$, and $`c_{22}`$ are antisymmetric. Because $`m`$ is antisymmetric in $`h_{}`$ and $`\sigma `$ is an even function of $`m`$ (see Section 3), from the mole fractions equations (2.12) to (2.14) it follows that $`X_{AB}`$ is symmetric in $`h_{}`$, while $`X_{AA}X_{BB}`$ is antisymmetric. These properties determine the symmetry of the phase diagrams (discussed in Section 4) under reflections around the ($`X_{AA}=X_{BB}`$)-axis. Moreover, at $`h_{}=0`$ all $`c_{13}=c_{21}=c_{22}=0`$ and $`\sigma _{\text{3-12}}`$ decouples from $`m_{\text{3-12}}`$ \[see (2.15) and (2.16)\].

Another interesting limit is when $`KK_{\mathrm{max}}`$ \[see (2.20)\]. In this case $`\lambda \mathrm{}`$, $`\omega `$ goes to infinity as $`\lambda ^2`$, and from (A.5) we get $`c_{11}=1`$, $`c_{12}=c_{13}=0`$; consequently, $`\sigma _{\text{3-12}}`$ given by (2.15) is equal to one. From (2.12) and (2.13) it then follows that $`X_{AB}=0`$. This property is used in Section 4 to show that in the limit $`KK_{\mathrm{max}}`$ the coexistence and the spinodal curves always reach the bottom side of the triangular diagram.