Martial Mazars
email: martial.mazars@th.u-psud.fr

The Chiral Dipolar Hard Sphere Model

Keywords: Chirality, Monte-Carlo simulations, Liquid Crystals

The Chiral Dipolar Hard Sphere Model.

Introduction

Chirality is used to describe a lack of symmetry in compounds : it refers to objects, molecules or structures that are not superposable on their mirror image. The separation of the two optical active forms of the tartaric acid, based on the observation of the symmetry of crystals by Pasteur in 1847, has been a major step in the understanding of optical properties of crystals and liquids [START_REF] Pasteur | Études des Phénomènes Relatifs à la Polarisation Rotatoires des Liquides[END_REF]. Later, in 1884, Kelvin introduced the term chirality, but it has been extensively and systematically used only eighty years later with the introduction of the Cahn-Ingold-Prelog priority rules [2]. Since Arago, Biot and Pasteur, the ideas and concepts on optical activity of substances and on chirality have had and still have a lot of applications in chemistry, biology and physics.

In chemistry and biochemistry, chirality is an important tool in stereochemistry and for the determination of the functions of molecules. The interactions between biological molecules are frequently a consequence of the chiral interactions between them. For instance, as first noted by Pasteur, the two forms of asparagine (dextro and levo) tasted differently ; it was found later that this is a frequent property the two optical active forms of amino acids, because of the chiral structure of the taste bud receptor site (several others interesting examples may be found in [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]).

The helical structures of DNA, RNA, proteins, etc. are also important to permit biological recognition and interactions between biological structures [START_REF] Kornyshev | [END_REF]5].

In condensed matter physics many properties depend on the chirality of constituants or structures. More specifically, in liquid crystal physics, many phases as cholesteric, some smectics phases (S * A , S * C , etc.) and blue phases possess an optical rotatory power [6][7][8][START_REF] Chandrasekhar | Liquid Crystals[END_REF][START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Kamien | [END_REF][12] ; in most cases, a helical arrangement of the molecules is responsible for this property. Obviously, these properties are of a particular interest for technological applications. It is not necessary for the molecules to be chiral to observe chiral liquid crystal phases [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]13,14] ; however, some particular phases, such as the ferroelectric smectic-C * phase, may be observed only when some chiral molecules are present in the sample (see section 5.10 of ref. [START_REF] Chandrasekhar | Liquid Crystals[END_REF] and refs. [15][16][17]).

To describe the optical activity of liquids and solids, and the discriminatory interaction between chiral molecules several molecular models have been built. Among the most known are Kuhn models that are based on coupled-oscillators arranged on an asymmetric tetrahedron [18,[START_REF] Caldwell | The Theory of Optical Activity[END_REF], the optical activity tensor or gyration polarizability may be computed in these models [START_REF] Osipov | [END_REF][21][22]. Models based on chirality of helix have also been used [21,23,24]. These helical models are important for the study of the properties of carbon nanotubes [25][26][27] and biological macromolecules [START_REF] Kornyshev | [END_REF]28,29].

For the study of chirality in fluids and solids, the discriminatory interactions between enantiomers are important too [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]. The difference in the interaction between enantiomeric molecules may have several origins, it may stems from steric, electrostatic, magnetic or dispersion interactions [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]31,32]. These chirodiastaltic interactions [START_REF] Craig | 3] have chosen the term chirodiastaltic ('diastaltic' ≡ 'serving to distinguish') because of its explicit reference to chirality. Another terminology diastereotopic has been introduced[END_REF] are responsible for some differences in thermodynamical and structural properties between enantiopure systems (homochiral systems -systems that are composed with exclusively one enantiomer) and racemic systems (systems composed in equal proportions of the two enantiomers). For instance, the melting point of crystals of the enantiomeric pure (+)-tartaric acid is 170 o C, while, for the racemic crystals, the melting temperature is 204 -6 o C. Some differences are also reported for boiling points, but since the chirodiastaltic interactions are quite small in general, these difference in boiling points are small and do not occur systematically. In liquid crystals, the phase diagrams depend also on the fractional concentrations of enantiomers [33].

Several simple molecular models have been proposed and used to study the thermodynamical and structural properties of fluids and solids of chiral molecules in computer simulations and integral equation theories. In these simple molecular models, the chirodiastaltic interactions are represented either via steric interactions [34][35][36][37][38] or electrostatic interactions [39][40][41][42] ; a chiral Gay-Berne model has also been built [43]. In the present work, we describe another simple model of chiral molecules, the chiral dipolar hard sphere model. In this model, the chirodiastaltic interactions is represented by electrostatic (or magnetic) dipoles-dipoles interactions, it consists of hard sphere with two point dipoles located inside the molecule. As shown by Craig and coworkers, electrostatic interactions between enantiomers may be discriminatory at most if the dipole-quadrupole interaction between two enantiomers differs [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]44,45] ; this implies that if the quadrupole moment of the two form of enantiomers differs then the interaction may be discriminatory. Therefore, a chiral model based on electrostatic interactions must have at least four point electric charges arranged on an asymmetric tetrahedron geometry, as the models in refs. [39][40][41][42], two non-colinear point dipoles [START_REF] Camp | An example of an achiral model of dipolar hard spheres possessing two colinear dipoles is given in the reference[END_REF] or one quadrupole with no planar symmetry [44]. The present paper is organised as follow. In section 2, we describe the model of chiral dipolar hard sphere in detail and give some of its elementary properties. More precisely, in section 2.1, we define the model, the nomenclature, its pertinent parameters and we compute all its multipole moments ; it is shown that the knowledge of only three multipole spherical components (namely : Q determination of its electrical properties and also, in view of the determination of the optical activity of the model, for the computation of multipole-multipole polarizabilities [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF][START_REF] Caldwell | The Theory of Optical Activity[END_REF]44,45,[47][48][49][50]. In section 2.2, we compute interaction energies between two enantiomers in some particular configurations ; it is shown that, despite the simplicity of the molecular model, the energy landscape of the interaction between two enantiomers is quite complicated. In section 3, we report some preliminary results of Monte Carlo simulations performed in the canonical (N V T ) and isobaric (N P T ) ensembles ; in these computations, all dipole-dipole interactions between molecules are explicitly taken into account by using the Ewald method for dipolar interactions [START_REF] Weis | Advanced Computer Simulation Approaches for Soft Matter Sciences II[END_REF]. The paper ends with a discussion on the perspectives that this model offers to the study of chiral fluids and solids. For completeness, in appendix A, we give general formulas for the computation of multipole moments for any dipole distribution.

Last, but not least, this paper is dedicaded to my colleague Jean-Jacques Weis.

2.

Elementary properties of chiral dipolar hard sphere model.

In this section, we describe the chiral dipolar hard sphere model and we give some of its elementary properties. The subsection 2.1 is devoted to the definition of the model and to the computation of the multipole moments due to the two permanent dipoles moments ; all multipole moments are computed with the center of the sphere taken as origin.

In subsection 2.2, we compute interaction energies between two enantiomers for few particular configurations ; in this subsection, a particular attention is pay to the chirodiastaltic interactions for colinear aligned total dipole moments and antiparallel dipole moment configurations. The discriminatory interaction energies are often computed by using the multipole expansion [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]22,24,39,41,44,47,49,50] ; thus, to facilitate the comparison with these works, the discriminatory interaction energies between chiral dipolar hard sphere have also been computed analytically by using the multipole expansion when the distance between both enantiomers is large in comparison to their diameter. However, if both enantiomers are at contact or if the distance between them is too small to use with accuracy the multipole expansion, then the interaction energies must be computed by summing all dipole-dipole interactions ; it is done so in subsection 2.2 for all configurations at contact and in section in all Monte-Carlo computations done with the Ewald method.

2.1. Definition of the model and its multipole moments.

The Chiral Dipolar Hard Sphere model consist of hard spheres of diameter σ with two point dipoles µ 1 and µ 2 located at points O 1 and O 2 inside the sphere. We define the molecular axis û of a chiral dipolar hard sphere as O 1 O 2 = 2Lû with 2L < σ. The dipoles µ 1 and µ 2 are chosen such that µ a .û = 0 and we set

   µ 1 = µ μ1 µ 2 = λµ (cos α μ1 + sin α ( μ1 × û)) (1)
where the parameters α ∈ [-π, π] and 0 ≤ λ ≤ 1 define the second dipole µ 2 from the dipole µ 1 and the axis û. We name Rectus (R-) bidipolar hard sphere (R α -HS) for 0 < α < π and Sinister (S-) bidipolar hard sphere (S α -HS) for -π < α < 0 (see figure 1 ). In this nomenclature, it is clear that R (-α) -HS ≡ S α -HS ; molecules R α -HS and S α -HS are enantiomers. Chiral parameters or chiral indices have been introduced to give a quantitative measure of the chirality of a molecule [START_REF] Osipov | [END_REF][START_REF] Harris | [END_REF]. In particular, A.B. Harris and co-workers [START_REF] Harris | [END_REF] have defined chiral parameters as pseudoscalars constructed from the structure of molecules described via the representation theory of the threedimensional rotation group O(3) ; they have also related chiral parameters to the cholesteric pitch of generic liquid crystal models. For the chiral dipolar hard sphere model, we may defined a pseudoscalar from the three vectors û, µ 1 and µ 2 as

χ = Lû.(µ 2 × µ 1 ) = Lλµ 2 sin α (2)
and take χ as the chiral parameter of this model ; we have χ(S α -HS) = -χ(R α -HS). If χ = 0, then both enantiomers are equivalent and the bidipolar hard sphere model is achiral. This occurs if one of the following conditions is verified : (a) L = 0, then both point dipoles are located at the center of the sphere and the model reduces to the dipolar hard sphere model ; (b) λ = 0, in this case, there is only one point dipole in the molecule and (c) sin α = 0 then the dipoles and the molecular axis are all in the same plane which is a plane of symmetry for the molecule.

A physical interpretation to | χ * |= 1 is less obvious ; to obtain | χ * |= 1 one needs to have λ = 1, thus a symmetry is restored and then the chirality is perhaps less marked for these particular cases with λ = 1 (for instance, see below subsection 2.2 and Fig. 7). A trivial interpretation of values | χ * |= 1 would be to say that for these value a maximum discrimination is achieved ; however, this is perhaps not fully correct.

We define the molecule-fixed local frame as 2)). Enantiomers Rα-HS and Sα-HS have equal spherical components when an even symmetry occurs, that is to say for P and Q 21 . One may note that for λ = 1 one has Q 21 = 0. where P is the total dipole of the molecule. In this molecule-fixed frame, the projections of the dipoles onto the frame axis are given by It is worthwhile to note that the total dipole of enantiomers R α -HS and S α -HS in the molecule-fixed frame are the same, we have P (R α -HS) ≡ P (S α -HS). One should note also that P = 0 can be obtained only if λ = 1 and α = π, then χ = 0 and the model is achiral ; thus, all chiral bi-dipolar hard sphere have a non-zero net total dipole. The dipole distribution of the chiral hard sphere may be written as

             êz = µ 1 + µ 2 | µ 1 + µ 2 | = P | P | êx = û êy = êz × û (3)
           µ 1 = µ(1 + 2λ cos α + λ 2 ) -1/2 [-λ sin α êy + (1 + λ cos α) êz ] µ 2 = µ(1 + 2λ cos α + λ 2 ) -1/2 [λ sin α êy + λ(cos α + λ) êz ] P = µ(1 + 2λ cos α + λ 2 ) 1/2 êz
