New solid phase of dipolar systems

The systems of molecules with a permanent dipole moment have solid phases with various crystal symmetries. In particular, the solid phases of the simplest of these systems, the dipolar hard sphere model, have been extensively studied in the literature. The article presents Monte Carlo simulation results which, at low temperature, point to the stability of a polarized solid phase of dipolar hard spheres with the unusual number of eleven nearest neighbors, the so-called primitive tetragonal packing or tetragonal close packing.


Introduction
This work is dedicated to Jean-Pierre Badiali and our friendly collaboration. The solid phases of dipolar spherical molecules have been studied both by numerical simulations [1][2][3][4] and theoretical approaches [5,6]. At low temperatures, these phases are polarized and ferroelectric. They have crystal symmetries seemingly induced by the strong dipolar head-to-tail interaction which, in polarized phases at these temperatures, modifies the unit cells of the body centered cubic (bcc), face centered cubic (fcc) or hexagonal close packed (hcp) lattices which are expected to be the stable solid phases of particles with a spherical symmetry. Hence, on the basis of the density functional theory, following the density, temperature and dipolar moment values, the body centered tetragonal (bct), body centered orthorhombic (bco), face centered tetragonal or orthorhombic lattices have been estimated to be the stable solid phases of the Stockmayer system as well as the fundamental states of the dipolar soft sphere system [7]. Furthermore, the ground state of electrorheological fluids has been investigated; in [8,9], supposing that the ground state lattice is one among simple cubic, fcc, bcc, hcp or bct, it was established that the bct lattice has a lower energy.
The solid phases of the dipolar hard sphere (DHS) system have been studied in references [10][11][12] by numerical simulations. However, at the solid phase density, in Monte Carlo (MC) simulations realized in the canonical (NVT) ensemble, the constant value and shape of the volume with periodic boundary conditions and the hard core pair potential strongly preclude in the MC sampling the initial configuration, of a given crystal symmetry, to evolve towards the configurations of a lattice with a different symmetry and lower free energy. The MC simulations in the isothermal constant pressure (N pT) ensemble where the size and shape of the volume with periodic boundary conditions change, seem to make easier the transition between the initial configuration and the configurations of a thermodynamically more stable lattice. In practice, a very tight packing of dipolar hard spheres in solid polarized phases slows or thwarts such an evolution. Furthermore, a relative stability of phases estimated from their free energies is not obtained from the NVT or N pT simulations.
In order to circumvent this last shortcoming of the MC sampling in the NVT and N pT ensembles and then to remove the uncertainty on the relative thermodynamic stability of solid polarized phases, recent works [11][12][13] have estimated the free energy differences between solid phases of specified symmetry for the DHS system. These studies rely on the thermodynamic integration scheme [14] and on the Jarzynski relation [15,16]. In the two approaches, the free energy difference ∆F BA = F B − F A between two equilibrium states A and B of a system at density ρ and temperature T is computed. ∆F BA is given in the thermodynamic integration scheme by: Here, . . . λ is the canonical average of the partial derivative with respect to λ of the internal energy U λ of the system evolving through a set of equilibrium states from the state A to state B when the parameter λ varies from λ A to λ B . Such a computation of ∆F BA supposes that the transition from U λ A to U λ B can be made by using a parametrisation with a unique parameter and that the number of degrees of freedom stays identical in the two states A and B. For the computation of solid phase free energies in [11,12], the internal energy U λ A is that of a solid harmonic model while U λ B is that of a system of particles interacting by a pair potential, the crystal symmetries of the A, intermediate and B states being supposed identical. However, since in the harmonic solid model (state A), the particles do not have rotational degrees of freedom, also, in state B, the particles do not have such degrees of freedom. Such a constraint, if the state B corresponds to a system of polar particles, implies that the particle dipole moments remain fixed (cf. [12]). The Jarzynski relation writes: in which the average exp(−W/k B T) , where k B is the Boltzmann constant, is estimated from the works W effectuated along a set of paths inducing a transition between the states A and B, paths which can involve stable and metastable states. However, clearly, the absolute free energy of state B will be obtained only when that of the state A is known. In [13], ∆F BA is estimated from equation (1.2) by computing the work W needed to progressively modify the volume shape of systems with periodic boundary conditions from an initial shape compatible with the lattice symmetry of state A to a shape compatible with the lattice symmetry of state B.
In both computation methods of ∆F BA , defects can appear in the particle arrangements when λ or the volume shape vary. In particular, in the ∆F BA evaluation scheme based on the Jarzynski relation, it is needed to check if the state B has the expected symmetry. It is one of the aims of this work to discuss to what extent this identification is unambiguous. Section 2 is devoted to presenting the details of the tests which can be used to proceed to this identification. In section 3, it is shown on the basis of these tests that at high density and high dipole moment values, a DHS stable solid phase can be of a tetragonal close packing (tcp) type.

