Polarity Formation in Molecular Crystals as a Symmetry Breaking Effect

The transition of molecular crystals into a polar state is modeled by a one-dimensional Ising Hamiltonian in multipole expansion and a suitable order parameter. Two symmetry breakings are necessary for the transition: the translational and the spin flip invariance—the former being broken by geometric constraints, the latter by the interaction of the first non-zero multipole with the next order multipole. Two different behaviors of the thermal average of the order parameter as a function of position are found. The free energy per lattice site converges to a finite value in the thermodynamic limit showing the consistency of the model in a macroscopic representation.


Introduction
Polar molecules located at lattice sites as found in molecular crystals represent a system to discuss the effects of symmetry breaking in the statistical calculation of a property such as electrical polarity.Given the point symmetry of the lattice and its finite size, the system will thermalize into a bi-polar state featuring zero overall polarity [1].
Starting from our previous work [2], we consider a generalized Ising model whose Hamiltonian is given as an infinite series of multipoles.The well known spin flip symmetry of the Ising Hamiltonian is found to be broken by the interaction of the first non-zero multipole with the next order multipole, when such an interaction is present.This is analogous to the effect of an external magnetic field acting on a paramagnetic system where the symmetry breaking in the Hamiltonian generates an induced magnetization in the material.As a result of the coupling of the symmetry breaking in the Hamiltonian with the translational symmetry breaking due to the free boundary conditions, a spatial ordering in the system is achieved.

Ising Model for Bi-Polar Transition
The simplest model one can imagine to reproduce bi-polar state formation in media is an Ising-like one.It consists in a three-dimensional cubic lattice Γ and the binary variables S i = ±1 defined on each of its site i.Only interaction between neighboring variables are included.The energy of a given configuration {S i } is given by the Ising Hamiltonian where the symbol i, j indicates that the sum is between the next neighboring S i only.Since the Hamiltonian (1) is invariant under the transformation S i → −S i , S j → −S j (global spin flip symmetry) the thermal average S i = 0, unless this symmetry is broken.Thus, Equation (1), in this simple form, cannot reproduce a net average polarity.In order to analyze the issue in more detail, consider two non-overlapping charge distributions centered around R A and R B , respectively.The electrostatic energy in multipole expansion can be written as: where R AB = R B − R A , k C is the Coulomb constant whose value depends on the units, and where for each l A , l B ∈ (0, ∞) the (2l A + 1) × (2l B + 1) matrix B l A ,l B is defined as where I m l (r) are the irregular solid harmonics, are the Clebsch-Gordan coefficients.and are the multipole moments of distribution A (resp.B) given by where r i is the position of the charge q i of the distribution considered, and R m l (r i ) are the regular solid harmonics.Since the Q m l have parity (−1) l ∀m, it turns out that each term in the sum of Equation ( 2) has a well definite parity.The interactions of multipoles such that l A + l B is odd are, therefore, the origin of the symmetry breaking in Equation ( 1).This condition is, however, not sufficient for the onset of a bi-polar state because, as we will show in the following, the translational symmetry of the system should also be broken in one direction at least.This is not surprising as in the thermodynamic limit we need to have a "free" (of interaction from the outside) surface where the bi-polar state shows up.
According to the above discussion, if the system has a net charge, the lowest order terms of the series Equation (2) with broken (parity) symmetry are those with l A = 0, l B = 1 and l A = 1, l B = 0, which correspond to charge-dipole interactions.If the system is electrically neutral but polar, the symmetry is broken by the dipole-quadrupole interactions, at lowest order.
Let us consider a one-dimensional array of N identical molecules located at the sites i of the lattice and assume that each molecule i can be in one of the two possible states S i differing by an overall rotation of 180 • around the center of mass, i.e., an inversion of the dipolar direction.The multipole moments will, in general, be dependent on the state of the molecule according to some relation of the kind: where the q m l are the multipole moments of molecule i in a reference state arbitrarily chosen for all molecules, S l i is the usual effective "spin" variable associated to molecule i to the power of l.Due to the parity symmetry, only l-odded multipoles are affected by molecule flip.
The Ising Hamiltonian (1) can be re-written as where A l i ,l j are given by It is straightforward to see that the transformation {l i , m i } ←→ {l j , m j } gives A l i ,l j = (−1) l i +l j A l j ,l i , and, therefore, the matrix of interactions has parity l i + l j .The diagonal terms A l i ,l i can be expressed in the more compact form where, for each l ∈ (0, ∞), By using the symmetry properties of the Clebsch-Gordan coefficients, it can easily be seen that B l is symmetric ∀l ∈ (0, ∞).For one-dimensional lattices, the matrices B l are of a particular simple form.Let us assume that the lattice extends over the z-axis, from the symmetry of the spherical harmonics and of the Clebsch-Gordan coefficients, we then get which, for each l, is an anti-diagonal matrix.
To highlight the properties of the Hamiltonian let us separate the sums in Equation ( 6) over odd and even l i and l j where the superscript g (resp.u) indicates that the corresponding sum runs over the even (resp.odd) indices.Using the definition of S i , the previous Equation simplifies to give and after swapping the index in the third sum and using the rule of transformation for A l i ,l j , we finally get In respect to the S i variables, the Hamiltonian is therefore decomposed into the sum of three terms: one is a constant, one anti-symmetric and one symmetric where Note that H g is the usual Ising Hamiltonian with coupling constants given by and, therefore, it is the same Hamiltonian as an Ising ferromagnetic or anti-ferromagnetic depending on the sign of the J ij s.