µ(r) = µ 1 δ(r + Lû) + µ 2 δ(r -Lû) (5)
and from multipole expansion for any dipolar distribution (see appendix A), the cartesian components of the quadrupole moment are given by

[q] = 3 2 Lµ √ 1 + 2λ cos α + λ 2   0 2λ sin α (λ 2 -1) 2λ sin α 0 0 (λ 2 -1) 0 0   (6)
and spherical components by

                   Q 22 = i 15 2π Lλµ sin α √ 1 + 2λ cos α + λ 2 Q 21 = 15 8π Lµ(1 -λ 2 ) √ 1 + 2λ cos α + λ 2 Q 20 = 0 (7)
For enantiomers, we have

Q 22 (R α -HS) ≡ -Q 22 (S α -HS) and Q 21 (R α -HS) ≡ Q 21 (S α -HS). The magnitude of the quadrupole is [53] Q2 2 = 4π 5 m | Q 2m | 2 = 3 2 L 2 µ 2 4λ 2 sin 2 α + (1 -λ 2 ) 2 1 + 2λ cos α + λ 2 (8)
Obviously, one has Q2 2 (R α -HS) = Q2 2 (S α -HS). In Eq.( 8), it might be surprising to find Q2 2 = 0 for λ = 0, since there is only a single dipole in the hard sphere in these cases. This stems from dependence on the choice of the origin in the computation of the quadrupole (the single dipole is shifted by -Lû from the center of the hard sphere -see also ref. [START_REF] Gray | Theory of molecular fluids[END_REF], p.69). This origin dependence highlight the fact that when setting λ = 0, the chiral bidipolar hard sphere model do not reduce to the classical dipolar hard sphere model. On Figs. (2)[START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF], we show the dipole and quadrupole components as functions of the chiral parameter χ and angle α ; on these representations, the influence of chirality is related to odd or even properties of the function. As shown on Fig. 23, a difference between enantiomers R α and S α is found only in Q 22 . The spherical components of the octopole are given by

                           Q 33 = 0 Q 32 = 1 2 105 8π L 2 µ(1 + 2λ cos α + λ 2 ) 1/2 = 1 2 35 2 L 2 Q 10 Q 31 = 0 Q 30 = - 3 2 7 4π L 2 µ(1 + 2λ cos α + λ 2 ) 1/2 = - √ 21 2 L 2 Q 10 (9)
and its magnitude is

Q2 3 = 4π 7 m | Q 3m | 2 = 33 8 L 4 µ 2 (1 + 2λ cos α + λ 2 ) = 11π 2 L 4 Q 2 10 ( 10
)
It is apparent also on Eqs.(9-10) that the octopole moments are the same for both enantiomers. More generally, by using Eq.( 5) and the definition of the spherical component of multipoles Eq.(A4), we found

Q ll µL (l-1) = il λ sin α √ 1 + 2λ cos α + λ 2 (1 + (-1) l )Y ll (ê x ) (11)
and, for 0 ≤ m ≤ l -1, we have

Q lm µL (l-1) = -(l + m + 1)(l -m) (λ 2 -(-1) l ) + (1 -(-1) l )λ cos α √ 1 + 2λ cos α + λ 2 Y l(m+1) (ê x ) +im λ sin α √ 1 + 2λ cos α + λ 2 (1 + (-1) l )Y lm (ê x ) (12) The value of Y lm (ê x ) may be easily computed, it comes Y lm (ê x ) = 2l + 1 4π (-1) (l+m)/2 [(l -m)!(l + m)!] 1/2 (l -m)!!(l + m)!! for (l + m) even 0 for (l + m) odd (13)
therefore, only one of the two contributions in the right handed side of Eq.( 12) gives Q lm . From this computation, we see that all spherical components can be computed only from Q 10 , Q 21 and Q 22 ; more precisely, from Eqs. [START_REF] Kamien | [END_REF](12), we have

                                     Q ll µL (l-1) = l(1 + (-1) l ) 2π 15 Q 22 µL Y ll (ê x ) Q lm µL (l-1) = (1 + (-1) l ) 8π 15 (l + m + 1)(l -m) 2 Q 21 µL Y l(m+1) (ê x ) +m Q 22 µL Y lm (ê x ) - (1 -(-1) l ) 2 (l + m + 1)(l -m) 4π 3
Q 10 µ Y l(m+1) (ê x ) (14)
Since R (-α) -HS ≡ S α -HS, we have

Q lm (R α ) -Q lm (S α ) = 2(1 + (-1) l )mL (l-2) 2π 15 Q 22 Y lm (ê x ) (15)
thus, all spherical components of multipoles of enantiomers R α -HS and S α -HS are equal, but those with l and m even and m = 0. It is worthwhile to note that Q lm is proportional to only one of the three components Q 10 , Q 21 or Q 22 , because of Eq.( 13). Eq.( 14) shows that the behavior of all components of multipoles are given (up to a numerical multiplicative factor) by Figs. (2 -3) as functions of the model parameters ; in particular the difference As shown in Eq.( 14) all components of the multipoles dependent only on the three components Q 10 , Q 21 or Q 22 , therefore it is advantageous to define the three reduced geometrical parameters χ * , P * and r * as

Q lm (R α ) -Q lm (S α ) is
                   χ * = χ Lµ 2 = λ sin α P * = P µ = (1 + 2λ cos α + λ 2 ) 1/2 r * = r Lµ = (1 -λ 2 ) (16)
2.2. Interaction energies between two enantiomers.

The interaction energy between two enantiomers can be computed in two ways. First, by summing up all interaction energies due to dipole-dipole interactions and/or, second, by using the multipole expansion. In the next section, we give some preliminary numerical results obtained with Monte Carlo Metropolis sampling of the phase space ; in these computations, all dipole-dipole interactions are taken into account explicitely and long ranged interactions are computed with the Ewald method [START_REF] Weis | Advanced Computer Simulation Approaches for Soft Matter Sciences II[END_REF][START_REF] De Leeuw | Proc. R. Soc. Lond. A[END_REF].

In the present subsection, we give some analytical results for the interaction energy between two enantiomers by using the multipole expansion when r σ > 2L and we give also some analytical and numerical results obtained by computing explicitly all dipole-dipole interactions between two enantiomers at contact (i.e. From the multipole expansion, the interaction energy between the multipoles of orders l 1 and l 2 for two molecules of general shape is given by [START_REF] Gray | Theory of molecular fluids[END_REF]

U l1l2 = A l1l2 n1,n2 Q l1n1 Q l2n2 r l+1 m1,m2,m C(l 1 l 2 l; m 1 m 2 m) ×D l1 m1n1 (ω 1 ) * D l2 m2n2 (ω 2 ) * Y lm (ω) * (17)
where l = l 1 +l 2 and A l1l2 is a numerical constant, Q ln the components of multipoles in the body-fixed frame, C(l 1 l 2 l; m 1 m 2 m) the Clebsch-Gordan coefficients, D l mn are the rotation matrices, the Euler angles ω a ≡ (φ a θ a χ a ) define the orientation of molecule a and ω ≡ r.

The first contribution to interaction energy is given by the total dipole-total dipole interaction and it is given by

U 11 (a ≡ X α , b ≡ Y α ) = 1 r 3 [P a P b -3(P a .r)(P b .r)] (18)
where X α , Y α ≡ R α or S α . This contribution is not discriminatory between enantiomers (see Eq.( 4)) and it is the longest ranged interaction between two chiral dipolar hard spheres. Therefore, many of the properties of this model will coincide with the properties of dipolar hard sphere systems that have been extensively studied by Jean-Jacques Weis in the last decade [START_REF] Weis | Advanced Computer Simulation Approaches for Soft Matter Sciences II[END_REF][START_REF] Weis | [END_REF][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72].

The discriminatory interaction between enantiomers (Chirodiastaltic interaction [START_REF] Craig | 3] have chosen the term chirodiastaltic ('diastaltic' ≡ 'serving to distinguish') because of its explicit reference to chirality. Another terminology diastereotopic has been introduced[END_REF]), may be defined by the difference ∆U

(LU ) = U (R α , R α ) -U (R α , S α ) (
where LU stands for Like -U nlike), we note also ∆U

(LU ) l1l2 = U l1l2 (R α , R α ) -U l1l2 (R α , S α
) as the contribution of multipoles of orders l 1 and l 2 to the discriminatory chirodiastaltic interaction. In the following, we choose the molecule at origin to be an R α enantiomer and all coordinates of vectors and positions are given in its moleculefixed local frame given by Eq.(3) -this is equivalent to choose ω 1 ≡ (000) in Eq. (17). From equation (15), we have ∆U (LU ) l1l2 = 0 only if l 2 even, thus, the first non zero contribution to ∆U (LU ) is given by the dipole-quadrupole interaction (l 1 = 1, l 2 = 2) and after some algebra, it comes

∆U (LU ) 12 = Lµ 2 r 4 λ sin α [15 M 2 (φ 2 θ 2 χ 2 ; θφ) + 12 M 1 (φ 2 θ 2 χ 2 ; θφ) + 9 M 0 (θ 2 χ 2 ; θ)] (19) with                            M 2 (φ 2 θ 2 χ 2 ; θφ) = sin 2 θ cos θ (1 + cos 2 θ 2 ) sin 2χ 2 cos 2(φ 2 -φ) +2 cos θ 2 cos 2χ 2 sin 2(φ 2 -φ)] M 1 (φ 2 θ 2 χ 2 ; θφ) = sin θ(5 cos 2 θ -1) [cos 2χ 2 sin(φ 2 -φ) + cos θ 2 sin 2χ 2 cos(φ 2 -φ)] M 0 (θ 2 χ 2 ; θ) = (5 cos 3 θ -3 cos θ) sin 2 θ 2 sin 2χ 2 (20)
where ω 2 ≡ (φ 2 θ 2 χ 2 ) and ω = (φθ) are respectively the orientations of the second enantiomer and of the bond vector r in the molecule-fixed local frame of the first enantiomer. The quadrupole-quadrupole interaction contributes also to ∆U (LU ) ; ∆U (LU ) 22 may be written as

∆U (LU ) 22 = - 3 2
L 2 µ 2 r 5 f 22 (λ, α) A 0 + 20A 1 + 30 √ 2A 2 + 140A 3 + 35A 4 -6 L 2 µ 2 r 5 f 21 (λ, α) B 0 + 5 √ 2 B 1 -5 √ 2B 2 - 35 4 B 3 (21)
with the functions f 22 and f 21 describing the dependence of ∆U

(LU ) 22
on the model parameters ; these two functions are given by

           f 22 (λ, α) = λ 2 sin 2 α (1 + 2λ cos α + λ 2 ) = χ * 2 P * 2 f 21 (λ, α) = (1 -λ 2 )λ sin α (1 + 2λ cos α + λ 2 ) = r * χ * P * 2 (22)
