- Email: [email protected]
Stability of the hexagonal lattice of charged colloids
Journal of Molecular Liquids 131–132 (2007) 173 – 178 www.elsevier.com/locate/molliq
Stability of the hexagonal lattice of charged colloids J. Dobnikar a,b,⁎, P. Ziherl b,c a
Institute of Chemistry, Karl-Franzens Universität Graz, Heinrichstrasse 28, A-8010 Graz, Austria b J. Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia c Department of Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia Available online 25 September 2006
Abstract We discuss the stability of the hexagonal lattice in charged colloidal systems using (i) an effective Yukawa potential characterized by a cutoff introduced to model many-body effects and (ii) the full Poisson–Boltzmann cellular theory. We analyze the hexagonal-to-square structural transition in terms of a simple shear mode, and we scan the phase diagram at T = 0. We find that a short enough cutoff destabilizes the hexagonal lattice, while the more complete Poisson–Boltzmann theory does not show this feature. By expanding the energy in terms of the order parameter of the transition, we identify the main structural differences between the two lattices and we show that the hexagonal lattice could be destabilized by a shoulder-like potential; however, electrostatic interaction is not expected to give rise to an effective potential of this kind. © 2006 Elsevier B.V. All rights reserved. Keywords: Charged colloids; Crystal lattices; Phase diagram; Colloidal interactions
1. Introduction Understanding crystal structure is difficult. Even when the interparticle potential is as simple as the Coulomb interaction, finding the minimal-energy distribution of particles is far from trivial. The notoriously hard reduced-dimensionality examples include the well-studied Wigner crystal of electrons on a plane and the Thomson problem of the equilibrium distribution of point charges on unit sphere (solved for 2–8 particles, still unsolved for 9 and 11 particles). In three dimensions, things are even more complex but more interesting and important at the same time. In soft matter and biological physics, electrostatic interaction is of paramount importance: For example, many colloids are stabilized against flocculation and coagulation by charge alone. In these systems, the analysis is simplified because they are purely classical but more complicated by the finite temperature. The key feature of the Coulombic interaction that makes the theoretical analysis of charged systems difficult is its long range. Thus it seems that systems with short-range interaction are more transparent as there each particle has fewer neighbors to interact with. However, for short-range interaction potentials the situation ⁎ Corresponding author. J. Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia. E-mail address: [email protected] (J. Dobnikar). 0167-7322/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2006.08.034
is more complicated because in this case the relative weight of the entropic part of the free energy becomes larger and larger as the temperature is increased. At large enough temperatures, the system then behaves essentially as an assembly of hardcore particles. One would expect that some insight into the phase diagram could be obtained by T = 0 calculation. More precisely, such an analysis would yield the candidate structures that could but need not appear in the full (T, p) phase diagram. This was the motivation for several T = 0 studies of the phase behavior of monodisperse particles characterized by particularly simple interactions. Most of the studies were limited to 2D systems, which already exhibited a number of very interesting phases. For example, in a hard-core–linear-ramp potential [1], square, oblique, and more complicated phases were found. In a related study, a family of short-range potentials was studied and the crystalline phases numerically observed were all very interesting [2,3]. The common denominator of all these studies is that a weak, finite-range repulsive potential tail attached to the hard core can render the hexagonal phase unstable. At the same time, it is known that for long-range potentials such as the Coulomb potential, the minimal-energy configuration is hexagonal. So if one would plot the phase diagram as a function of some parameter characterizing the range of the potential, one or more non-hexagonal phases would be expected in between the hexagonal phase at very short and at very long range potentials.
174
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
Needless to say, transitions between these phases should be accompanied by changes of thermodynamic, elastomechanical, and optical properties of the system. At fixed density and temperature, such a transition could be most easily induced in charged particles in an electrolyte. In this paper, we would like to explore the possibility that non-closepacked lattices be stable in such a system on top of the hexagonal phase. This question is not new, but so far it has only been addressed within the Yukawa approximation [4]. By now, it is well-established that this formalism does not account properly for all the details of the true screened interaction. As the density is increased, many-body effects gradually become more and more important. Here we re-address the symmetry of two-dimensional crystal lattices of a charged system by describing the pair interaction using the truncated Yukawa potential, which has been shown to successfully capture some many-body effects [5], and the more accurate Poisson–Boltzmann theory. We adopt a conservative approach and compare the candidate lattices at T = 0. We find that within the Poisson–Boltzmann theory, the hexagonal lattice is stable at all screening lengths, thereby reproducing the qualitative results obtained for Yukawa pair potential [4].
Lennard–Jones interaction [7] and some features of the more complicated non-close-packed lattices are also present in systems with attractive rather than repulsive short-range potentials [9]. Given that the true electrostatic potential does behave as if it had a well-defined cutoff and that potentials with a cutoff can stabilize non-close-packed lattices, our main goal in this paper is to answer the question: Is there a parameter range for which the ground state of a charged two-dimensional colloidal system is a lattice other than hexagonal? To this end, we have performed a T = 0 analysis. The reason why we chose this approach is that it provides a good estimate for the symmetry of the candidate structures to be found at finite temperatures at reasonable computational costs. Although we cannot rule out the possibility that some phases appear at finite temperatures but not at T = 0, all the evidence on similar systems tells us that this would be very unusual [1]. We have systematically analyzed all lattices with 1 particle per unit cell and a few other lattices. If a pairwise additivity is assumed, the energy per particle in a given lattice configuration is
2. Non-close-packed lattices in two dimensions
where U(r) is the pair interaction betweenPthe particles, rij the distance between lattice points i and j and L the sum over the lattice points of a given lattice L. In order to analyze the phase diagram it suffices to know the relative energies of the different lattices. A suitable choice of the unit of energy is the energy per particle in a hexagonal lattice with nearest-neighbor interactions qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi only, EH/N = 3U(dNN); here dNN ¼ 2=q 3 is the lattice spacing in the hexagonal lattice. — The reduced energies of the lattices are then defined as eL = (EL⁎/N)/(EH/N). In what follows, we are going to use dNN as a unit of length.
The search for non-close-packed lattices in charged colloidal systems is mainly motivated by the well-known result that in three dimensions, an FCC crystal of point-like Yukawa particles undergoes a structural transition to the BCC phase if the screening is weak enough, i.e., if the screening length is large [6]. This result has been corroborated to include many-body effects [5] which do modify the phase boundaries to some extent, but the qualitative conclusions remain the same. In two dimensions, this question has been addressed almost as early as in three dimensions [4], and it was found that for Yukawa potential the hexagonal lattice is the most stable structure at all screening lengths [4]. Thus no structural transition can be induced solely by varying the salt concentration and thus the screening length. However, the Yukawa potential is only an approximation of the true interparticle potential, which is valid in the limit of infinitely small density. The true screened electrostatic potential can depart from the Yukawa potential considerably, and some of these effects can be described by model potentials such as the truncated Yukawa potential with a density-dependent cutoff [5]. The presence of a cutoff in this effective potential is a crucial point, as short-range potentials are known to stabilize non-closepacked lattices [1,2,3,7,8]. The numerical studies of both T = 0 and T N 0 systems cover a wide spectrum of potentials ranging from the simple hard-core–soft-shoulder [8] and hard-core–linear ramp potentials [1] to the more general hard-core–short-tail potentials [2]. In these systems, stripe, square, and oblique lattice as well as more complicated 2D lattices with more than 1 particle per unit cell were observed. Another often encountered phenomenon is clustering or aggregation where particles form dimers or trimers (i.e., colloidal molecules) [8]. Similar behavior was found in repulsive long-range but shoulder-like potential based on 6–12
L X L L EL4 1X 1X ¼ U ðrij Þ ¼ U ðrij Þ; N i jpi 2 jpi N
ð1Þ
3. Truncated Yukawa potential The interaction between two isolated charged colloids in an electrolyte suspension is well-described by the DLVO theory [10]. This theory predicts screening of the electric charge on the colloids by counterions and effective interactions of Yukawa type, U(r) = U0exp(−κ⁎r)/r, where 1/κ⁎ is the effective screening length and U0 is the magnitude of the potential. In the present context we do not need to include the van der Waals term – which is also a part of the theory – because it is unimportant if only equilibrium properties of colloidal crystals are studied. Due to the strong electrostatic repulsion the mean separation between colloids is namely always much larger than the range of the van der Waals interaction. When there are more than two colloids in the system, i.e., when we analyze a dense colloidal suspension, pairwise additivity cannot be generally assumed due to the presence of many-body interactions. It has been observed in Ref. [11] and in Refs. [5,12] that the effective potential between colloids in a dense system deviates from the simple Yukawa form and that, due to the many-body interactions, it depends on the density and other state parameters of the system. In Refs. [5,12] the solid– liquid melting line in the three-dimensional colloidal suspension
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
Fig. 1. Construction of the lattices.
has been calculated by combining the solutions of the nonlinear Poisson–Boltzmann equation with Brownian dynamics simulations. The Poisson–Boltzmann theory correctly accounts for all many-body interactions. The resulting phase diagram has been compared to the phase diagram for a pairwise interacting Yukawa system [13] and systematic deviations have been observed. This effect of many-body interactions has been relatively successfully modeled by introducing a density-dependent cutoff of the effective Yukawa potentials such that the potential is simply truncated at a given distance dC. If dC is between the first and second neighbor shell, the truncated Yukawa potential only includes the nearest neighbor interactions. Using such effective pair interactions, the many-body phase diagram has been reproduced reasonably well. This was our motivation to study the phase diagram of 2D colloidal crystals interacting via the truncated Yukawa interaction. We have already fixed the unit of length to dNN. The parameter κ⁎ rescales to κD = κ⁎dNN. Additionally, the unit of energy is chosen to be the energy per particle in a hexagonal lattice with nearest neighbor interactions only, EH/N = 3U0exp (− κD). The energy per particle in a 2D colloidal crystal interacting via the truncated Yukawa interaction is eL ðjD ; dC Þ ¼
L 1 X expð−jD ðrij −1ÞÞ : 6 j; rij b d rij
175
unit cell are generated. For example, a hexagonal phase is described by a =ffiaH = 1 and Δx = 1/2, while for square phase it is qp ffiffiffiffiffiffiffiffiffiffiffi ffiffiffi 3=2c0:93 and Δx = 0. We have systematically a ¼ aS ¼ varied both a and Δx and we have evaluated the energy (Eq. (2)) for a wide range values of κD and dC in each lattice. At small values of κD, the screening length is large and the interaction in relatively long-ranged. When κD is large, the interaction is short-ranged and at very high values the hard-core limit is approached. We find the hexagonal phase to be energetically most favorable at large κD values for all values of the cutoff dC. The value of dC can be varied between 1 (which corresponds to nearest-neighbor interactions only) and ∞ (which represents the simple Yukawa limit). Anything in between is simply an interaction terminated at a certain distance. Again, an infinite cutoff reproduces the simple Yukawa system where the hexagonal lattice is the stable phase for all κD values [4]. However, at the cutoff value around 1, where only the nearest neighbor interactions are considered, we find other phases to be stable at small values of κD. The calculated phase diagram is shown in Fig. 2. We see two regions where the square phase has the lowest energy and many regions where a deformed hexagonal phase was found to have the smallest energy among the candidate lattices considered. The deformed hexagonal phases are denoted by the letter H followed by a number denoting the value of a, while Δx is equal to 0.5 for all H-phases. We have not studied the mechanical stability of the lowest energy lattices in the phase diagram, so it could well be that the phases shown here are in fact not stable and therefore not the real ground states of the system. What we consider the main message of the text is, however, that the hexagonal phase is destabilized by the cutoff in the effective interactions. One must be aware that due to the discrete structure of the lattice distances the phase diagram is not a smooth function of the parameter dC. If we scanned the phase diagram in more detail, it would exhibit a more elaborate structure, possibly
ð2Þ
C
The energy eL depends on the lattice type L and on two parameters, κD and dC. 3.1. Phase diagram The general crystal lattice in 2D with a single-particle per unit cell is described by three parameters. Once the density is fixed, the remaining two parameters are the angle between the base vectors and the ratio of their lengths. We have constructed the lattices as illustrated in Fig. 1. The parameters of a lattice are the particle distance along the row a and the horizontal displacement of two consecutive rows Δx. The vertical spacing pffiffiffiffiffiffiffiffiffiffi between the rows, Δy, is fixed by density, Dy ¼ 3=2a. By varying a and Δx all possible 2D lattices with one particle per
Fig. 2. Truncated Yukawa potential: Phase diagram. The symbols H followed by apnumber denote a deformed hexagonal lattice (see main text), CR represents the ffiffiffi 3R30- triangular commensurate lattice [14]. We see the regions of square and CR lattices having smaller energy than the hexagonal lattice.
176
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
including other distinct phases on top of those plotted in Fig. 2. However, it was not our intention to deliver a precise phase diagram for truncated Yukawa interaction as it depends on the exact shape of the cutoff which is to a large extent arbitrary although the existence of the cutoff itself is not. Given this reservation, the phase diagram includes all lattices generated by the described procedure, which have one particle per unit cell, and a few lattices with more ffithan one particles per unit cell, like pffiffiffiffiffiffiffiffiffiffiffiffi the CR structure ( 3R30- triangular commensurate lattice [14]). As a word of caution, we should stress that it is possible that some more complicated and exotic lattices with several particles per unit cell could also intervene and replace some of the phases that are present in the phase diagram. But even if this would be the case, what we find to be most important is that for truncated Yukawa potential, we do find phases other than hexagonal to be the ground state for a range of cutoffs. 4. Poisson–Boltzmann cellular model The energy of a particle in a two-dimensional colloidal crystal can be obtained by solving the nonlinear Poisson– Boltzmann equation within a unit cell with periodic boundary conditions at the cell boundaries such that the normal derivative of the potential at the boundary is zero. Let us denote the unit cell of lattice L by GL. GL is the volume of the unit cell outside the colloid, i.e., the solvent volume. The grand canonical Poisson–Boltzmann equation reads j2 w ¼ j2 sinhðwÞ
for
Y raGL ;
¼−
ZkB R2
for
Y raAGc ;
ð4Þ
where R is the radius of the colloid. The second boundary condition should account for the interaction of the central particle with all other particles in the crystal. It is not difficult to see that this is taken care of by demanding that the normal derivative of the potential at the cell boundary ∂GL be 0, which is nothing but the method of images applied in the reverse direction. Eq. (3) together with the two boundary conditions represents a complete description of the electrostatic problem. Once the potential ψ(r) is known, the grand potential energy can be computed as follows 1 ZkB 4pkB bX ¼ 2 R
Z AG Zc
wdr
þ ðjRÞ2
ðw sinhw−2 coshw þ 2Þdr GL
In order to systematically study the hexagonal and the square lattice and the lattices in between them, we partition the crystal into rectangular cells containing 2 rather than 1 particle. We denote the sides of the rectangle pffiffiffiffiffiffiffiffi by a and b, the latter being fixed by the density, b ¼ 3=a; all distances are measured in units of nearest neighbor distance in the hexagonal lattice. The positions of the two particles within the cell are denoted by r1 and rq2.ffiffiffiffiffiffiffiffiffiffiffiffi Let us put δ = Δx/a. If a = aS + 2δ(1 − aS) where pffiffiffiffiffiffiffiffiffi aS ¼ 3=2, then the shear mode
að1 þ dÞ=2 ; b=4 að1−dÞ=2 r2 ¼ 3b=4
r1 ¼
ð6Þ
continuously transforms the hexagonal phase (δ = 0.5) to the square phase (δ = 0); δ clearly plays the role of the order parameter. We have solved the Poisson–Boltzmann equation in the appropriate unit cell for several values of δ between 0 and 0.5 and we have evaluated the energy (Eq. (5)). We did this for a whole range of the relevant parameters characterizing the solution of the Poisson–Boltzmann equation (the parameters are ZλB, κa and the volume fraction of the colloids). The hexagonal phase always had the lowest energy and in contrast to the results from the truncated Yukawa interactions, we did not see any other phases to be stable.
ð3Þ
where ψ = eϕ/κBT is the reduced electrostatic potential and κ2 = 4πλBcS; here λB = e2/4πϵϵ0κBT is the Bjerrum length and cS is the reservoir monovalent salt concentration. The first boundary condition is determined by the fixed surface charge of the colloids. By denoting the surface of the colloid by ∂Gc, and its normal vector pointing into the solvent by Y nˆ p, this boundary condition can be written as Y nˆ p djw
4.1. Hexagonal to square lattice
ð5Þ
5. Stability Another point of view is to analyze the stability of the hexagonal lattice with respect to shear deformations. This is usually done by parameterizing the transition from the initial to the final crystal configuration by a suitably chosen order parameter, such as the Bain strain which transforms a BCC lattice into an FCC lattice via a continuous sequence of BCT lattices. Of course, there are infinitely many shear modes whose complexity grows with the size of the final configuration, which may contain more than 1 particle per unit cell. But for the simple lattices, say a square lattice, it is plausible to expect that the mode will also be simple and amenable to analysis. The parameterization of the shear modes is not unique; the only constraint they must satisfy is that the density is fixed. This program, used in Ref. [1], consists of a harmonic expansion of the total energy in terms of the order parameter ϵ starting with an expansion of the interparticle distances Ari f ri ðeÞ ¼ ri ðe0 Þ þ ðe−e0 Þ Ae e¼e0 1 A2 ri ðe−e0 Þ2 þ : : : ð7Þ þ 2 Ae2 e¼e0 where r i(ϵ) is a suitable parameterization of the structural transformation and ϵ0 corresponds to the configuration
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
under consideration. To second order in ϵ, the energy reads 1X Ari /ðri ð0ÞÞ þ / Vðri ð0ÞÞ ðe−e0 Þ U ðaÞ ¼ 2 i Ae e¼e0 2 1 A ri þ/ Vðri ð0ÞÞ ðe−e0 Þ2 2 Ae2 e¼e0 2 1 W Ari þ / ð r i ð 0Þ Þ ðe−e0 Þ2 2 Ae e¼e0
½
ð8Þ
where i runs over all neighbors of the central particle and ϕ is the pair potential; ϕ′ = dϕ/dr, etc. By assumption, the initial lattice corresponds to an extremum of the energy, so terms linear in ϵ cancel out and we have U ðeÞ ¼ U ðe0 Þ þ
ðe−e0 Þ2 a; 2
ð9Þ
where U(0) = (1/2)∑ i ϕ(r i) is the energy of the initial state and
A2 U au Ae2 ¼
1X 2
i
"
e¼e0
2 2 # A r Ari i / Vðri ðe0 ÞÞ þ/W ðri ðe0 ÞÞ : Ae2 e¼e0 Ae e¼e0 ð10Þ
If α b 0, the lattice with ϵ = ϵ0 is unstable with respect to the shear mode studied. We stress again that α depends on the parameterization of the shear mode, and so this analysis can but need not indicate that a certain lattice is unstable. A general analysis of the structure of the two terms in Eq. (10) is difficult, and we will limit our discussion to the square and the hexagonal lattice and the shear mode introduced in Section 4. Within this framework, the derivatives can be calculated analytically but the results themselves are not very
177
instructive. However, a numerical evaluation of the derivatives shows that for the two lattices in question, all shell subsums C1 = ∑i∈shell∂2ri/∂ϵ2 are positive. For repulsive potentials, ϕ′ b 0 for all r and the first sum in Eq. (10) will always be negative, and the lattice will only be stable if the second term in Eq. (10) is dominant. Shell subsums C1 and C2 = ∑i∈shell(∂ri/∂ϵ)2 for the hexagonal and the square lattice are shown in Fig. 3. The two lattices do have much in common but there are also some important differences. In general, the magnitudes of C1 and C2 are rather similar, and C1 and C2 both increase with separation when appropriately averaged such that the discreteness of the diagrams is washed out. However, the short-range behavior of the shell subsums in the hexagonal lattice is very different than in the square lattice. One notes immediately that in the hexagonal lattice, C2 is almost constant and rather large within the first six or so shells, whereas in the square lattice it goes to 0 with decreasing separation. This implies that the stability of the hexagonal lattice largely depends on the short-range behavior of the pair potential, and the sign of α (i.e., the stability of the lattice) can be estimated by considering the first shell alone. In the square lattice, many more shells must be taken into account to compute the exact value α and to see whether this lattice is stable. Another consequence of the large relative magnitude of C2/C1 in the first few shells of the hexagonal lattice (which is ≈5 when averaged over the first few shells) is that it can be easily destabilized by a shoulder-like, concave pair potential where ϕ″ b 0 at small separations. On the other hand, the square lattice would not be affected by such potentials so easily as in the first lattice, the magnitude of C2/C1 is about 1. This is, we believe, the mechanism that stabilizes the numerous non-close-packed lattices observed in numerical simulations of the so-called core-softened particles [6]. We stress that the same effect can also be achieved with potentials with a sharp cutoff where the energies of the various lattices differ primarily because of the different number of shells covered by the potential. Although limited in scope, this illustration does point out that the effective pair potential could render the hexagonal lattice
Fig. 3. Shell subsums for the hexagonal lattice (left) and the square lattice (right). As far as C1 is concerned, the two lattices are not very different, but within the first few shells C2 of the hexagonal lattice is much larger than in the square lattice. This gives rise to the robustness of the hexagonal lattice in case of convex repulsive pair interactions such as the Yukawa potential. Note that the scales on the top and bottom panels are not the same.
178
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
unstable not only if it were characterized by an abrupt cutoff but also if it were concave at small separations. Whether this can happen in charged colloidal systems is presently unexplored: The Yukawa potential is strictly convex at all separations, and the many-body corrections to the pair potential should be rather strong to make the curvature of the effective potential negative. 6. Conclusions We have presented the T = 0 energy calculations for 2D crystal structures in systems interacting via the truncated Yukawa pair interaction. We have evaluated the energies for phases continuously varying from the hexagonal to the square lattice. The results show that the cutoff introduced with the intention to capture many-body effects destabilizes the hexagonal lattice at certain parameters values. The destabilization is also consistent with the additionally performed stability analysis in Section 5, where we show that the interaction potential needs to be concave in order to destabilize the hexagonal phase. Although we have only analyzed a single shear mode chosen for its simplicity, we believe that other possible ways of parameterization of the hexagonal-to-square transformation should produce rather similar results. However, we have found the hexagonal phase to be stable at all parameters in the full Poisson–Boltzmann description, which means that the truncated Yukawa description does not correctly capture many-body effects in the system under investigation, although previous studies in 3D [5,12] claim the opposite. This finding shows – once again – that the implementation of effective
interactions in the genuine many-body systems is subtle and should be applied with considerable precautions. Acknowledgments We kindly acknowledge useful discussions with H.H. von Grünberg and O. Glatter. This work was supported by a Marie Curie fellowship (MEIF-CT-2003-501789). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
E.A. Jagla, Phys. Rev., E 58 (1998) 1478. E.A. Jagla, J. Chem. Phys. 110 (1999) 451. G. Malescio, G. Pellicane, Phys. Rev., E 70 (2004) 021202. F.M. Peeters, X. Wu, Phys. Rev., A 35 (1987) 3109. J. Dobnikar, R. Rzehak, H.H.v. Grünberg, Europhys. Lett. 61 (2003) 695. K. Kremer, M.O. Robbins, G.S. Grest, Phys. Rev. Lett. 57 (1986) 2694. P.J. Camp, Phys. Rev., E 68 (2003) 061506. G. Malescio, G. Pellicane, Nat. Mater. 2 (2003) 97. A. Skibinsky, S.V. Buldyrev, A. Scala, S. Havlin, H.E. Stanley, Phys. Rev., E 60 (1999) 2664. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions, Cambridge University Press, 1989. M. Brunner, C. Bechinger, W. Strepp, V. Lobaskin, H.H.v. Grünberg, Europhys. Lett. 58 (2002) 926. J. Dobnikar, Y. Chen, R. Rzehak, H.H. von Grünberg, J. Chem. Phys. 119 (9) (2003) 4971. M.O. Robbins, K. Kremer, G.S. Grest, J. Chem. Phys. 88 (1988) 3286. S. Bleil, H.H. von Grünberg, J. Dobnikar, R. Castañeda-Priego, C. Bechinger, Europhys. Lett. 73 (3) (2006) 450.
Stability of the hexagonal lattice of charged colloids J. Dobnikar a,b,⁎, P. Ziherl b,c a
Institute of Chemistry, Karl-Franzens Universität Graz, Heinrichstrasse 28, A-8010 Graz, Austria b J. Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia c Department of Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia Available online 25 September 2006
Abstract We discuss the stability of the hexagonal lattice in charged colloidal systems using (i) an effective Yukawa potential characterized by a cutoff introduced to model many-body effects and (ii) the full Poisson–Boltzmann cellular theory. We analyze the hexagonal-to-square structural transition in terms of a simple shear mode, and we scan the phase diagram at T = 0. We find that a short enough cutoff destabilizes the hexagonal lattice, while the more complete Poisson–Boltzmann theory does not show this feature. By expanding the energy in terms of the order parameter of the transition, we identify the main structural differences between the two lattices and we show that the hexagonal lattice could be destabilized by a shoulder-like potential; however, electrostatic interaction is not expected to give rise to an effective potential of this kind. © 2006 Elsevier B.V. All rights reserved. Keywords: Charged colloids; Crystal lattices; Phase diagram; Colloidal interactions
1. Introduction Understanding crystal structure is difficult. Even when the interparticle potential is as simple as the Coulomb interaction, finding the minimal-energy distribution of particles is far from trivial. The notoriously hard reduced-dimensionality examples include the well-studied Wigner crystal of electrons on a plane and the Thomson problem of the equilibrium distribution of point charges on unit sphere (solved for 2–8 particles, still unsolved for 9 and 11 particles). In three dimensions, things are even more complex but more interesting and important at the same time. In soft matter and biological physics, electrostatic interaction is of paramount importance: For example, many colloids are stabilized against flocculation and coagulation by charge alone. In these systems, the analysis is simplified because they are purely classical but more complicated by the finite temperature. The key feature of the Coulombic interaction that makes the theoretical analysis of charged systems difficult is its long range. Thus it seems that systems with short-range interaction are more transparent as there each particle has fewer neighbors to interact with. However, for short-range interaction potentials the situation ⁎ Corresponding author. J. Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia. E-mail address: [email protected] (J. Dobnikar). 0167-7322/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2006.08.034
is more complicated because in this case the relative weight of the entropic part of the free energy becomes larger and larger as the temperature is increased. At large enough temperatures, the system then behaves essentially as an assembly of hardcore particles. One would expect that some insight into the phase diagram could be obtained by T = 0 calculation. More precisely, such an analysis would yield the candidate structures that could but need not appear in the full (T, p) phase diagram. This was the motivation for several T = 0 studies of the phase behavior of monodisperse particles characterized by particularly simple interactions. Most of the studies were limited to 2D systems, which already exhibited a number of very interesting phases. For example, in a hard-core–linear-ramp potential [1], square, oblique, and more complicated phases were found. In a related study, a family of short-range potentials was studied and the crystalline phases numerically observed were all very interesting [2,3]. The common denominator of all these studies is that a weak, finite-range repulsive potential tail attached to the hard core can render the hexagonal phase unstable. At the same time, it is known that for long-range potentials such as the Coulomb potential, the minimal-energy configuration is hexagonal. So if one would plot the phase diagram as a function of some parameter characterizing the range of the potential, one or more non-hexagonal phases would be expected in between the hexagonal phase at very short and at very long range potentials.
174
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
Needless to say, transitions between these phases should be accompanied by changes of thermodynamic, elastomechanical, and optical properties of the system. At fixed density and temperature, such a transition could be most easily induced in charged particles in an electrolyte. In this paper, we would like to explore the possibility that non-closepacked lattices be stable in such a system on top of the hexagonal phase. This question is not new, but so far it has only been addressed within the Yukawa approximation [4]. By now, it is well-established that this formalism does not account properly for all the details of the true screened interaction. As the density is increased, many-body effects gradually become more and more important. Here we re-address the symmetry of two-dimensional crystal lattices of a charged system by describing the pair interaction using the truncated Yukawa potential, which has been shown to successfully capture some many-body effects [5], and the more accurate Poisson–Boltzmann theory. We adopt a conservative approach and compare the candidate lattices at T = 0. We find that within the Poisson–Boltzmann theory, the hexagonal lattice is stable at all screening lengths, thereby reproducing the qualitative results obtained for Yukawa pair potential [4].
Lennard–Jones interaction [7] and some features of the more complicated non-close-packed lattices are also present in systems with attractive rather than repulsive short-range potentials [9]. Given that the true electrostatic potential does behave as if it had a well-defined cutoff and that potentials with a cutoff can stabilize non-close-packed lattices, our main goal in this paper is to answer the question: Is there a parameter range for which the ground state of a charged two-dimensional colloidal system is a lattice other than hexagonal? To this end, we have performed a T = 0 analysis. The reason why we chose this approach is that it provides a good estimate for the symmetry of the candidate structures to be found at finite temperatures at reasonable computational costs. Although we cannot rule out the possibility that some phases appear at finite temperatures but not at T = 0, all the evidence on similar systems tells us that this would be very unusual [1]. We have systematically analyzed all lattices with 1 particle per unit cell and a few other lattices. If a pairwise additivity is assumed, the energy per particle in a given lattice configuration is
2. Non-close-packed lattices in two dimensions
where U(r) is the pair interaction betweenPthe particles, rij the distance between lattice points i and j and L the sum over the lattice points of a given lattice L. In order to analyze the phase diagram it suffices to know the relative energies of the different lattices. A suitable choice of the unit of energy is the energy per particle in a hexagonal lattice with nearest-neighbor interactions qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi only, EH/N = 3U(dNN); here dNN ¼ 2=q 3 is the lattice spacing in the hexagonal lattice. — The reduced energies of the lattices are then defined as eL = (EL⁎/N)/(EH/N). In what follows, we are going to use dNN as a unit of length.
The search for non-close-packed lattices in charged colloidal systems is mainly motivated by the well-known result that in three dimensions, an FCC crystal of point-like Yukawa particles undergoes a structural transition to the BCC phase if the screening is weak enough, i.e., if the screening length is large [6]. This result has been corroborated to include many-body effects [5] which do modify the phase boundaries to some extent, but the qualitative conclusions remain the same. In two dimensions, this question has been addressed almost as early as in three dimensions [4], and it was found that for Yukawa potential the hexagonal lattice is the most stable structure at all screening lengths [4]. Thus no structural transition can be induced solely by varying the salt concentration and thus the screening length. However, the Yukawa potential is only an approximation of the true interparticle potential, which is valid in the limit of infinitely small density. The true screened electrostatic potential can depart from the Yukawa potential considerably, and some of these effects can be described by model potentials such as the truncated Yukawa potential with a density-dependent cutoff [5]. The presence of a cutoff in this effective potential is a crucial point, as short-range potentials are known to stabilize non-closepacked lattices [1,2,3,7,8]. The numerical studies of both T = 0 and T N 0 systems cover a wide spectrum of potentials ranging from the simple hard-core–soft-shoulder [8] and hard-core–linear ramp potentials [1] to the more general hard-core–short-tail potentials [2]. In these systems, stripe, square, and oblique lattice as well as more complicated 2D lattices with more than 1 particle per unit cell were observed. Another often encountered phenomenon is clustering or aggregation where particles form dimers or trimers (i.e., colloidal molecules) [8]. Similar behavior was found in repulsive long-range but shoulder-like potential based on 6–12
L X L L EL4 1X 1X ¼ U ðrij Þ ¼ U ðrij Þ; N i jpi 2 jpi N
ð1Þ
3. Truncated Yukawa potential The interaction between two isolated charged colloids in an electrolyte suspension is well-described by the DLVO theory [10]. This theory predicts screening of the electric charge on the colloids by counterions and effective interactions of Yukawa type, U(r) = U0exp(−κ⁎r)/r, where 1/κ⁎ is the effective screening length and U0 is the magnitude of the potential. In the present context we do not need to include the van der Waals term – which is also a part of the theory – because it is unimportant if only equilibrium properties of colloidal crystals are studied. Due to the strong electrostatic repulsion the mean separation between colloids is namely always much larger than the range of the van der Waals interaction. When there are more than two colloids in the system, i.e., when we analyze a dense colloidal suspension, pairwise additivity cannot be generally assumed due to the presence of many-body interactions. It has been observed in Ref. [11] and in Refs. [5,12] that the effective potential between colloids in a dense system deviates from the simple Yukawa form and that, due to the many-body interactions, it depends on the density and other state parameters of the system. In Refs. [5,12] the solid– liquid melting line in the three-dimensional colloidal suspension
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
Fig. 1. Construction of the lattices.
has been calculated by combining the solutions of the nonlinear Poisson–Boltzmann equation with Brownian dynamics simulations. The Poisson–Boltzmann theory correctly accounts for all many-body interactions. The resulting phase diagram has been compared to the phase diagram for a pairwise interacting Yukawa system [13] and systematic deviations have been observed. This effect of many-body interactions has been relatively successfully modeled by introducing a density-dependent cutoff of the effective Yukawa potentials such that the potential is simply truncated at a given distance dC. If dC is between the first and second neighbor shell, the truncated Yukawa potential only includes the nearest neighbor interactions. Using such effective pair interactions, the many-body phase diagram has been reproduced reasonably well. This was our motivation to study the phase diagram of 2D colloidal crystals interacting via the truncated Yukawa interaction. We have already fixed the unit of length to dNN. The parameter κ⁎ rescales to κD = κ⁎dNN. Additionally, the unit of energy is chosen to be the energy per particle in a hexagonal lattice with nearest neighbor interactions only, EH/N = 3U0exp (− κD). The energy per particle in a 2D colloidal crystal interacting via the truncated Yukawa interaction is eL ðjD ; dC Þ ¼
L 1 X expð−jD ðrij −1ÞÞ : 6 j; rij b d rij
175
unit cell are generated. For example, a hexagonal phase is described by a =ffiaH = 1 and Δx = 1/2, while for square phase it is qp ffiffiffiffiffiffiffiffiffiffiffi ffiffiffi 3=2c0:93 and Δx = 0. We have systematically a ¼ aS ¼ varied both a and Δx and we have evaluated the energy (Eq. (2)) for a wide range values of κD and dC in each lattice. At small values of κD, the screening length is large and the interaction in relatively long-ranged. When κD is large, the interaction is short-ranged and at very high values the hard-core limit is approached. We find the hexagonal phase to be energetically most favorable at large κD values for all values of the cutoff dC. The value of dC can be varied between 1 (which corresponds to nearest-neighbor interactions only) and ∞ (which represents the simple Yukawa limit). Anything in between is simply an interaction terminated at a certain distance. Again, an infinite cutoff reproduces the simple Yukawa system where the hexagonal lattice is the stable phase for all κD values [4]. However, at the cutoff value around 1, where only the nearest neighbor interactions are considered, we find other phases to be stable at small values of κD. The calculated phase diagram is shown in Fig. 2. We see two regions where the square phase has the lowest energy and many regions where a deformed hexagonal phase was found to have the smallest energy among the candidate lattices considered. The deformed hexagonal phases are denoted by the letter H followed by a number denoting the value of a, while Δx is equal to 0.5 for all H-phases. We have not studied the mechanical stability of the lowest energy lattices in the phase diagram, so it could well be that the phases shown here are in fact not stable and therefore not the real ground states of the system. What we consider the main message of the text is, however, that the hexagonal phase is destabilized by the cutoff in the effective interactions. One must be aware that due to the discrete structure of the lattice distances the phase diagram is not a smooth function of the parameter dC. If we scanned the phase diagram in more detail, it would exhibit a more elaborate structure, possibly
ð2Þ
C
The energy eL depends on the lattice type L and on two parameters, κD and dC. 3.1. Phase diagram The general crystal lattice in 2D with a single-particle per unit cell is described by three parameters. Once the density is fixed, the remaining two parameters are the angle between the base vectors and the ratio of their lengths. We have constructed the lattices as illustrated in Fig. 1. The parameters of a lattice are the particle distance along the row a and the horizontal displacement of two consecutive rows Δx. The vertical spacing pffiffiffiffiffiffiffiffiffiffi between the rows, Δy, is fixed by density, Dy ¼ 3=2a. By varying a and Δx all possible 2D lattices with one particle per
Fig. 2. Truncated Yukawa potential: Phase diagram. The symbols H followed by apnumber denote a deformed hexagonal lattice (see main text), CR represents the ffiffiffi 3R30- triangular commensurate lattice [14]. We see the regions of square and CR lattices having smaller energy than the hexagonal lattice.
176
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
including other distinct phases on top of those plotted in Fig. 2. However, it was not our intention to deliver a precise phase diagram for truncated Yukawa interaction as it depends on the exact shape of the cutoff which is to a large extent arbitrary although the existence of the cutoff itself is not. Given this reservation, the phase diagram includes all lattices generated by the described procedure, which have one particle per unit cell, and a few lattices with more ffithan one particles per unit cell, like pffiffiffiffiffiffiffiffiffiffiffiffi the CR structure ( 3R30- triangular commensurate lattice [14]). As a word of caution, we should stress that it is possible that some more complicated and exotic lattices with several particles per unit cell could also intervene and replace some of the phases that are present in the phase diagram. But even if this would be the case, what we find to be most important is that for truncated Yukawa potential, we do find phases other than hexagonal to be the ground state for a range of cutoffs. 4. Poisson–Boltzmann cellular model The energy of a particle in a two-dimensional colloidal crystal can be obtained by solving the nonlinear Poisson– Boltzmann equation within a unit cell with periodic boundary conditions at the cell boundaries such that the normal derivative of the potential at the boundary is zero. Let us denote the unit cell of lattice L by GL. GL is the volume of the unit cell outside the colloid, i.e., the solvent volume. The grand canonical Poisson–Boltzmann equation reads j2 w ¼ j2 sinhðwÞ
for
Y raGL ;
¼−
ZkB R2
for
Y raAGc ;
ð4Þ
where R is the radius of the colloid. The second boundary condition should account for the interaction of the central particle with all other particles in the crystal. It is not difficult to see that this is taken care of by demanding that the normal derivative of the potential at the cell boundary ∂GL be 0, which is nothing but the method of images applied in the reverse direction. Eq. (3) together with the two boundary conditions represents a complete description of the electrostatic problem. Once the potential ψ(r) is known, the grand potential energy can be computed as follows 1 ZkB 4pkB bX ¼ 2 R
Z AG Zc
wdr
þ ðjRÞ2
ðw sinhw−2 coshw þ 2Þdr GL
In order to systematically study the hexagonal and the square lattice and the lattices in between them, we partition the crystal into rectangular cells containing 2 rather than 1 particle. We denote the sides of the rectangle pffiffiffiffiffiffiffiffi by a and b, the latter being fixed by the density, b ¼ 3=a; all distances are measured in units of nearest neighbor distance in the hexagonal lattice. The positions of the two particles within the cell are denoted by r1 and rq2.ffiffiffiffiffiffiffiffiffiffiffiffi Let us put δ = Δx/a. If a = aS + 2δ(1 − aS) where pffiffiffiffiffiffiffiffiffi aS ¼ 3=2, then the shear mode
að1 þ dÞ=2 ; b=4 að1−dÞ=2 r2 ¼ 3b=4
r1 ¼
ð6Þ
continuously transforms the hexagonal phase (δ = 0.5) to the square phase (δ = 0); δ clearly plays the role of the order parameter. We have solved the Poisson–Boltzmann equation in the appropriate unit cell for several values of δ between 0 and 0.5 and we have evaluated the energy (Eq. (5)). We did this for a whole range of the relevant parameters characterizing the solution of the Poisson–Boltzmann equation (the parameters are ZλB, κa and the volume fraction of the colloids). The hexagonal phase always had the lowest energy and in contrast to the results from the truncated Yukawa interactions, we did not see any other phases to be stable.
ð3Þ
where ψ = eϕ/κBT is the reduced electrostatic potential and κ2 = 4πλBcS; here λB = e2/4πϵϵ0κBT is the Bjerrum length and cS is the reservoir monovalent salt concentration. The first boundary condition is determined by the fixed surface charge of the colloids. By denoting the surface of the colloid by ∂Gc, and its normal vector pointing into the solvent by Y nˆ p, this boundary condition can be written as Y nˆ p djw
4.1. Hexagonal to square lattice
ð5Þ
5. Stability Another point of view is to analyze the stability of the hexagonal lattice with respect to shear deformations. This is usually done by parameterizing the transition from the initial to the final crystal configuration by a suitably chosen order parameter, such as the Bain strain which transforms a BCC lattice into an FCC lattice via a continuous sequence of BCT lattices. Of course, there are infinitely many shear modes whose complexity grows with the size of the final configuration, which may contain more than 1 particle per unit cell. But for the simple lattices, say a square lattice, it is plausible to expect that the mode will also be simple and amenable to analysis. The parameterization of the shear modes is not unique; the only constraint they must satisfy is that the density is fixed. This program, used in Ref. [1], consists of a harmonic expansion of the total energy in terms of the order parameter ϵ starting with an expansion of the interparticle distances Ari f ri ðeÞ ¼ ri ðe0 Þ þ ðe−e0 Þ Ae e¼e0 1 A2 ri ðe−e0 Þ2 þ : : : ð7Þ þ 2 Ae2 e¼e0 where r i(ϵ) is a suitable parameterization of the structural transformation and ϵ0 corresponds to the configuration
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
under consideration. To second order in ϵ, the energy reads 1X Ari /ðri ð0ÞÞ þ / Vðri ð0ÞÞ ðe−e0 Þ U ðaÞ ¼ 2 i Ae e¼e0 2 1 A ri þ/ Vðri ð0ÞÞ ðe−e0 Þ2 2 Ae2 e¼e0 2 1 W Ari þ / ð r i ð 0Þ Þ ðe−e0 Þ2 2 Ae e¼e0
½
ð8Þ
where i runs over all neighbors of the central particle and ϕ is the pair potential; ϕ′ = dϕ/dr, etc. By assumption, the initial lattice corresponds to an extremum of the energy, so terms linear in ϵ cancel out and we have U ðeÞ ¼ U ðe0 Þ þ
ðe−e0 Þ2 a; 2
ð9Þ
where U(0) = (1/2)∑ i ϕ(r i) is the energy of the initial state and
A2 U au Ae2 ¼
1X 2
i
"
e¼e0
2 2 # A r Ari i / Vðri ðe0 ÞÞ þ/W ðri ðe0 ÞÞ : Ae2 e¼e0 Ae e¼e0 ð10Þ
If α b 0, the lattice with ϵ = ϵ0 is unstable with respect to the shear mode studied. We stress again that α depends on the parameterization of the shear mode, and so this analysis can but need not indicate that a certain lattice is unstable. A general analysis of the structure of the two terms in Eq. (10) is difficult, and we will limit our discussion to the square and the hexagonal lattice and the shear mode introduced in Section 4. Within this framework, the derivatives can be calculated analytically but the results themselves are not very
177
instructive. However, a numerical evaluation of the derivatives shows that for the two lattices in question, all shell subsums C1 = ∑i∈shell∂2ri/∂ϵ2 are positive. For repulsive potentials, ϕ′ b 0 for all r and the first sum in Eq. (10) will always be negative, and the lattice will only be stable if the second term in Eq. (10) is dominant. Shell subsums C1 and C2 = ∑i∈shell(∂ri/∂ϵ)2 for the hexagonal and the square lattice are shown in Fig. 3. The two lattices do have much in common but there are also some important differences. In general, the magnitudes of C1 and C2 are rather similar, and C1 and C2 both increase with separation when appropriately averaged such that the discreteness of the diagrams is washed out. However, the short-range behavior of the shell subsums in the hexagonal lattice is very different than in the square lattice. One notes immediately that in the hexagonal lattice, C2 is almost constant and rather large within the first six or so shells, whereas in the square lattice it goes to 0 with decreasing separation. This implies that the stability of the hexagonal lattice largely depends on the short-range behavior of the pair potential, and the sign of α (i.e., the stability of the lattice) can be estimated by considering the first shell alone. In the square lattice, many more shells must be taken into account to compute the exact value α and to see whether this lattice is stable. Another consequence of the large relative magnitude of C2/C1 in the first few shells of the hexagonal lattice (which is ≈5 when averaged over the first few shells) is that it can be easily destabilized by a shoulder-like, concave pair potential where ϕ″ b 0 at small separations. On the other hand, the square lattice would not be affected by such potentials so easily as in the first lattice, the magnitude of C2/C1 is about 1. This is, we believe, the mechanism that stabilizes the numerous non-close-packed lattices observed in numerical simulations of the so-called core-softened particles [6]. We stress that the same effect can also be achieved with potentials with a sharp cutoff where the energies of the various lattices differ primarily because of the different number of shells covered by the potential. Although limited in scope, this illustration does point out that the effective pair potential could render the hexagonal lattice
Fig. 3. Shell subsums for the hexagonal lattice (left) and the square lattice (right). As far as C1 is concerned, the two lattices are not very different, but within the first few shells C2 of the hexagonal lattice is much larger than in the square lattice. This gives rise to the robustness of the hexagonal lattice in case of convex repulsive pair interactions such as the Yukawa potential. Note that the scales on the top and bottom panels are not the same.
178
J. Dobnikar, P. Ziherl / Journal of Molecular Liquids 131–132 (2007) 173–178
unstable not only if it were characterized by an abrupt cutoff but also if it were concave at small separations. Whether this can happen in charged colloidal systems is presently unexplored: The Yukawa potential is strictly convex at all separations, and the many-body corrections to the pair potential should be rather strong to make the curvature of the effective potential negative. 6. Conclusions We have presented the T = 0 energy calculations for 2D crystal structures in systems interacting via the truncated Yukawa pair interaction. We have evaluated the energies for phases continuously varying from the hexagonal to the square lattice. The results show that the cutoff introduced with the intention to capture many-body effects destabilizes the hexagonal lattice at certain parameters values. The destabilization is also consistent with the additionally performed stability analysis in Section 5, where we show that the interaction potential needs to be concave in order to destabilize the hexagonal phase. Although we have only analyzed a single shear mode chosen for its simplicity, we believe that other possible ways of parameterization of the hexagonal-to-square transformation should produce rather similar results. However, we have found the hexagonal phase to be stable at all parameters in the full Poisson–Boltzmann description, which means that the truncated Yukawa description does not correctly capture many-body effects in the system under investigation, although previous studies in 3D [5,12] claim the opposite. This finding shows – once again – that the implementation of effective
interactions in the genuine many-body systems is subtle and should be applied with considerable precautions. Acknowledgments We kindly acknowledge useful discussions with H.H. von Grünberg and O. Glatter. This work was supported by a Marie Curie fellowship (MEIF-CT-2003-501789). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
E.A. Jagla, Phys. Rev., E 58 (1998) 1478. E.A. Jagla, J. Chem. Phys. 110 (1999) 451. G. Malescio, G. Pellicane, Phys. Rev., E 70 (2004) 021202. F.M. Peeters, X. Wu, Phys. Rev., A 35 (1987) 3109. J. Dobnikar, R. Rzehak, H.H.v. Grünberg, Europhys. Lett. 61 (2003) 695. K. Kremer, M.O. Robbins, G.S. Grest, Phys. Rev. Lett. 57 (1986) 2694. P.J. Camp, Phys. Rev., E 68 (2003) 061506. G. Malescio, G. Pellicane, Nat. Mater. 2 (2003) 97. A. Skibinsky, S.V. Buldyrev, A. Scala, S. Havlin, H.E. Stanley, Phys. Rev., E 60 (1999) 2664. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions, Cambridge University Press, 1989. M. Brunner, C. Bechinger, W. Strepp, V. Lobaskin, H.H.v. Grünberg, Europhys. Lett. 58 (2002) 926. J. Dobnikar, Y. Chen, R. Rzehak, H.H. von Grünberg, J. Chem. Phys. 119 (9) (2003) 4971. M.O. Robbins, K. Kremer, G.S. Grest, J. Chem. Phys. 88 (1988) 3286. S. Bleil, H.H. von Grünberg, J. Dobnikar, R. Castañeda-Priego, C. Bechinger, Europhys. Lett. 73 (3) (2006) 450.