First Order Expansion
Even if, as a matter of fact, crystals are grown with neutral molecules, we consider here the first order expansion of charged molecules because, in spite of its simplicity, it contains the fundamental ingredients accounting for symmetry breaking.The more general case of second order will be treated in the next section.The first order coefficients Equation ( 7) are In Cartesian coordinates, with the z-axis assumed along the lattice, it is straightforward to recognize that A 0,0 3 are the energies of a charge-charge, charge-dipole, dipole-charge and dipole-dipole Coulomb interactions, respectively.In the previous formulas, we have denoted by q tot the total charge of the molecule and with p z the z-component of the dipole moment (up to the sign) of the molecule.
Finally, the Hamiltonian (13), at the first order approximation, can be expressed as:

Second Order Expansion
Let us turn ourselves to the more frequent situation of neutral molecules with non-zero first order moment.In addition to Equation (17), we need the second order interaction coefficients given by: The second order Hamiltonian is given by Equation ( 18) for neutral molecules reduces to H (1) = A 1,1 ∑ i,j S i S j .Note that the neutrality of the molecule restores the parity invariance of the Hamiltonian at a first order approximation.By inserting H (1) in the expression for H (2) , and after minor manipulations, we finally get This expression is formally equivalent to the Hamiltonian introduced by Bebie and Hulliger [2], with the position where E AA , E DD , and E AD are the three possible longitudinal energies of two neighboring polar molecules of a chain (see [2] for details).
The partition function of the model can be expressed by the matrix transfer method [3] as where T is the transfer matrix given by the symmetrical (for global flip) part of the Hamiltonian with β = 1/k B T, where k B is the Boltzmann constant, T the absolute temperature, and Since the real matrix T is symmetric, it can be diagonalized.Equation ( 23) can, thus, be considerably simplified.If U is an unitary matrix such that where λ 1 and λ 2 are the eigenvalues of T , given by the partition function reads where the unitary matrix U is given by the partition function can be written as and using Equation ( 27) In order to point out the asymptotic behavior of the partition function for large systems, we notice that the extensive quantities in Equation (31) can be re-written as since λ 2 /λ 1 < 1.Since the last term decays exponentially with N, this approximation is expected to hold even for systems as small as the size of a crystal seed.The asymptotic partition function reads The same arguments give the asymptotic behavior of Equation ( 30) The expectation value of a spin can now be calculated as or in the transfer matrix formalism where σ z is the Pauli matrix given by Using the diagonalized expression for T , we get and, using Equation (30), we obtain Let us consider a system with A 2,1 = 0.In order to untangle the behavior of S k , we notice that its functional dependence on k is only through the ratio of the two eigenvalues of the transfer matrix and can thus be written as where Since ρ = 0 at every finite temperature, the two regimes of S k can be singled out −1 < ρ < 0, and 0 < ρ < 1 corresponding to A 1,1 > 0, and A 1,1 < 0, respectively, giving It is straightforwardly seen that the following general properties of the average order parameter hold: As a result of the broken symmetries, a bi-polar state appears with opposite states at the ends of the chain.2. S N has the sign of A 2,1 if ρ > 0, opposite if ρ < 0.