Symmetry tests
The MC simulations for estimating exp(−W/k B T) were realized for systems of parallelepipedic volume V with periodic boundary conditions and N ∼ 1000 dipolar hard spheres interacting by the potential is a hard sphere potential of diameter σ and the second term on the r.h.s. of equation (2.1) the dipolar interaction. r i and s i are, respectively, the vector position of the sphere i and a unit vector associated with the orientation of its dipole moment µ i = µ s i , µ = |µ i | and r i j = |r j − r i |. The reduced density ρ * = ρσ 3 and the reduced dipole moment µ * = µ/ k B T σ 3 characterize the states of the DHS system. Furthermore, the value of N and the periodic boundary conditions must be chosen such that they are compatible with the crystal symmetry of the initial state A and that expected for the final state B. From the pressure tensor associated with the interparticle interactions v(r i j , s i , s j ) which includes both the contributions of v hs (r i j ) and dipolar interaction, it is straightforward, as described in [17][18][19],  to compute the work needed to change the shape of a DHS system by a series of small homogeneous deformations. These deformations generate the transition between states A and B of the same V , N and T values but with different periodic boundary conditions. Obviously, there are many types of such paths allowing one to generate a transition between the volumes of the same size but with different shapes. In [13], to maintain the system in a solid phase and to minimize the formation of defects, the transition paths between the states A and B at constant T prevented the system from melting by keeping the density variation in the limit of a few percent with respect to the identical density of the initial states A and final states B.
The identification of the lattice symmetry for a configuration relies on the pair distribution function g(r), the structure factor, the positions and numbers of particles present in the neighboring shells, as well as on the types of Voronoi polyhedra associated with the nearest neighbors of a particle. Clearly, these properties are not independent, but their simultaneous use allows one to reduce the ambiguities resulting from the fluctuations in the local arrangement of the particles.
For instance, to determine the crystal symmetry of DHS configurations at ρ * = 1.06−1.26 and µ * = 1.0−3.0, a simple way is to compare their two-body distribution functions g(r) with those of solid HS systems, computed by MC simulations in the canonical ensemble for solid states of the known symmetry. Indeed, as mentioned above, at these densities, the constant volume and hard core packing effects seem a priori to preclude any symmetry change. Figure 1 shows, at these densities, such results for µ * = 0.0, i.e., a HS system, obtained after 40 · 10 6 MC trial moves in the NVT ensemble, when, in the initial configurations, the HS are located on a bcc, bct, fcc or hcp lattice. The functions g(r) of equilibrium bcc configurations at both densities clearly differ from those computed for the other crystal symmetries, while the positions of the two first peaks in g(r) at r < 1.6 for the bct, fcc and hcp configurations are almost identical, the third peaks at r = 1.8−1.9 differing only by their amplitudes. It is only for r > 1.9 that, in particular, the g(r) functions of the fcc solid are clearly distinct from those of the bct and hcp solids, whereas the bct and hcp g(r) functions stay more similar.
Additional information to interpret these g(r) data can be obtained by computing the number of neighbors n k involved in the g(r)'s k-th peak and, also, the average locations of the neighbors around one HS in these equilibrium solid configurations. n k is defined by