In Figs. 4 (a,b), we represent these two functions versus the chiral parameter χ and, in Figs. 4 (c,d) as functions of α. The functions A and B depend on orientations ω 2 and ω. More precisely, we have

                           A 0 = (35 cos 4 θ -30 cos 2 θ + 3)[(1 + cos 2 θ 2 ) sin 2χ 2 sin 2φ 2 -2 cos θ 2 cos 2χ 2 cos 2φ 2 ] A 1 = sin θ (7 cos 3 θ -3 cos θ) sin θ 2 [cos 2χ 2 cos(φ 2 + φ) -cos θ 2 sin 2χ 2 sin(φ 2 + φ)] A 2 = sin 2 θ (7 cos 2 θ -1) sin 2φ sin 2 θ 2 sin 2χ 2 A 3 = sin 3 θ cos θ sin θ 2 [cos 2χ 2 cos(φ 2 -3φ) -cos θ 2 sin 2χ 2 sin(φ 2 -3φ)] A 4 = sin 4 θ[2 cos θ 2 cos 2χ 2 cos(2φ 2 -4φ) -(1 + cos 2 θ 2 ) sin 2χ 2 sin(2φ 2 -4φ)] (23) and
                               B 0 = (35 cos 4 θ -30 cos 2 θ + 3) sin θ 2 [cos 2χ 2 cos φ 2 + cos θ 2 cos 2χ 2 cos φ 2 ] B 1 = sin θ (7 cos 3 θ -3 cos θ) (1 + cos 2 θ 2 ) sin 2χ 2 sin(2φ 2 -φ) +2 cos θ 2 cos 2χ 2 cos(2φ 2 -φ) - 3 2 √ 2 sin 2 θ 2 sin 2χ 2 cos φ] B 2 = sin 2 θ (7 cos 2 θ -1) sin θ 2 [cos 2χ 2 sin(φ 2 -2φ) + cos θ 2 sin 2χ 2 sin(φ 2 -2φ)] B 3 = sin 3 θ cos θ[(1 + cos 2 θ 2 ) sin 2χ 2 sin(2φ 2 -3φ) -2 cos θ 2 cos 2χ 2 cos(2φ 2 -3φ)] (24
) Obviously, if the distance between two enantiomers is too small the multipole expansion ceases to be accurate, however, the interaction energy between enantiomers may still be computed by summing all dipole-dipole interactions. The total interaction energy between two enantiomers a and b has thus to be written as where E (ia,jb) denotes the interaction energy between the dipole i of enantiomer a with dipole j of enantiomer b ; it is given by

E (a,b) = E (1a,1b) + E (1a,2b) + E (2a,1b) + E (2a,2b) (25)
E (ia,jb) = 1 | r ij ab | 3 µ (a) i .µ (b) j -3(µ (a) i .r ij ab )(µ (b) j .r ij ab ) (26)
and the distance vectors between point dipoles are

                   r 11 ab = r + L(û a -ûb ) r 12 ab = r + L(û a + ûb ) r 21 ab = r -L(û a + ûb ) r 22 ab = r -L(û a -ûb ) (27)
Even for simple relative configurations, as the ones shown on Fig. 5(a-b), the total interaction energy between two enantiomers is a bit complicated when computed by summing all dipole-dipole interactions. For instance, the interaction energy between two R α enantiomers in a configuration given by Fig. 5(a) is

E (a,b) = µ 2 | z | 3 × 1 (1 + 2 L 2 z 2 (1 -cos φ)) 5/2 P * 2 + r * 2 P * 2 L z 2 (1 -cos φ) -1 -6 L z χ * sin φ + 2 χ * 2 P * 2 cos φ - L z 2 (2 -2 cos φ + sin 2 φ) + 1 (1 + 2 L 2 z 2 (1 + cos φ)) 5/2 P * 2 - r * 2 P * 2 L z 2 (1 + cos φ) -1 +6 L z χ * sin φ -2 χ * 2 P * 2 cos φ - L z 2 (2 + 2 cos φ + sin 2 φ) (28
) where cos φ = ûa .û b and (χ * , P * , r * ) given by Eq.( 16). The chirodiastaltic interaction between enantiomers R α and S α , computed by summing the four dipole-dipole interaction energies, for the configuration given on Fig. 5(a) is For | z | L, a Taylor expansion of Eq.( 29) gives

∆E (LU ) = E (Rα,Rα) -E (Rα,Sα) = µ 2 | z | 3 × 1 (1 + 2 L 2 z 2 (1 -cos φ)) 5/2 4 χ * 2 P * 2 cos φ -6 L z χ * sin φ -4 L 2 z 2 χ * 2 P * 2 (sin 2 φ + 2 -2 cos φ) + 1 (1 + 2 L 2 z 2 (1 + cos φ)) 5/2 -4 χ * 2 P * 2 cos φ + 6 L z χ * sin φ -4 L 2 z 2 χ * 2 P * 2 (sin 2 φ + 2 + 2 cos φ) (29)
∆E (LU ) = µ 2 | z | 3 24 χ * 2 P * 2 L 2 z 2 sin 2φ -60χ * L 3 z 3 sin φ + o L 4 z 4 = ∆U (LU ) 22 (θ 2 = 0, θ = 0) + ∆U (LU ) 32 (θ 2 = 0, θ = 0) + o L 4 µ 2 z 7 (30)
where we have set φ = φ 2 + χ 2 in Eqs. (23)(24) to obtain ∆U (LU ) 22

. Substracting Eq.( 30) from Eq.( 29), may allow to check the degree of accuracy of the multipole expansion for the configuration of Fig. 5 (a). On Figs.6, we represent interaction energies for enantiomers in the colinear aligned total dipole moment configuration (Fig. 5-a). These energies have been computed by summing all dipole interactions between the two enantiomers. To allow simple comparison between interaction energies, for different values of α, we have substracted to all curves the average, over the orientation ûb , of the interaction energy. For all λ, a minimum in the interaction energy between two different enantiomers is always obtained for φ = 0 (û a and ûb parallel) ; the anti-parallel configuration of molecular axis ûa and ûb (φ = π) is a relative minimum for the interaction energy for all λ < 1 (see Fig. 6(b-d) -if λ = 1, the configurations for φ = 0 and φ = π have the same energy) ; one should note also that these minima are relative, because the configurations like in Fig. 5 (a) do not permit to reach the smallest distance between two point dipoles. The minimum interaction energy between two enantiomers of same kind is not obtained for parallel or anti-parallel configurations of the molecular axis. The value φ m of the angle for which a minimum interaction energy is obtained, depends on parameters λ and α (see Fig. 6(a-c)), several values are given in Table 1 ; average values of < E XY > φ are also reported in this Table .  Therefore, in homochiral systems (systems that are composed with exclusively one enantiomer) with a large enough total dipole moment and in thermodynamical conditions where dipolar systems form chains [START_REF] Weis | [END_REF]58], we may observe formation of chains too, and the molecular axis of the enantiomers may have helical arrangement. Addition of chiral enantiomers in the system will reduce the range of the helical arrangement in the chains. For λ = 1, the minimum at φ = φ m has the same value than the one found at φ = φ m -π, therefore right and left-handed helical arrangement of the molecular axis in the chains will compete in homochiral systems. For the configuration given in Fig. 5(b) and for r L (θ = π/2 and θ 2 = π), the multipole-multipole chirodiastaltic interaction is found as

∆E (LU ) = -3 µ 2 r 3 4 L r χ * sin(φ 2 -φ -2χ 2 ) + L 2 r 2 χ * 2 P * 2 (3 cos 2(φ 2 -χ 2 ) -35 cos 2(φ 2 -χ 2 -2φ)) + o L 3 µ 2 z 6 = ∆U (LU ) 12 (θ 2 = π, θ = π 2 ) + ∆U (LU ) 22 (θ 2 = π, θ = π 2 ) + o L 3 µ 2 z 6 (
31) On Figs. 7 and8, we represent interaction energies for enantiomers in the antiparallel total dipole configurations (Fig. 5-b) for φ = 0 and π/4. On Figs. 7 (a,c) and 8 (a,c), we represent interactions energies between two R α enantiomers for λ = 1 and 0.75 and on Figs. 7 (b,d) and 8 (b,d) we represent interactions energies between enantiomers R α and S α . The coupling between φ and φ 2 renders quite complicated the energy landscape of antiparallel configurations. The relative minima for antiparallel configurations (Fig. 5-b) are smaller than the minima found for colinear configurations (Fig. 5-a), therefore the colinear and antiparallel configurations are in drastic competition by comparison to dipolar hard sphere systems. For two identical enantiomers, when φ = nπ/2 (with n = -1, 0, 1 or 2), the location of minimum interaction energies are independent on α. For instance, if φ = 0, the relative minima for two R α enantiomers are always obtained for φ 2 = 0 or π ; if λ = 1, these two minima have the same value, but if λ < 1 the two minima are not equivalent and their value dependent on α. On Fig. 9, we give the value of the minimum energy for φ = 0 and φ 2 = 0 or π for λ = 0.75 and 0.5 ; as shown on this figure for values of α > α m (λ) the configuration φ 2 = π is more stable than the configuration φ 2 = 0. The amplitude of variation of interaction energies with φ 2 are larger for antiparallel configurations than variation with φ for aligned colinear configurations. The locations of relative minimum interaction energies between enantiomers R α and S α depend not only on φ and φ 2 , but also on α. For instance on Fig. 7(b) (φ = 0) and for α = 5π/12 (75 o ), we found four different minima located at φ 2 ±141 o and ±39 o . The computations of interaction energies between two chiral dipolar hard spheres show that a lot of different configurations may be in competition, even for simple configurations as the ones shown on Figs.5. In the next section, we give some preliminary results obtained with Monte Carlo simulations for several cases.

Preliminary Monte Carlo results

As shown in the previous section, interaction energies between two enantiomers as function of their relative orientation may have a lot of relative minima ; therefore, despite the apparent simplicity of the model, the thermodynamical properties chiral dipolar hard sphere systems are quite complicated. In particular, because the energy landscape of the interaction between two enantiomers shows large variations, one may predict that frustration phenomena will play an important role in phases at low densities, low temperatures and also in solid phases. We define the reduced quantities as follow : for volume of the simulation box V * = V /σ 3 , density ρ * = ρσ 3 = N σ 3 /V , pressure p * = pσ 2 /kT and dipole moment µ * = (µ/kT σ 3 ) 1/2 . For a given model of enantiomer, defined by (λ, α) and L, the reduced total dipole moment P * of the molecule is related to µ * according to Eq.( 4). In all computations, we have chosen L * = L/σ = 1/4 as in the numerical computations done in the previous section (see . At the end of this section, a short discussion on the influence of L on structures is presented ; this influence is illustrated on two cases : L * = 0.025 and L * = 0.425. For notational convenience the asterisks will be drop in the following. Mixtures of enantiomers R α -HS and S α -HS are also of interest, we define densities ρ R = N R /V and ρ S = N S /V respectively for enantiomers R α -HS and S α -HS and the ratio

X = ρ S /ρ R = N S /N R .
The ratio X is related to the enantiomeric excess (ee) in enantiomer R α -HS by ee = 1 -X. In the following, we have considered only two mixtures of enantiomers : homochiral R α systems (X = 0) and racemic R α -S α systems (X = 1). Monte Carlo simulations have been performed in canonical (N V T ) and isobaric (N P T ) ensembles with system size N = N R + N S = 512 or 1000. Periodic boundary conditions are used and dipolar interaction energies between two enantiomers is computed with Eqs. (25)(26)(27) where long ranged contributions are taken into account by using the Ewald method [START_REF] Weis | Advanced Computer Simulation Approaches for Soft Matter Sciences II[END_REF].

From a technical point of view, as a consequence of the numerous number of relative minima in interaction energies, the convergence of the Monte Carlo algorithm and the sampling of the phase space will require a large number of MC cycles (one MC cycle corresponds, on average, to one trial move per particle for N V T simulations ; for N P T simulation, an additional trial move of the volme of the box is done every MC-cycle). For most results presented in this section, between 5 × 10 4 and 10 5 MC cycles have been achieved for the relaxation and convergence from an initial condition and averages have been accumulated over about 10 5 and 2 × 10 5 MC cycles, respectively for systems with 1000 and 512 particles. Periodic boundary conditions are used and interaction energies are computed with the Ewald method [START_REF] Weis | Advanced Computer Simulation Approaches for Soft Matter Sciences II[END_REF] by summing all dipole-dipole interactions between enantiomers. The average energy U is computed as

U = a =b E (a,b) (32)
where E (a,b) is the interaction energy between enantiomers a and b Eqs. (25)(26)(27) and < . > is the average of the MC sampling. Some computation have been done with an external electric field E 0 = E 0 êz , the energy in the external field is computed as a) .E 0 = -