Since S +
k has opposite values at the boundary, it must be zero in one point internal to the interval (1, N).It is readily seen that if N is odd, S N/2 = 0.

The point of coordinate (k, S +
k ) with k = (N − 1)/2 is an inflection point.

S −
k is an oscillating function of k, bounded by the two sequences S + k with opposite values of A 2,1 (also seen in Monte Carlo simulations).7. The sequence {S + k } converges to 0 as k → ∞ ∀ρ : |ρ| < 1 and to 1 if ρ = 1 (convergent sequence bounded theorem [4]).The existence of these limits makes a sound point for the consistency of the model in the thermodynamic limit (cf.Section 4).
The behavior of S k is sketched in Figure 1.
Version February 19, 2016 submitted to Symmetry 13 of 13 Figure 1.General behavior of the average order parameter vs. distance on the chain.For clarity the curve corresponding to A 1,1 > 0 and A 2,1 < 0 has been omitted as it is simply the oscillating curve plotted with reversed sign.The values ρ = ±0.9have been arbitrarily chosen.General behavior of the average order parameter vs. distance on the chain.For clarity, the curve corresponding to A 1,1 > 0 and A 2,1 < 0 has been omitted as it is simply the oscillating curve plotted with reversed sign.The values ρ = ±0.9have been arbitrarily chosen.

Higher Order Expansion
The general treatment of Section 2 is valid at every order.It is easily seen that by replacing A 1,1 with A u , A 2,1 with A gu , and A 2,2 with A g all the results obtained are still valid, for instance where In the following, to keep the notation light, we keep considering the second order expansion.

Molecular Symmetry and Symmetry Breaking in the Hamiltonian
In Section 2, we have seen that, as long as the molecules on the chain have a non-zero odd-ranked multipole l together with a non-zero l + 1 multipole, a bi-polar state may appear as a result of the broken spin flip symmetry of the Hamiltonian ( 6).This occurrence is of course closely connected with the molecule symmetry.We stress that no other assumption has been made on the molecules forming the system.The bi-polar state is, in that respect, a universal feature.From a standpoint of group symmetry, we notice that molecules with a center of inversion, i.e., of the C i point group, cannot have an l-ranked quadrupole with l odd and, therefore, no symmetry breaking can occur.On the contrary, C n , C nν , and C s have a non-zero dipole and quadrupole, and, therefore, the symmetry breaking occurs and a bi-polar state may show up.
In the next subsections, we will examine three more significant point groups and calculate the expectation value of the order parameter.The method outlined can, however, be employed for all point groups using the corresponding character tables [5,6].Attention should be paid considering the appropriate orientation of the molecules in the chain to match the symmetry of the molecule.
We have already noticed that, roughly speaking, symmetry breaking cannot occur if the building blocks do not have an appropriate symmetry, such as C i .On the other side, symmetries can greatly simplify calculation and allow major generalization.
The C n with n > 2 turns out to be the simplest molecular symmetry required to break the symmetry in the Hamiltonian, in the sense that this point group ensures the presence of a dipole and a quadrupole, both of them having only one component.
Let us consider a frame of reference with the z-axis along the chain and assume also the molecule principal axis parallel to the z-axis.In this configuration, both the dipole and quadrupole have only the z-component non-zero, i.e., q 0 1 ,q 0 2 respectively.The relevant terms A h,k in Equations ( 17) and (19) for C n , n > 2 molecules are Since A 1,1 is negative, the order parameter should have the S + behavior.

Point Group C 2v
The symmetry of this point group requires Q ±1 1 = 0 and Q ±1 2 = 0, and, therefore, the matrices of interactions are given by Since A 1,1 is negative, the order parameter shows again the S + behavior.