33601-3
where r inf k and r sup k are the positions of the g(r) local minima surrounding its k-th peak. For a configuration of N HS, the average positionr i of the i-th neighbor can be defined as: where r j is the position of the HS particle j, while r j i is that of its i-th neighbor sorted into the increasing values of |r j − r j i |, so i = 1 is the nearest neighbor, i = 2 the second nearest neighbor, . . . . These average positions indicate whether the numbers of particles in the g(r) peaks are in agreement with those expected in the neighboring shells for the symmetry of the considered solid phase. Indeed, due to the density fluctuations, the neighbor average positions in a shell present a dispersion reflecting the g(r) peak width and, in ther i graphs, the neighboring shells are indicated by a gap or a variation of the curves slope.
Results ofr i computations are presented in figures 2 and 3 for the equilibrium HS configurations at ρ * = 1.06 and 1.26 obtained from HS initial configurations located on a bcc, bct, fcc or hcp lattice. At ρ * = 1.06, figure 2 clearly shows the important dispersion of the first, second and third neighbor shells for the noncompact initial bcc and bct configurations. From ther i plot for the bcc equilibrium configurations, up to 2.5, it appears that the two first shells overlap and the 4 to 6 th shells are largely spread. Furthermore, the n k values indicate that the g(r)'s first peak involves 14 neighbors and the second peak involves 46 neighbors. Similarly, in the bct configuration, the two first shells overlap, the g(r) first peak involving 12 neighbors and the second, third and fourth peak 6, 20 and 18 neighbors, respectively. In the compact fcc and hcp configurations, the two first peaks in g(r) correspond to 12 and 6 neighbors in agreement with the values expected for a perfect crystal. This agreement holds for the third, fourth and fifth shells in the fcc configuration in spite of the shell broadening. In the hcp configuration, the third to sixth shells overlap to form the g(r) third peak. For this solid low density, figures 1 and 2 give clear evidence that in the bct configuration, the local structure of the three first neighboring shells rearranges towards that found for the hcp configuration.  At ρ * = 1.26 compared to ρ * = 1.06, the shell spreading is reduced. In the fcc g(r) plot, the peak positions and widths correspond to the positions and expected numbers of neighbors of a perfect fcc solid up to r > 3.0. For the hcp configuration, the g(r) third and fourth peaks result from the overlaps of the third and fourth shells, and of the fifth and sixth shells, respectively. Clearly, the HS local arrangement in the equilibrium bct configuration is very similar, almost identical to that obtained in the equilibrium hcp configurations. In the equilibrium bcc configurations, the g(r) first and second peaks involve together ∼ 14 neighbors, about 10−11 in the first and 3−4 in the second peak, these neighbor numbers and the position of the second peak being similar to those expected in a bct crystal.
The Voronoi tessellation of the HS considered solid configurations give little additional information on their lattice symmetry. In the fcc and hcp configurations, at both densities, the face, edge and vertice numbers of the Voronoi polyhedra significantly differ for most HS from the values characterising the polyhedra in the fcc and hcp perfect lattices. A similar result is obtained for the bct configurations. It is only for the bcc configuration that the 14 HS neighbors, present in the g(r) first peak at ρ * = 1.06 and the first and second peaks at ρ * = 1.26, delimit, for more than 50% of HS, a Voronoi polyhedra with face, edge and vertice numbers in agreement with the polyhedra of a bcc perfect lattice.
The analysis of HS solid configurations with regard to the identification of their crystal symmetry, shows that in solid noncompact configurations, the HS local arrangement, compared to that of perfect lattices, can be strongly modified by a large overlap of the neighboring shells. These modifications can induce an evolution of the local arrangement towards that of more stable configurations, as seems to happen for the bct configurations, where the local arrangement becomes of hcp type at the two considered densities. Noticeably, in spite of the nearest neighbor shell spreading, in the tessellation of the bcc configurations, most Voronoi polyhedra have 6 faces with 4 edges and 8 faces with 6 edges disposed as in the polyhedra of the bcc and bct perfect lattices.
In conclusion, the association of g(r) data with computed values of n k andr i and, possibly, the Voronoi tessellation allows one to obtain, at the equilibrium in a solid phase, a definite characterization of the HS local arrangement, but, of the crystal symmetry only for HS compact configurations.