V = - N a=1 P (

Molecular Physics Mazars˙TMPH˙2009˙0027

Table 2. Average energies and eigenvalues S+, S0 and S-for homochiral-Rα and racemic systems in external electric field E = Eêz (E = 10) obtained in canonical (N V T ) Monte Carlo simulations. The numbers in brackets give the accuracy on the last digit of the averages. For all systems N = 512, µ = 1.0 and ρ = 0.05 (or ρ = 0.5). The parameters α and λ define the molecular structure of enantiomers as given by Eq.( 1). a For these computations, the density is ρ = 0.5. Possible orientational order can be established by computing the eigenvalues of the tensor where xi α is the cartesian component of a unit vector x of molecule i. The largest eigenvalue S + of Q αβ is used as an order parameter for the isotropic-nematic phase transition in liquid crystals. It is also worthwhile to note that if, all unit vectors x are perpendicular to a given fixed direction then, the eigenvalues are S + = S 0 = 1/4 and S -= -1/2. In the present work, we have restricted the computation of Q αβ the total dipole moment : (O, P ) ; the second set is the site O 1 and the orientation of the dipole µ 1 : (O 1 , μ1 ). To study the influence of the chirodiastaltic interaction on structure of mixtures of both enantiomers, it will be necessary to compute projections of the pair distribution functions on rotational invariants of higher order since the discriminatory interaction involves dipole-quadrupole interactions. Table 2, summarizes average energies and eigenvalues of Q αβ Eq.( 34) for homochiral-R α and racemic systems in external electric field E = E ẑ (E = 10), for several values of α and λ. On Table 3, we report some preliminary results of Monte-Carlo computations in the isobaric (N P T ) ensemble of homochiral-R α systems and on Fig. 10, we show g 000 (r) obtained with these computations. On Fig. 11, we give few snapshots of homochiral-R α and racemic systems in an external field. The results reported on table 2 show that energies and the nematic order parameters for the molecular axis are almost independent on the composition of systems. The eigenvalues S + , S 0 and S -show that on average the molecular axis of enantiomers are perpendicular to the direction of the electric field, this results is quite obvious since the total dipole moment, as defined in Eqs. [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF][START_REF] Kornyshev | [END_REF], is perpendicular to the molecular axis this tendency decreases as α increases, since the total dipole moment decreases. As shown on the snapshots in Fig. [START_REF] Kamien | [END_REF], if α 60 o molecules tend to condense in a large columnar cluster, for both homochiral-R α and racemic systems in an external electric field (cf. Fig. 11 (a-d)) ; if α 70 o the columnar cluster does not form (cf. Fig. 11 (e-f)). There are at least two facts that may help to understand qualitatively the disappearance of the columnar clusters as α increases, while λ and µ are fixed. First, as α increases, the total dipole of each molecule decreases and thus formations of chains in the direction of the external field is less favoured. Second, for large enough α, configurations with two dipoles belonging to two different molecules may be more easily found, on average, in an antiparallel configuration in a plane perpendicular to the field, while the total dipole of both molecules can still be parallel to the field. At present, it is not clear how the columnar clusters form and which mecanism favour their formation ; further investigations are needed to be conclusive. Homochiral-R α and racemic R α -S α systems in an external electric field exhibit the same properties, however, a close inspection of the snapshots of Fig. 11 (a-d) shows that the internal structure of the columnar clusters is not the same in homochiral and racemic systems. In homochiral systems, the molecules in the columnar cluster have an helical arrangement (in agreement with the results shown on Fig. 6) ; in racemic systems, these helical arrangement in the columnar clusters are less marked because of the Like-Unlike interactions. To obtain more quantitative descriptions of the internal structure of the columnar clusters, some particular order parameters and distribution functions have to be defined. If there is no external electric field, columnar clusters are not formed instead some complicated structures are formed at low density and large enough µ ; this has occured for all values of α that have been considered in these prelyminary results. On Fig. 12, snapshots of homochiral-R α and racemic R α -S α systems are represented, for these systems α = 120 o and no external electric fileds are applied. Again, in these computations homochiral and racemic systems behave very similarly, however, an inspection of the pair distribution functions show that there is small differences between pair distributions functions for Like-Like pairs of molecules and those for Like-Unlike pairs (see for instance Fig. 13). The comparison between pair distribution functions computed for λ = 1.0 and 0.75 shows that this parameter have an important influence on the local structure (cf. Fig. 13). In particular, when λ = 1.0, interactions between the dipoles µ 1 of two enantiomers are favoured and so an asymmetric arrangement of the molecular axis.

Homochiral

Q αβ = 1 2N N i=1 (3x i α xi β -δ αβ ) (34)
There are three important parameters in this chiral model : λ, α and L. So far in this section in the simulations reported (and also in numerical computations done in subsection 2.2), we have set always L/σ = 1/4 and showed how energy minima and structures depend on λ and α. They will depend also on the ratio L/σ : for instance, if L/σ is small, chains with colinear aligned total dipole moment will be favoured, while, if L/σ is close to 1, structures with colinear aligned molecular axis will be favoured, this behaviour is illustrated on Fig. 14.

Perspectives

A full determination of the structures and clusters at low densities and temperatures will require a larger sampling of the phase space ; this can be achieved with longer runs and by implementing trial cluster moves [75][76][77]. Nematic and helical-cholesteric order parameter and also the bond orientational order parameters [78] will be useful to investigate the properties of structure at low density as function of the enantiomeric excess and model parameters, homochiral systems have interest by themselves. Same order parameters and simulations in the isobaric ensemble (N P T ) will serve to study the crystal phases of homochiral systems and also as functions of ee. Works and computations in these directions are already in progress.

Another chiral model may be defined from the computation of the multipole moments done in section 2, and in agreement with earlier studies on chiral discrimination due to multipole-multipole interactions [START_REF] Craig | Discriminating Interactions Between Chiral Molecules in Topic in current Chemistry[END_REF]44,45]. According to Eqs. (4,7), the three components, Q 10 , Q 21 and Q 22 are independent, thus an hard sphere model with dipole-dipole, dipole-quadrupole and quadrupole-quadrupole interactions [79,80] can be built by using these three multipole components as parameter for the model instead of (λ, α, L). Such model will be slightly different from the bi-dipolar hard sphere model (and also from the four charges model), since higher order multipole-multipole interactions will not be taken into account in this model. Such systems has not been studied in the present work. With the bidipolar model, one may easily build a model of chiral liquid crystal by changing the short ranged steric repulsion by replacing the hard sphere potential by a spherocylinder [81,82] or ellipsoid hard potential [83,84], or even by a soft Gay-Berne potential [85,86]. For these chiral liquid crystal models with two dipoles, all the results of section 2 hold and in particular : the computations of the multipole components and interaction energies between to enantiomers. The modifications of the phase diagram of these systems with the introduction of these chirodiastaltic interactions can be of particular interest for the physics of liquid crystals.

F

4 Figure 1 .

 41 Figure 1. Representation the molecule-fixed local frame for the Sinister bidipolar hard sphere enantiomer (Sα-HS). The frame is defined from vectors û, µ 1 and µ 2 via Eqs.(1-3) ; | O 1 O 2 |= 2L < σ. The thick dashed vector is the projection of µ 2 in the plane (O 1 , êz, μ1 ).

Figure 2 .

 2 Figure 2. Representation of the spherical components of the dipole and quadrupole moments as function of the chiral parameter χ. The spherical components, represented as function of χ, are multivalued because of the proprety χ(α) = χ(π -α) (cf.Eq.(2)). Enantiomers Rα-HS and Sα-HS have equal spherical components when an even symmetry occurs, that is to say for P and Q 21 . One may note that for λ = 1 one has Q 21 = 0.

Figure 3 .

 3 Figure 3. Representation of the spherical components of the dipole and quadrupole as function of α. As shown below, the knowledge of P , Q 21 and Q 22 as function of the parameters of the model determine fully all components of the spherical multipole tensor of an arbitrary order l (see below).

  given in Fig.2(b), as function of χ, and in Fig.3(b), as function of α.

r

  = σ > 2L).

Figure 4 .

 4 Figure 4. Representation of the two functions f 21 and f 22 as function of the reduced chiral parameter χ * = χ/Lµ 2 and as functions of α.

11 Figure 5 .

 115 Figure 5. Configurations of particular interest in dipolar system. (a) Colinear aligned total dipole moments, these kind of configurations between two dipolar particle is responsible for chains formation at low temperature and low density : for this configuration, one has θ = 0 and θ 2 = 0. (b) Antiparallel total dipole moment : for this configuration, one has θ = π/2 and θ 2 = π.

Figure 6 .

 6 Figure 6. Interaction energies for two enantiomers in configurations given by Fig.5 (a) at contact z = σ. The angle φ is defined by cos φ = ûa. ûb . In all computations reported in these figures, µ = 1.0, L = 0.25σ and the curves correspond to an increase in the value of α by 15 o ; the red lines correspond to α = 0 o , 45 o , 90 o and 135 o , these values are indicated at edge of curves in figures. The energy < E XY > φ is the average energy over the orientation ûb of enantiomers Y . (a) X and Y ≡ Rα and λ = 1.0 ; (b) Y ≡ Sα and λ = 1.0 ; (c) X and Y ≡ Rα and λ = 0.75 and (d) Y ≡ Sα and λ = 0.75.

Figure 7 .Figure 8 .

 78 Figure 7. Interaction energies for two enantiomers in configurations given by Fig.5 (b) at contact r = σ for λ = 1. In all computation, µ = 1.0, L = 0.25σ and notations are the same as in Fig.6. For the configuration given by Fig.5 (b), the angle φ is defined by cos φ = ûa. r and the angle φ 2 by cos φ 2 = ûa. ûb . (a-b) : φ = 0 o and (c-d) : φ = 45 o .

Figure 9 .

 9 Figure 9. Relative minimum of interaction energies between two enantiomers Rα in antiparallel total dipole configurations with φ = 0 as functions of α. Open circles correspond to λ = 0.75 and filled squares to λ = 0.5, the black solid lines are for φ 2 = 0 and the red dashed lines for φ 2 = π. For λ = 0.75, one has αm 166 o and for λ = 0.5, αm 150 o .

Figure 10 .

 10 Figure 10. Representation of g 000 (r) for several value of the pressure obtained by Monte-Carlo sampling of the isobaric (N P T ) ensemble of homochiral-Rα systems for N = 1000, µ = 0.5, α = 60 o and λ = 1.0.