Oscillating Behavior: The C s Point Group
The general behavior of S k is, at second order, eventually determined by the sign of A 1,1 , specifically A 1,1 > 0 corresponds to an oscillating behavior.We shall look at this possibility more in detail.Since I 0 2 (r) > 0 ∀r, the condition for oscillating behavior reads where * denotes complex conjugate and we have used the relation q −m l = (−1) l (q m l ) * .Thus, the S − behavior should be expected when This inequality is surely verified if q 0 1 = 0. Let us consider a chain of molecules with C s symmetry, oriented with the mirror plane orthogonal to the chain direction.The z-axis is assumed, as usual, parallel to the chain.In this case, q 0 1 = 0 and the oscillation behavior should be expected.

Macroscopic Systems
A d-dimensional Ising lattice model is defined on a finite subset Γ L of integers Z d , the number of spin variable configurations {S i } i∈Γ is therefore finite.In order to extrapolate results for macroscopic system, we consider the limit of Γ L with free boundary condition for L → ∞.Note that averaged quantities depend in general on the nature of the boundary conditions.We demonstrate that in the canonical ensemble, the thermodynamic limit of the free energy per site exists and is finite for the model considered.The free energy is given by F = −k B T ln Z.Using Equation (34), the free energy per site at the second order can be expressed in the limit of large systems as At every finite temperature, therefore, the limit lim exists and is finite.Analogously, the averaged energy of the system E , which is the thermodynamical internal energy U, is given by Note the explicit appearance of the extensivity in the energy due to the factor N. The energy fluctuations are given by The specific heat is given by We notice that the non-diagonals terms of the coupling parameters A i,j , i = j, which are essential for the symmetry breaking in the Hamiltonian and, ultimately, to the onset of spontaneous polarization do not enter the fundamental extensive thermodynamic quantities.
We conclude this section by showing that the free energy per site is a concave function of the absolute temperature.This is easily accomplished by observing that, since f ∈ C 2 (R), In consideration that the thermodynamic properties at equilibrium are determined by the free energy, this last statement ensures the thermodynamic stability of the model [7].

Conclusions
Starting from a generalized Ising model with the electrostatic interactions given in multipole expansion, we have shown that the onset of polarity in a molecular system can be described as the result of two symmetry breaking effects: the translational and the global spin flip invariance.The former is of geometrical origin and is reflected by the presence of free boundary conditions, the latter is explained as the effect of specific multipole interactions.The spontaneous polarization is assessed by the thermal average of the order parameter S analogous to the average magnetization of the standard Ising model.The analysis of S as a function of the length inside the system shows two different behaviors depending on the sign of the asymmetric interaction of same odd order multipoles (e.g., dipole-dipole): an S-shaped odd function characteristic of spin correlations fading out inside the system, and an oscillating function corresponding to an alternate ordering of plus/minus spins.
The thermodynamic limit is shown to exist under very general conditions, and the concavity of the free energy per lattice site demonstrates the thermodynamic stability of the model.At a second order treatment, the asymmetric interactions (e.g., the dipole-quadrupole one), which are crucial for the Hamiltonian symmetry breaking, do not enter any of the macroscopic thermodynamic quantities such as the free energy and specific heat.These latter quantities appear to be only determined by the symmetric (e.g., dipole-dipole) interactions.Asymmetric interactions seem, therefore, to only trigger the onset of polarity by a symmetry breaking process, leaving other thermodynamic properties, including the general behavior of S , unaffected.
The extension of the previous analysis to higher spatial dimensions and to off-lattice models is a priority for us, and, to this goal, works are already on the way.Along with theoretical calculations, large scale Monte Carlo and molecular dynamics simulations with more complex force fields are running and will be the object of forthcoming publications.

1 < 0 cFigure 1 .
Figure 1.General behavior of the average order parameter vs. distance on the chain.For clarity, the curve corresponding to A 1,1 > 0 and A 2,1 < 0 has been omitted as it is simply the oscillating curve plotted with reversed sign.The values ρ = ±0.9have been arbitrarily chosen.