Symmetry of DHS polarized solids
The previous analysis made on the HS solid configurations is used now to study the symmetries of the DHS configurations at ρ * = 1.26 with µ * 0, the aim being to achieve a better characterization of these symmetries than in [13]. In the latter work, as mentioned above, the symmetries were supposed to be those associated with the periodic boundary conditions of the volume V enclosing the DHS system in the final states B of the transition paths used to compute ∆F BA by equation (1.2). In these computations, N, the numbers of DHS in V and periodic boundary conditions were chosen in such a way that the B states can be bct, bco, fcc or hcp configurations without defects. By comparing the ∆F BA values obtained with respect to an identical initial A state, namely a bcc configuration at equilibrium, it was possible to estimate the relative stability of the B states of different symmetries.
Hence, it was obtained that the more stable configuration at ρ * = 1.26 at µ * = 1.0 was bcc and at µ * = 2.0 and 3.0 bco with a ratio of the volume edges L y /L x along the y and x directions ∼ 1.10. The g(r) functions of such configurations are given in figure 4. In figure 5, similarly to figures 2 and 3, there are given the numbers of neighbors present in the peaks of these functions. At µ * = 1.0, as in the HS system at ρ * = 1.26, the first shells are modified. The first peak involves 11 neighbors and the second, third and fourth peaks 3, 6 and 12 neighbors. This important modification of the local arrangement with respect to that of the bcc crystal remains, as for the HS configuration, compatible with a tessellation made of Voronoi polyhedra with 14 faces (8 with 6 edges and 6 with 4 edges) typical of the tessellation of bcc and bct crystals. These polyhedra around a specific DHS in the configuration are strongly distorted, but when the N polyhedra are averaged, the average polyhedra are close to that of a bcc perfect crystal. On the basis of this latter result, it seems still reasonable to label these unpolarized DHS configurations as bcc.
For    [20][21][22]. In this type of configuration, the nine first shells correspond to 11 The fifth and sixth shells and eighth and ninth shells, considering their respective very close positions, are expected, at finite temperature, to overlap and, hence, leading to the 14 and 10 neighbors in the g(r) sixth and eighth peaks. The g(r) computed by MC simulation for DHS system in tcp configurations and plotted in figure 5, unambiguously confirms a full agreement between the g(r) functions of a tcp solid and those of the states B with bco periodic boundary condition obtained at the end of the transition paths.
In [13], it was concluded from ∆F BA computations that these configurations with local tcp symmetry were the most stable among the equilibrium configurations with bct, bco, fcc and hcp periodic boundary conditions obtained from transition paths starting from bcc configurations. In addition, in this reference, this result was supported by the fact that in N pT simulations, at µ * = 3.0 and ρ * = 1.26, the bct initial configurations also evolved to configurations with g(r) functions similar to those of a tcp solid (cf. figure 4 of [13]).
From this analysis of the local arrangements made on the configurations estimated to be the most stable at µ * = 3.0 and ρ * = 1.26, it is possible to conclude that in this domain of density and µ * value, the DHS solid phase has a tcp symmetry.

Conclusion
In section 2, it was shown that the identification of the lattice symmetry of noncompact HS solid phases made only on the basis of the g(r) functions is not obvious. The local arrangement in such phases around a given HS, presents large overlaps of the nearest neighbor shells, so that the n k values of the g(r) peaks for r < 2.0 do not correspond to those expected in bcc or bct solid phases, but rather those of fcc  or hcp phases. However, to determine if these HS noncompact configurations eventually would evolve towards compact configurations without defects of fcc or hcp type, prohibitively long simulations seem required. A similar uncertainty is encountered to characterise the lattice symmetry of DHS bct or bco configurations, as shown above for the DHS systems at µ * = 1.0 and 2.0 and ρ * = 1.26. Thus, as already mentioned, in [13], these configurations with the lower values of ∆F BA were labelled according to the periodic boundary conditions of the states terminating the transition paths. Then, it is remarkable that, for µ * = 3.0 and ρ * = 1.26, on the basis of the excellent agreement between the g(r), n k andr i values with those of tcp solid phase, it can be concluded almost unambiguously that this phase is a stable phase of the DHS system.