Figure 11 .

 11 Figure 11. Snapshots of homochiral-Rα and racemic systems in external electric field. For all systems N = 512, µ = 1.0, ρ = 0.05, λ = 1.0 and E = Eêz (E = 10) ; the solid vertical line indicates the direction of the field. (a-b): homochiral-Rα with α = 35 o ; (c-d): racemic Rα-Sα system with α = 35 o and (e-f)homochiral-Rα with α = 150 o . Rα-HS enantiomers are represented in blue and Sα-HS enantiomers in red. In (b), (d) and (f), we represent, for the same configurations, molecular axis and segment in the direction of dipoles for each enantiomer ; dipoles of Rα-HS enantiomers are represented in blue and dipoles of Sα-HS enantiomers in red.

Figure 13 .

 13 Figure 13. Representation of projections of molecular pair distribution functions computed by Monte Carlo sampling of the (N V T ) ensemble of racemic Rα-Sα systems. In these computation, N = 1000, µ = 1.5, ρ = 0.1 and α = 30 o ; (a-b) λ = 1.0 and (c-d) λ = 0.75. The functions represented in (a) and (c) are computed for the set (O, P ) and the functions represented in (b) and (d) for the set (O 1 , μ1 ). In all figures, pair distributions functions for Like-Like pairs of molecules are represented in black and those for Like-Unlike pairs in red.

Figure 14 .

 14 Figure 14. Snapshots of chiral dipolar hard sphere systems for homochiral-Rα system with N = 512, µ = 1.5, ρ = 0.05, α = 30 o and λ = 1.0 and for two value of L. The notations are the same as those in Figs.11. (a-b): 2L = 0.05σ ; (c-d): 2L = 0.95σ. In (a) and (c) hard spheres are represented as balls and molecular axis O 1 O 2 by a black segment. No external electric field are applied to the systems.

  10 , Q 21 and Q 22 ) allows to compute all multipole spherical components of the model. The knowledge of all multipole moments is important in the definition of the model and in the

	February 18, 2009 19:25 Molecular Physics Mazars˙TMPH˙2009˙0027
	3
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1

 1

		< E XY > φ		φm
	α	λ = 1	λ = 0.75 λ = 1	λ = 0.75
	0 o	-5.794	-4.436	0 o	0 o
	15 o	-5.695	-4.362	11.2 o	10.4 o
	30 o	-5.506	-4.145	22.4 o	20.8 o
	45 o	-4.945	-3.800	33 o	31 o
	60 o	-4.345	-3.350	43 o	40.2 o
	75 o	-3.647	-2.825	51.6 o	48.4 o
	90 o	-2.897	-2.263	59.2 o	55.4 o
	105 o	-2.147	-1.701	65.8 o	61 o
	120 o	-1.448	-1.177	71.4 o	65.2 o
	135 o	-0.848	-0.727	76.6 o	67 o
	150 o	-0.388	-0.381	81.2 o	62.8 o
	165 o	-0.098	-0.164	85.6 o	11 o

. Average energy over the orientation ûb of enantiomers Y and location φm of the orientation of ûb , for which the minimun energy is reached, when Y ≡ Rα for configurations given by Fig.

5

(a) at contact z = σ.

Acknowledgement

The author acknowledges computation facilities provided by the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) under projects 0682104 and 0992104. In the middle of nineties Jean-Jacques was my PhD supervisor ; I have to say that I have been very lucky to benefit of his broad knowledge in the physics of liquids and in computer simulations. Therefore, I am extremely glad to dedicate this paper to Jean-Jacques Weis in this special issue of Molecular physics in his honour.

Table 3. Preliminary results for Monte-Carlo computations in the isobaric (N P T ) ensemble of homochiral-Rα systems. These computations are performed with N = 1000, µ = 0.5, α = 60 o and λ = 1.0. Notations are the same as in Table 2, p is the reduced pressure and < ρ > is the average value of the density. for the molecular axis û of enantiomers.

To study the structure of the systems we have computed the angular projections on rotational invariants of the molecular pair distribution function [59,61,73,74] ; in the present work, we have computed only g 000 (r), h 110 (r) and h 112 (r) for two sets of points and vectors. The first set is the molecule center and the orientation of In this appendix, we give general formulas for the multipole expansion of a dipole distribution µ(r) located in a finite region Ω of the space. From classical electromagnetism theory [START_REF] Jackson | Classical Electrodynamics[END_REF][START_REF] Gray | Theory of molecular fluids[END_REF] and choosing the origin inside Ω, the electrostatic potential Φ(r), due to the dipole distribution in Ω, is given by

Since we are interested in the potential at r outside the dipole distribution, we set r > r and r and r, the unitary vectors that define directions of r and r. The relation

is reported in Eq.(A1) and the multipole expansion for a dipole distribution is given by

with Q lm the m th component of the spherical multipole moment tensor of order l defined by

Eqs. (A3-A4) are the multipole expansion for a dipole distribution. In the following, we give explicitly the Q lm up to l = 3 for any dipole distribution. From Eq.(A4), and since any dipole distribution can also be represented as a real charge density, it is easy to show that

A.1. The total charge in Ω : Q 00 .

From Eq.(A4), we have

Of course, this result is obvious for any dipole distribution.

A.2. The dipole moment tensor : Q 1m .

The components Q 11 and Q 10 are given by

thus, we find

Where P is the total dipole in Ω given by

This result is exactly the same as the one obtained for a charge distribution ρ(r) with the total dipole defined by

The same derivation gives for the spherical components of the quadrupole

The relations between spherical and cartesian components are given by [START_REF] Gray | Theory of molecular fluids[END_REF],

(q xx -q yy + 2iq xy )

and the traceless quadrupole moment tensor in cartesian coordinates can be written as

The spherical components of the octopole moment may also be computed from Eq.(A4) ; after some algebra, we found


John Jairo Molina
email: john.molina@etu.upmc.fr

Carlo Pierleoni

Barbara Capone

Jean-Pierre Hansen

Igor Oliveira

Oliveira

John Jairo Molina

Igor Saulo

Santos De Oliveira

Crystal Stability of Diblock Copolymer Micelles in Solution

Keywords: diblock copolymers, micelles, self-aggregation, effective potentials diblock copolymers, micelles, self-aggregation, effective potentials

published or not. The documents may come

Introduction

Diblock copolymers are known to self-assemble into ordered or disordered supramolecular aggregates or microphases, including micelles, lamellae, cylinders or bicontinuous structures, both in the melt [1,2] and in selective solvents [3,[START_REF] Hamley | Block Copolymers in Solution[END_REF]. Selfconsistent mean-field theory, which is well adapted to melts [START_REF] Matsen | [END_REF], does not apply to semi-dilute solutions of copolymers, where concentration fluctuations are dominant, as signalled by a large osmotic compressibility. Molecular simulation is well adapted to copolymer solutions, but due to the large number of copolymers required to observe micro-phase separation, simulations of fully atomistic, monomer-level models are restricted to rather short copolymer chains [6][7][8].

To overcome this limitation and model long chains, we have developed a systematic coarse-graining strategy which reduces the computational effort by orders of magnitude [9][10][11]. The basic idea generalises earlier work on homopolymer solutions, where each polymer is represented by a single blob, and interactions between polymer coils reduce to an effective pair potential between the centres of mass (CM) of the coils [12][13][14].

The effective pair potential v(r) is calculated from a fully microscopic, monomerlevel simulation of two isolated polymer chains (corresponding to the infinite dilution limit); v(r) is easily extracted as the potential of mean force associated with the measured CM-CM pair distribution function g(r), according to:

v(r) = -k B T ln g(r) (1)
where k B is the Boltzmann's constant and T the temperature. Even for polymers with excluded volume correlations between monomers, v(r) turns out to be soft and finite at full overlap (i.e. for r = 0) [13][14][15]. In fact v(r) turns out to be accurately represented by a gaussian function of amplitude v(r = 0) ≃ 2k B T and width practically equal to the polymer radius of gyration R g . At higher polymer concentration v(r) is extracted from g(r) by the hypernetted chain (HNC) integral equation [START_REF] Hansen | Theory of Simple Liquids[END_REF] and turns out to depend weakly on concentration [14]. This concentration dependence of the effective interaction may be eliminated by switching to a multi-blob representation of long polymer chains [START_REF] Pierleoni | [END_REF].

Returning to copolymers, the minimal coarse-grained model of an AB diblock copolymer is a "soft dumbbell" representation, where each of the two blocks is represented by a single effective blob, and the two blobs are tethered by an entropic spring [9][10][11]. In references [10] and [11], this representation was applied to a simple athermal copolymer model, mimicking AB copolymers in a selective solvent, with ideal (random walk) A-blocks (corresponding roughly to θ-solvent conditions) and self-avoiding walk (SAW) B-blocks (corresponding to good solvent conditions); A and B blocks are taken to be mutually avoiding, and the corresponding copolymer solution is referred to as the ISS model. The model is characterised by three effective intermolecular blob-blob pair potentials v AA (r), v AB (r), and v BB (r), and by an intramolecular tethering potential φ AB (r), which are determined from the low density intermolecular pair distribution functions g αβ (r) and the intramolecular distribution function s AB (r) by an exact inversion procedure [9]. The resulting soft-dumbbell representation was explored by extensive Monte Carlo (MC) simulations in refs. [10] and [11] for several values of the size ratio:

f = L A L A + L B (2)
where L A and L B are the lengths of the A and B strands. It was shown that at an f -dependent critical micellar concentration (cmc), the copolymers self-assembled into polydisperse spherical micelles with A-cores and B-coronae, and typical aggregation numbers n ≃ 100. The micelles first arranged into a disordered, fluid-like structure, but upon further increasing the copolymer concentration, the micelles are found to undergo a disorder-order transition to a crystalline structure. The symmetry of the crystal which is observed is strongly dependent on the number of self-assembled micelles which form in the simulation box and the compatibility of the structure with the imposed periodic boundary conditions [10,11]. The crystals which form may hence be merely metastable, and it would clearly be of interest to determine the most stable (lowest free energy) structure at a given micelle density as determined by the copolymer concentration and the mean aggregation number n. This can be achieved by going one step further in the coarse-graining strategy, i.e. by determining the effective pair interaction between micelle CM's, and then by carrying out MC simulations of large arrays of micelles to determine the relative stability of standard crystal structures, including FCC, BCC, and diamond structures, as well as the less familiar A15 structure [18,19], which has been shown to be the most stable micellar crystal in the melt by a geometric argument and a self-consistent mean-field calculation [19][20][21]. The same phase has been shown to be stable for star polymer solutions at very high densities, at least at zero temperature [22]. Since these arguments do not, a priori, apply to diblock copolymer solutions, the main objective of the present paper is to investigate the stability of the above crystal lattices by numerical simulation. The inversion procedure used to extract the effective micelle-micelle pair potential is detailed in Section 2. MC results for the micelle pair distribution functions in the ordered phases and for the Lindemann ratios are presented in Section 3 while free energy calculations are detailed in Section 4. The predictions for the relative stabilities of the four crystal structures and the disordered (fluid) phase are summarised in Section 5 while a discussion and conclusions are presented in Section 6.

Extracting effective pair potentials from structural data

In references [10] and [11] we reported the results of extensive MC simulations of systems of 5000 "soft dumbbells" representing ISS diblock copolymers over a range of size ratios 0.2 ≤ f ≤ 0.8 and of copolymer densities ρ/ρ * , where ρ is the number of copolymers per unit volume, while ρ * = 3/4πR 3 g is the overlap density, i.e. the density at which neighbouring copolymers start to overlap; R g is the radius of gyration of an isolated diblock copolymer, and will henceforth be our unit of length. Note that for sufficiently long copolymers (L = L A + L B ≫ monomer size), their size uniquely determines R g , so that any copolymer size dependence scales with R g . For copolymer concentrations beyond the cmc, when most copolymers have aggregated into micelles, the MC simulations yield the pair distribution function g(r) of micelle CM's. This g(r) must be understood as being averaged over all micelles, irrespective of their size (i.e. their instantaneous aggregation number, which fluctuates) and over all configurations of the remaining copolymers.

According to a uniqueness theorem [23] there is a unique effective pair potential associated with a measured g(r), irrespective of the microscopic Hamiltonian governing the underlying physical system. We are hence faced with an inverse problem, namely how to extract the effective micelle-micelle pair potential from the pair distribution function provided by the MC simulations of self-assembling soft dumbbells representing the diblock copolymers. Several inversion methods, usually based on fluid integral equations, have been implemented in the past with varying degrees of success [24][25][26][27][28]. The most accurate of the inversion schemes involves an iterative predictor-corrector procedure based on MC simulations [26,27]. All inversion procedures require the input g(r) to be known over a wide range of inter-particle distances, up to values of r where g(r) has converged to its asymptotic value 1. This is unfortunately not the case for the micelle-micelle g(r) generated in the MC simulations of references [10,11] because of limitations on system size; even with 5000 copolymers, the number of micelles which form is of the order of 50, and for such a system size, g(r) still differs significantly from 1 at densities close to freezing, even at the largest accessible distance r = L/2, where L is the size of the cubic simulation cell. Hence, before attempting an inversion, the MC data for g(r) must be extended to distances r > L/2. Such an extension is easily achieved, provided the underlying pair potential v(r) is known [29], which is precisely not the case in the present situation.

We have overcome this problem by devising two independent extrapolation schemes. The first scheme is based on the observation that the static structure factor S(k), which is the Fourier transform (FT) of the pair correlation function h(r) = g(r) -1, is expected to behave smoothly at small wavenumbers k ≃ k min = 2π/L, which is the smallest wavenumber accessible in direct MC estimates of S(k) calculated via [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

S(k) = 1 N ρ k ρ -k = 1 + 1 N i =j e i k•( ri-rj) (3)
where N is the mean number of micelles in the simulation box. If the (unknown) pair potential is assumed to be of sufficiently short range, then S(k) admits a Taylor expansion in even powers of k at small k:

S(k) = S 0 + S 2 k 2 + S 4 k 4 + S 6 k 6 + S 8 k 8 + O(k 10 ) (4)
where S 0 = S(k = 0) is related to the osmotic compressibility of the micelles. If S mod (k) refers to the above Taylor expansion, and the MC data for k ≥ k are referred to as S dat (k), then the total structure factor is taken to be:

S(k) = S mod (k), k < k S dat (k), k ≥ k (5)
with k chosen to ensure continuity of S(k). The corresponding pair correlation function is then given by the FT of (5)1 , i.e.:

ρh(r) = 4π (2π) 3 e k 0 dk k sin (kr) r [S mod (k) -1] + ∞ e k dk k sin (kr) r [S dat (k) -1] (6)
If h dat (r) are the MC data for the pair correlation function, available up to r = L/2, we define the two functions:

Φ dat (r) = ρh dat (r) - 4π (2π) 3 ∞ e k dk k sin (kr) r [S dat (k) -1] (7a) Φ mod (r) = e k 0 dk k sin (kr) r [S mod (k) -1] (7b)
In view of eq. ( 6), we now require that: This is achieved by minimising the following weighted mean square deviation with respect to the expansion coefficients S 2l in eq. ( 4):

Φ dat (r) = Φ mod (r) (8
χ 2 = kmax k=kmin [S dat (k) -S mod (k)] 2 δS 2 dat (k) + rmax r=rmin [Φ dat (r) -Φ mod (r)] 2 δΦ 2 dat (r) (9)
where δS 2 dat (k) and δΦ 2 dat (r) are statistical errors of the MC data; k min = 2π/L, k max = k, r max = L/2 and r min is the distance at which the MC-generated g(r) differs significantly from zero.

An alternative extrapolation procedure is based on the assumption that the (unknown) v(r) is of finite range, i.e. v(r) = 0, r > r c . If r c is the node of the MCgenerated correlation function h dat (r) closest to L/2, then we adopt the following (approximate) closure relation:

h(r) = h dat (r) ; r < r c c(r) = 0 ; r > r c ( 10
)
where c(r) is the direct correlation function related to h(r) by the familiar Ornstein-Zernike (OZ) relation [START_REF] Hansen | Theory of Simple Liquids[END_REF] (in k-space):

h(k) = c(k) + ρ c(k) h(k) (11)
Eqs. ( 10) and ( 11) form a closed set, which can be solved by iteration to yield h(r) for r > r c and c(r) for r < r c . The two extrapolation procedures yield similar results. Yet another extrapolations scheme, based on fitting the available MC data to an exponentially damped sinusoidal function, was found to be inadequate because the simulations yield too few oscillations in g(r), due to the small system size. Once h(r) and c(r) are known for all r, the effective potential between micelles can be extracted from the HNC closure:

v(r) = k B T [h(r) -c(r) -ln (1 + h(r))] (12)
This closure is only approximate, and has proved inadequate in earlier attempts to extract pair potentials from structural data of atomic liquids, like liquid metals. However the effective interaction between micelles is expected to be much softer than atom-atom potentials and to stay finite, even at full overlap, since the pair interactions of the self-assembled soft dumbbells are gaussian-like, with a range ≃ R g and an amplitude ≃ k B T . The HNC closure has been shown to be very accurate for such soft, bounded potentials [30], and hence we believe that eq. ( 12) is sufficient to extract reliable estimates of v(r) without resorting to more advanced inversion schemes which are necessary for strongly repulsive potentials between atoms [26,27].

The accuracy of the HNC inversion procedure is illustrated in Fig. 1, where the pair distribution function g(r) determined in the dumbbell level MC simulations, which is the starting point of our inversion procedure, is compared to the g(r) calculated with the resulting effective micelle-micelle potential v(r) (shown in Fig. 3). The agreement is seen to be excellent (a similar conclusive comparison for f = 0.5 and ρ/ρ * = 5 was made in ref. [10] and for f = 0.6 and ρ/ρ * = 3.5 in ref. [11]). It is worth stressing that micelle-level MC simulations are considerably faster (by two orders of magnitude) than the effective dumbbell-level simulations.

A comparison between the v(r) extracted from MC data for g(r) via eqn. ( 12), based on the two extrapolation procedures described earlier is made in Fig. 2. The agreement is seen to be quite satisfactory, and we will henceforth only present results obtained from MC simulations using the effective potentials derived from the S(k) extrapolation scheme. As expected, the effective potential is seen to be quite soft. Nonetheless there remain large uncertainties for micelle CM-CM distances r 2R g , because the pair distribution function g(r) is practically zero when micelles overlap strongly, due to a potential barrier of several k B T . This is not a major problem for our subsequent results, because strong micelle overlap is a very rare event due to the strong (but finite) repulsion between the coronae.

It is important to stress that the extrapolation and inversion schemes of the pair structure only apply to disordered (fluid) states of the micelle solutions, since they explicitly assume continuous translational invariance, i.e. a homogeneous oneparticle density, which does not hold for ordered (crystalline) states. In the remainder of this paper we will use the micelle-micelle g(r) generated in the MC simulations of references [10] and [11] at the highest f -dependent densities where no spontaneous crystallisation of the micelles was observed, to determine the effective pair potential v(r) for three values of f (f = 0.4, 0.5, and 0.6). The corresponding effective pair potentials are plotted in Fig. 3, together with a v(r) at a lower density (for f = 0.4), to illustrate the strong density dependence of v(r). For r 3R g , where statistical uncertainties become very large, the potentials are extrapolated linearly to r = 0, yielding values at full overlap of a few tens of k B T . It is these potentials which were used to obtain the data reported in the following sections. We have checked that different extrapolations do not lead to any significant differences in the results.

Micellar crystal structures and Lindemann ratios

The effective micelle-micelle potentials determined in Section 2 have served as input in extensive MC simulations of samples of several hundred micelles confined to a periodically repeated cubic simulation cell. In order to explore the four crystal structures FCC, BCC, A15, and diamond (2 interpenetrating FCC lattices displaced along the body diagonal of the unit cell by 1/4 of its length), we have placed N effective micelles on the equilibrium lattice sites for each of these crystal structures to construct initial configurations. Thereafter new configurations were generated and accepted or rejected according to the usual Metropolis algorithm [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. For convenience we use the reduced copolymer density ρ/ρ * = 4πρR 3 g /3 as density variable. The corresponding reduced micelle density follows by dividing the reduced copolymer density by the mean aggregation number n. For the original system of soft dumbbells the initial configuration is irrelevant to study the fluid phase of self-assembled micelles. However the spontaneous crystallisation of micelles upon increasing density is difficult to trigger and in general provides imperfect structures due to finite size effects. In order to study the properties of crystalline states of different symmetry we have prepared the system in those states as follows. We have placed the system in a periodic external field, acting only on the A blob of the dumbbells, which is attractive and has the symmetry of the required crystalline structure. Starting from a random configuration of dumbbells in the presence of such a field we have observed, during a standard Monte Carlo run, a rapid formation of polydisperse micelles around the lattice sites for all structures considered except the diamond structure. In the latter case the system prefers to occupy only the sites of one of the two interpenetrating FCC lattices rather than forming smaller micelles on each lattice site. This is strong indication that the diamond structure is not a stable one for our system. Upon releasing the external field, crystalline arrays of dumbbell micelles in FCC, BCC and A15 structures remain stable if the density is large enough (see discussion below). At lower density the prepared crystals are found to melt rapidly leading the fluid equilibrium state. Fig. 4a shows the micelle-micelle pair distribution function in the effective micelle representation calculated with an initial FCC configuration (represented by the vertical segments of length equal to the coordination numbers). Clearly the initial crystal has melted and g(r) is typical of a dense fluid. Identical g(r)'s (within the very small statistical uncertainties) are obtained starting from each of the perfect crystal configurations.

Fig. 4b shows the corresponding structure factor S(k). The amplitude of the main peak is close to 2.85 which suggests, according to the Hansen-Verlet criterion [START_REF] Hansen | [END_REF], that the liquid is very close to freezing. Fig. 4b also shows S(k) obtained by a full soft dumbbell simulation. Because the latter is much more computer intensive, the statistical uncertainties are much larger, but the agreement between the structure factor, obtained within the two representations, is excellent, vindicating once more our coarse-graining strategy and the inversion procedure. The good agreement between the two representations is observed throughout the micellar liquid phase. An example of a state (f = 0.4, ρ/ρ * = 5) close to the cmc is shown in Fig. 5. Near the cmc, S(k) exhibits a central (k = 0) peak associated with the incipient clustering of the dumbbells. It is interesting to note that the coarse-grained effective micelle representation clearly preserves this feature.

Pair distribution functions obtained, starting from the four crystal structures as initial configurations, for f = 0.4 but at the higher density ρ/ρ * = 7 are shown in Fig. 6. The peak positions for the FCC, BCC, and A15 crystals are seen to coincide with the corresponding positions of neighbouring shells in a perfect lattice. This is no longer true when the starting configuration is a perfect diamond lattice. The final structure is clearly no longer crystal like, but resembles that of a disordered, glassy system. The diamond structure is thus unstable, and this is confirmed by the time evolution of the mean-square displacement of the micelles from their initial positions, to be discussed later in this Section. Note that a similar instability has been found for the full dumbbell representation as mentioned above.

Very similar behaviour is observed for f = 0.5, ρ/ρ * = 6 and f = 0.6, ρ/ρ * = 4 (not shown). While the FCC, BCC, and A15 structures appear to be stable under those conditions, the diamond lattice always melts into a disordered structure. The reason for this appears to be the low (four-fold) coordination which forces the first shell of neighbours to be at a distance from a central micelle well inside the repulsive part of the effective pair potential v(r), while for the other three crystal structures, the corresponding distances fall into the attractive part of v(r).

The previous observations are confirmed by inspection of the static structure factor S(k), shown in Fig. 7 for f = 0.5, ρ/ρ * = 6 for both the soft dumbbell and the effective micelle representations. The agreement between both sets of data is seen to be reasonable, keeping in mind that the effective micelle representation ignores the polydispersity of micelles observed in the full dumbbell representation.

The S(k) of the FCC, BCC, and A15 crystals are dominated by Bragg peaks at their expected locations, while the structure factor resulting from an initial diamond structure no longer exhibits Bragg peaks, but resembles that of a dense disordered fluid or glass. Similar behaviour is observed for the other two cases studied, f = 0.4, ρ/ρ * = 7 and f = 0.6, ρ/ρ * = 4. Within harmonic theory, the amplitude of the Bragg peaks is governed by the Debye-Waller factor proportional to the mean-square displacement of particle positions r i from the equilibrium lattice sites r i,0 :

W = k 2 6 γ 2 = k 2 6 1 N N i=1 ( r i -r i,0 ) 2 (13)
namely [START_REF] Ashcroft | Solid State Physics[END_REF]:

S(k) = 1 + (N -1) e -2W (k) (14)
If first-order anharmonic corrections are included, the amplitude of the Bragg peaks is determined by [START_REF] Draeger | [END_REF]:

S(k) = 1 + (N -1) exp    -k 2 γ 2 3 + κk 4 12
γ 2 3 2    ( 15
)
where the non-gaussian parameter κ is:

κ = 3 γ 4 γ 2 2 -1 (16)
By fitting the MC S(k) data to the functional forms ( 14) or (15) we have extracted estimates of the mean-square displacements γ 2 , which can be compared to the values measured directly in the simulations. The latter estimates, plotted against the number of MC configurations, are shown in Fig. 8, in the form of the dimensionless Lindemann ratio,

L = γ 2 d ( 17
)
where d is the nearest neighbour distance in the crystal lattice at a given micelle density. In all three cases shown, L is seen to stabilise rapidly at a constant value for the FCC and A15 crystals, while it grows continuously with MC "time" when the initial configuration is the diamond structure, confirming the fact that the latter is unstable i.e. melts into a disordered structure. The BCC structure is somewhat intermediate; for f = 0.5 and ρ/ρ * = 6, the corresponding L stabilises rapidly, but at a value nearly 20% higher than for FCC and A15. In the case of f = 0.4, ρ/ρ * = 7, L BCC is more than twice as large as the corresponding FCC and A15 values, and may eventually diverge, suggesting that the BCC phase may only be metastable under those conditions. Finally, for f = 0.6 and ρ/ρ * = 4, L BCC first appears to stabilise at a value about 20% above those for FCC and A15, but when the MC run is extended to longer "times", L BCC is eventually seen to drift continuously to larger values, suggesting that the BCC phase is unstable under those conditions. To summarise the above results, we have found that in the three cases which we have examined in detail, the diamond structure is always highly unstable while the FCC and A15 phases are always stable. The BCC phase is stable for f = 0.5, ρ/ρ * = 6, metastable (and possibly unstable) for f = 0.4, ρ/ρ * = 7, and unstable for f = 0.6, ρ/ρ * = 4. According to the Lindemann melting criterion, a crystal will melt when L ≃ 0.15 (this value is not universal, but depends somewhat on the interatomic forces) [START_REF] Hansen | Theory of Simple Liquids[END_REF]. We may conclude accordingly that the states f = 0.4, ρ/ρ * = 7 and f = 0.5, ρ/ρ * = 6 are well inside the stable FCC and A15 crystal phases, or alternatively that these crystals will melt at somewhat lower micelle densities, while the state f = 0.6, ρ/ρ * = 4 is close to melting of the FCC and A15 crystals.

In Table 1 we compare the Lindemann ratios L obtained from direct MC estimates to those deduced from an analysis of the structure factors based on eqs. ( 14) and (15). Good agreement between the two estimates is only achieved when the anharmonic correction (eq. ( 15)) is taken into account. A reasonable agreement is also observed between the results of the effective micelle and the full dumbbell representations for all structures, confirming the validity of the coarse graining strategy in the symmetry broken phase.

Free Energy Calculations

In order to assess the relative stabilities of the FCC, A15, and BCC crystal phases, and of the disordered (liquid) micellar phase, we have computed the Helmholtz free energies of all competing phases for the three (f, ρ/ρ * ) states considered in the previous Section. Here we limit our study to the effective micelle representation. A free energy calculation for the full dumbbell representation, using a novel method for cluster crystals [35], would be much more demanding and is left for future work.

Since the free energy F at a given state point is not expressible as a Boltzmann weighted statistical average, it cannot be estimated in a single MC simulation. The most direct and robust method is to evaluate F by thermodynamic or Kirkwood coupling constant integration along a reversible path, starting from a reference state point for which the free energy is known [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Hansen | [END_REF]. For the free energy of the disordered (fluid) phase, it is natural to take a low density fluid as a reference, for which the free energy can be estimated from the virial expansion [START_REF] Hansen | Theory of Simple Liquids[END_REF]. Thermodynamic integration then links the free energies of the low (ρ II ) and high density (ρ I ) states via:

βF (ρ I , T ) N = βF (ρ II , T ) N + ρI ρII dρ ′ βP (ρ ′ , T ) ρ ′ 2 (18)
where P denotes the pressure, which is easily estimated by MC for each intermediate density ρ ′ using the virial theorem [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. However the effective potentials determined in Section 2 exhibit an attractive well of depth ǫ, such that ǫ/k B T 1 (cf. Fig. 3); it is hence to be expected that the thermodynamic path in eq. ( 18) will cross a vapour-liquid coexistence line, and will hence no longer be reversible. One way to overcome this difficulty is to replace the thermodynamic path (18) along a single isotherm, by two successive paths, one along a higher (supercritical) temperature T II to go from the low density state (ρ II , T II ) to the high density state (ρ I , T II ), and a second path along the isochore ρ = ρ I , to link the high temperature state (ρ I , T II ) to the final temperature state (ρ I , T I ). Along the latter path the thermodynamic variable to be integrated is the internal energy U = ∂ (βF ) /∂β, so that the final expression linking the free energies of the reference and actual states is:

β I F (ρ I , T I ) N = β II F (ρ II , T II ) N + ρI ρII dρ ′ P (ρ ′ , T II ) ρ ′ 2 + βI βII dβ ′ U (ρ I , β ′ ) N (19)
The free energies of the initial low density, high temperature state is easily evaluated using the second virial coefficient (B 2 ) approximation for the excess free energy:

β II F ex (ρ II , T II ) N = B 2 (T II )ρ II ( 20
)
where B 2 is calculated using the relevant effective pair potential,

β II v(r) = (β II /β I ) β I v(r), with β I v(r
) the dimensionless effective pair potential appropriate for the final state of interest. An alternative procedure, which short-cuts the first part of the reversible path (from (ρ II , T II ) to (ρ I , T II )) is to use the Widom particle insertion method, which yields the excess part of the chemical potential β II µ ex (ρ I , T II ), in a single MC run [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. The method works well at high temperatures, because of the weakness of the reduced effective potential β II v(r); it would be inadequate at the temperature of interest, when β I >> β II . The excess Helmholtz free energy β II F ex (ρ I , β II )/N follows then directly from the thermodynamic relation:

βF ex N = βµ ex - βP ex ρ (21)
To evaluate the free energies of the three stable or metastable crystalline phases (FCC, A15 and BCC), we have used the Frenkel-Ladd method [36], which takes a harmonic Einstein crystal as reference system. The total interaction energy of the Einstein crystal, made up of N independent harmonic oscillators is of the form:

U E ( r N ) = E M + U H ( r N ) = U ( r N 0 ) + α 2 N i=1 ( r i -r i,0 ) 2 (22)
where the r i,0 are the equilibrium lattice positions of the particles, E M = U ( r N 0 ) is the Madelung energy (the total energy of a perfect crystal), while U H ( r N ) is the harmonic energy, with α an appropriate spring constant. If,

U ( r N ) = i<j v (| r i -r j |) ( 23
)
denotes the total interaction energy of N micelles interacting by the effective pair potential v(r), one introduces a continuous set of systems intermediate between the Einstein crystal and the system of interest, with total potential energy:

U λ ( r N ) = U E ( r N ) + λ U ( r N ) -U E ( r N ) ( 24
)
where, by varying the dimensionless coupling constant λ from λ = 0 to λ = 1 we continuously switch from the reference system to the full system. The free energy of the latter is then given by the following Kirkwood λ-integration [START_REF] Hansen | Theory of Simple Liquids[END_REF]36]:

F
F = F E + 1 0 dλ U ( r N ) -E M -U H ( r N ) λ = F E + ∆F (25)
where the ensemble average < • • • > λ is weighted by the Boltzmann factor involving the interaction energy U λ ( r N ), defined in eq. ( 24). This average is computed by MC simulations for a discrete number of values of λ (typically 4 -8). The free energy of the Einstein crystal is known analytically and the integration over λ is calculated numerically, using gaussian quadrature.

To improve the numerical accuracy it is convenient to implement the simulations for a fixed CM of the total system of N micelles, thus preventing the crystal as a whole from drifting when λ → 1, and the harmonic restoring forces deriving from U H ( r N ) no longer operate [36]. Details of the method, together with the expression for the corresponding constrained Einstein crystal, can be found in reference [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. The excess free energy per micelle is given by

βF ex N = β∆F CM N + βE M N -1 - 1 N ln ρ 2π βα 3/2 - 2 N ln N - 1 N ln √ 2πn W S +1
(26) where ∆F CM is the ensemble average defined in ( 25) evaluated from MC simulations for a fixed CM and n W S is the number of micelles per Wigner-Seitz cell 1 .

The numerical results for the free energies based on the methods described in this Section are presented next.

Competing micellar phases

We now turn to numerical applications of the free energy expressions of the previous Section, in order to determine the relative stabilities of the disordered and ordered (FCC, A15, BCC) phases. As shown in Section 3, the diamond structure is unstable against melting, at least for the three states (f ; ρ/ρ * ) examined in the present work.

Consider first the free energy of the disordered (fluid) phase, which is calculated via eq. ( 19). The reference state is chosen to be at inverse temperature β II = β I /15 and density ρ II = ρ I /100, and only excess (non-ideal) contributions are considered; the ideal contribution to the free energy is trivially known. Detailed results for the various contributions to βF ex /N (eq. ( 19)) are given below for the state (0.5; 6). Using 33 and 14 equally spaced intervals for the integration over density and temperature, respectively, we found the following estimates for the dimensionless integrals: The excess free energy of the reference state is calculated from eq. ( 20); a simple quadrature yields the second virial coefficient for the state under consideration, with the result B 2 = -5.2695R 3 g , so that:

I ρ = ρI ρII dρ ′ β II P ex (ρ ′ , T II ) ρ ′ 2 = 0.1142(1) (27)
I β = βI βII dβ ′ U ex (ρ I , β ′ ) N = -12.957(1) (28)
β II F ex (ρ II , T II ) N = -0.0007 (29)
gathering results we find the total excess free energy of the liquid phase:

βF ex liq N = -12.843(1) (30)
As explained in Section 4, the free energy of the intermediate state (ρ I , β II ) may also be estimated by the Widom insertion method (eq. ( 21)), which yields:

β II F ex (ρ I , β II ) N = -0.123(2) (31)
and should be compared to the value -0.113(0) obtained by thermodynamic integration (I ρ + B 2 ρ I ). The agreement is seen to be satisfactory, and the small discrepancy on this small contribution leads to an error of less than 0.1% on the total excess free energy.

We have carried out similar calculations for the other (f, ρ/ρ * ) states, and the results will be summarised in Table 3.

We next turn to the free energies of three crystalline phases, which were calculated from eq. ( 26). Numerical values of various terms are provided in Table 2 for the state (f = 0.5; ρ/ρ * = 6). The total excess free energies of the disordered and ordered phases for all three (f, ρ/ρ * ) state points are summarised in Table 3. Note that the ideal contribution to the free energies, is not included, because βF id /N is identical for all phases at a given density. For all three states, the free energies of the FCC and A15 phases are significantly lower than those of the liquid and of the BCC phase; the latter melts for f = 0.6, ρ/ρ * = 4, and is nearly unstable for f = 0.4, ρ/ρ * = 7. In the latter case, the free energy of the A15 phase lies well below that of the FCC phase, and A15 is hence the stable phase. Note that the A15 phase has the lowest Madelung energy, but the rather large difference in Madelung energies of the FCC and A15 phases is partly cancelled by the larger entropy (the sum of the negative excess term listed in Table 3 and the positive ideal term) of the FCC phase. In the case f = 0.5, ρ/ρ * = 6, the free energy differences are smaller , but A15 still appears to be the most stable phase. In the case f = 0.6, ρ/ρ * = 4, the free energies of the FCC and A15 phases are practically identical and only slightly lower than the corresponding free energy of the liquid phase. Note that in this case the BCC crystal eventually melts (cf. Fig. 8), while the FCC and A15 phases are expected to be close to melting, in view of the value L ≈ 0.14 of the Lindemann ratio .

Discussion and Conclusion

We have investigated the stability of ordered phases of micelles resulting from the self-assembly of a simple model of diblock copolymers in a selective solvent. Our analysis is based on a coarse-grained representation of the micelles, which are treated as monodisperse soft, spherical particles interacting via central forces which derive from an effective, state-dependent potential v(r). The latter is deter- mined by an inversion of the pair distribution function of micelles, averaged over copolymer configurations. In order to cope with the relatively small numbers of micelles (typically 50) generated in copolymer-level simulations, we have developed a novel extrapolation/inversion scheme to extract reliable effective potentials from structural MC data in the disordered (fluid) phase. The v(r) turn out to depend sensitively on the copolymer size ratio f and on the overall copolymer density ρ/ρ * . This is not surprising, since the mean aggregation number n of the micelles varies significantly upon varying these two basic parameters of the athermal ISS model [11]. The soft effective repulsion between overlapping micelles may be rationalised in terms of the finite repulsive interaction between A and B blobs of the self-assembling copolymers. The physical origin of the attraction between micelles separated by more than their diameter is less obvious, and may be due to a combination of depletion of non-aggregated copolymers, and of exchange of copolymers between neighbouring micelles which leads to fluctuations in their instantaneous aggregation number.

The analysis of the micelle-micelle pair distribution functions, structure factors, mean square displacements, and free energies generated by extensive MC simulations based on the effective pair potentials v(r) for three values of the copolymer size ratio f lead to the following conclusions:

(1) For f = 0.4 (micelle core smaller than corona), the micelles remain in a disordered fluid state for ρ/ρ * = 6, but crystalline for ρ/ρ * = 7. While the diamond structure is unstable (melts), the FCC and A15 structures are stable, while the BCC phase is highly metastable ("overheated"). The A15 phase has a free energy significantly lower than the other phases, and is hence the most stable ordered phase.

(2) For f = 0.5, ρ/ρ * = 6, FCC, A15, and BCC phases are all stable relative to the disordered (liquid) phase. The excess free energies of the three ordered phases are close (within 1%), but the A15 structure has again the lowest free energy, and is hence the stable crystal phase.

(3) For f = 0.6 (micelle core larger than micelle corona), ρ/ρ * = 4, the BCC phase slowly melts, while the FCC and A15 phases are stable and have nearly identical free energies.

It must be underlined that for all three values of f , the A15 structure has the lowest Madelung energy, but this advantage is partially cancelled by the larger entropy of the FCC and BCC phases. The relative stability of the FCC and A15 phases varies with f , smaller size ratios (f < 0.6) apparently favouring the A15 phase. It would be worthwhile to check this trend in the future, by considering the case f = 0.8, for which micellar structure data are available [11]. The larger mean aggregation number (n > 100) observed for that case points to a "harder" micelle-micelle repulsion, and hence could favour FCC over A15.

The stability of the A15 phase (at least up to f ≃ 0.6) seems to confirm predictions of Kamien and co-workers for micellization in the solvent-free copolymer melt [20,21]. The somewhat unusual A15 structure (which is not a Bravais lattice) may hence be a common phase in self-assembling soft matter systems, although free energy differences with respect to other crystal structures could be very small, as suggested by the present calculation.

Returning finally to the copolymer level, the only limitation of the present coarsegraining approach (which replaces block copolymers by "soft dumbbells") is that the effective interactions between A and B blobs of an AB copolymer were initially determined in the dilute copolymer solution limit [10,11], while clearly local copolymer concentrations within micelles can be quite substantial. This limitation can be overcome by switching to a multi-blob representation of each block, similar to a recent proposal for homopolymers in concentrated solutions [START_REF] Pierleoni | [END_REF]. Work along these lines is in progress. Recently a novel method has been devised to compute the free energies of cluster crystals of soft core particles [35]. This method could allow the study of the relative stability of crystalline structures of micelles at the full dumbbell representation level and therefore more firmly establish the validity of the coarse grained effective micelle representation. a The Wigner-Seitz cell of the A15 lattice is a fourteen-faced polyhedron with 8 sites per cell [19]. 4) and ( 5)) while the points represent the potential obtained by assuming c(r) = 0 for r > rc (eq. ( 10)).

F

Figure 1 . 3 2

 13 Figure 1. Micelle-micelle pair distribution function g(r) calculated from dumbbell-level and micelle level MC simulations, for f = 0.4 and ρ/ρ * = 6. The former data have large statistical uncertainties for the range of micelle-micelle distances L 2 < r < L √ 3 2 , where the upper limit is the largest distance accessible in a cubic simulation box.

Figure 2 .

 2 Figure2. Effective micelle-micelle pair potentials v(r) extracted from data in ref.[11] using the HNC closure with the two extrapolation schemes detailed in Section 2, for f = 0.4 and ρ/ρ * = 5. The solid curve gives the (linearly extrapolated) potential obtained by fitting the structure factor S(k) at small values of k (eqs. (4) and (5)) while the points represent the potential obtained by assuming c(r) = 0 for r > rc (eq. (10)).

Figure 4 .Figure 5 .

 45 Figure 4. (a) MC-generated g(r) for f = 0.4, ρ/ρ * = 6, starting from a perfect FCC configuration of the micelles. The vertical segments indicate the positions of the first shells of nearest neighbours in the initial FCC lattice, and their height indicates the corresponding coordination number (right hand scale). (b) Corresponding structure factors S(k) (plusses). The structure factor calculated directly from dumbbell level MC simulations is shown as crosses.

Figure 6 .

 6 Figure 6. g(r) at f = 0.4, ρ/ρ * = 7, starting from perfect FCC (top left), BCC (bottom left), A15 (top right), and diamond (bottom right) lattice configurations. Vertical segments show positions of neighbouring shells in the corresponding initial perfect lattice configuration. Their heigh indicates the coordination number (right hand scale).

Figure 7 . 14 PageFigure 8 .

 7148 Figure 7. [S(k)/N ] vs k at f = 0.5, ρ/ρ * = 6, for FCC, A15, and BCC lattices in both the soft dumbbell (diamonds) and effective micelle (plusses) representations. The Bragg peaks are the large amplitude points decreasing with k according to eq. 14

Table 1 .

 1 Lindemann ratios estimated by direct MC calculations of ˙γ2 ¸(L), and by fitting eqs. 14 (L2) or 16 (L4) to the structure factor data.

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

a This system is highly metastable.

b This system is unstable.

Table 2 .

 2 Excess free energies for FCC, A15 and BCC structures; f = 0.5, ρ/ρ * = 6.

	Lattice	N	n W S	βα	βE M /N	No. Points	β∆F CM /N	βF ex /N
	BCC	432	1	10.28	-20.2833
						4	0.854	-13.419(1)
						6	0.855	-13.418(1)
						8	0.852	-13.421(1)
	FCC	500	1	13.70	-20.7511
						4	0.798	-13.508(1)
						6	0.798	-13.508(1)
						8	0.798	-13.508(1)
	A15	512	8 a	18.97	-21.0378
						4	0.555	-13.554(1)
						6	0.554	-13.555(1)
						8	0.5546	-13.5546(5)

Page 18 of 22 URL

	February 5, 2009 14:39 Molecular Physics crystal-paper-all February 5, 2009 14:39 Molecular Physics crystal-paper-all
					Molecular Physics
	18		J.J. Molina et al. Crystal Stability of Diblock Copolymer Micelles		19
		10 15				f=0.4 , ρ/ρ*=6 f=0.4 , ρ/ρ*=7 f=0.5 , ρ/ρ*=6 f=0.6 , ρ/ρ*=4
	β v(r)	5
		0
		-5	F o r 1 2 3 F o r	4	5	6	7	8	9
	P e Figure 3. Effective micelle-micelle pair potentials extracted from the HNC closure (12) and the S(k) r P extrapolation scheme (eqs. (4) and (5)) for 4 couples (f, ρ/ρ * ). e
				e r e r
					R R
					e e
					v i e v i e
							w w
							O n l O n l
									y y

: http://mc.manuscriptcentral.com/tandf/tmph

in order to perform the Fourier integrals in eq. (6), S dat (k) was determined on a uniform grid by a cubic spline interpolation of the coarser grid of MC data.

We take n W S = 1, which means that eq. (26) gives an upper bound to F ex for the A15 lattice (since it is not a Bravais lattice); this correction is negligible.

Acknowledgments

JJM would like to thank the EACEA for support provided by an Erasmus Mundus Scholarship during his Master project, BC is supported by an EU Marie Curie Fellowship, and CP thanks Schlumberger Cambridge Research and the Royal Society of London for support during visits to Cambridge.

Table 3. Per particle values for the excess free energy F ex /N (obtained from (26) using a 4 point Gaussian quadrature), internal energy U ex /N , entropy S ex /N = β (U ex -F ex ) /N , and Madelung energy βEM /N for the three ordered structures (FCC, BCC, and A15) together with the excess free energy of the disordered (fluid) state.

State

Structure