Article pubs.acs.org/JPCB
Dielectric and Phase Behavior of Dipolar Spheroids Lewis E. Johnson,* Stephanie J. Benight, Robin Barnes, and Bruce H. Robinson* Department of Chemistry, University of Washington, Seattle, Washington 98195, United States S Supporting Information *
ABSTRACT: The Stockmayer fluid, composed of dipolar spheres, has a well-known isotropic−ferroelectric phase transition at high dipole densities. However, there has been little investigation of the ferroelectric transition in nearly spherical fluids at dipole densities corresponding to those found in many polar solvents and in guest−host organic electro-optic materials. In this work, we examine the transition to ordered phases of low-aspect-ratio spheroids under both unperturbed and poled conditions, characterizing both the static dielectric response and thermodynamic properties of spheroidal systems. Spontaneous ferroelectric ordering was confined to a small region of aspect ratios about unity, indicating that subtle changes in sterics can have substantial influence on the behavior of coarse-grained liquid models. Our results demonstrate the importance of molecular shape in obtaining even qualitatively correct dielectric responses and provide an explanation for the success of the Onsager model as a phenomenological representation for the dielectric behavior of polar organic liquids.
■
INTRODUCTION Coarse-grained spheroidal models of molecules approximate the steric behavior of a wide variety of molecules, including small organic liquids,1−3 organic electro-optic chromophores,4−7 and liquid crystals (LCs)8 at low computational cost. Such models, typically implemented using the Gay−Berne potential,9 variants thereof,10,11 the Berne−Pechukas potential,12 or the method of Perram and Wertheim,13 introduce a directional dependence to the well depth (εLJ) and/or van der Waals diameter (σLJ) of the ubiquitous Lennard-Jones potential, breaking its spherical symmetry. The steric anisotropy can result in dramatically different behavior in comparison with purely spherical models such as the Stockmayer fluid14 (sphere containing a point dipole at the center). For uniaxial LCs, the molecular aspect ratio (length/width) is particularly crucial. Onsager demonstrated that smectic or nematic phases could be induced in rigid spherocylinders (cylinders with hemispherical ends) by increasing their aspect ratio in his pioneering 1949 paper15 on LC phases of colloidal materials. It was later shown that LC phases could also be formed by prolate spheroids16 and that steric anisotropy was a necessary condition but not a sufficient condition for forming LC phases.17,18 Sufficiently oblate (discotic) spheroids with axial dipoles were demonstrated to form LC phases by Frenkel,19 Patey,20,21 Zannoni,22 and others.23 Oblate spheroids with transverse dipoles24 have also been demonstrated. Further work by Baus,25 Groh,3,26 and Ayton21 demonstrated that dipolar spheroids could form ferroelectric LC phases. In addition, work by Berne2,8 and co-workers explored the effects of electronic polarizability on simple models of LC systems. Hall and co-workers have also examined the phase behavior of dipolar Janus particles.27 Development of simulation methods for these types of systems is crucial for the development of © XXXX American Chemical Society
coarse-grained protocols for simulating acentric ordering in electro-optic chromophores and other highly polar molecular systems. In addition to changes in centrosymmetric (LC) and acentric (ferroelectric) ordering at equilibrium, changes in the aspect ratio of a molecule can have dramatic effects on the induced acentric ordering of dipole moments with respect to an applied electric field. Both analytic models and Monte Carlo simulations have demonstrated that dipolar spheres tend to align more strongly with an electric field than prolate spheroids.4,28 The dramatic difference in ordering behavior has had significant implications for the design of chromophores for electro-optic applications; such chromophores are typically highly prolate due to their long π-systems.29 In electro-optics, the ordering is typically quantified by the acentric order parameter ⟨cos3 θ⟩, where θ is the angle between the molecular dipole moment and the poling field used to align the chromophores. This cubic average cannot currently be measured directly in amorphous materials and is generally estimated through simulation or measured indirectly based on electro-optic response30 or second harmonic generation.31 Acentric order may also be indirectly quantified through capacitance measurements32 of the dielectric constant ε, where
(ε − 1)E = 4π P
(1)
Here, E is the macroscopic electric field applied to the material, and P is the total polarization of the material. We further apply four simplifying conditions: (1) ε is isotropic (scalar). Received: January 1, 2015 Revised: March 23, 2015
A
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B (ε − n2)(2ε + n2) =y ε(n2 + 2)
(2) The material is composed of identical molecules with permanent dipole moment μ. (3) All polarization of the system is due to reorientation of permanent dipoles as opposed to formation of induced dipoles (negligible electronic polarizability). (4) The system is oriented such that the electric field is aligned along the z-axis of the lab frame. Under these conditions, eq 1 simplifies to ε−1=
4πρN μ⟨cos θ ⟩ 4πMz = VEz Ez
where n is the refractive index of the liquid, was developed to correct for the Debye37 model’s incorrect prediction of a ferroelectric transition at y = 1. However, while the Onsager model is capable of reasonably approximating the experimental relationship41 between the dipole moment and dielectric constant in liquids, it is qualitatively incorrect in modeling the many-body behavior of the system it was derived for (dipolar spheres) compared to results obtained from numerical simulations.14,42 The Onsager model does not predict a ferroelectric phase transition under any conditions. However, unlike the Stockmayer fluid, most real polar liquids exhibit little intermolecular correlation.36 The low intermolecular correlation in real liquids results in good phenomenological agreement between the Onsager model and real liquids, despite the Onsager model’s underestimation of correlation in the Stockmayer fluid. Here, we present work that examines the ordering behavior and degree of intermolecular correlation (as quantified by deviation from the Onsager model) for a range of dipolar spheroids. While previous work has often focused on wide ranges of aspect ratios (up to several orders of magnitude)17 and/or temperatures22 and densities,19 this work will focus on the effects of dipole moment and aspect ratio at constant temperature, density, and molecular volume. We consider a small range of different aspect ratios centered about unity and determine the degree of ordering of dipolar spheroids in a fluid compared to that of the Stockmayer fluid.
(2)
in the static (DC) limit. Here, V is the volume of the system, ρN is the number density of the dipoles, and θ is the angle between any individual dipole and the applied field, such that ⟨cos θ⟩ is the cosine of the angle between the total dipole moment M and the applied field. Within the assumptions of eq 2, ⟨cos θ⟩ is simply the normalized projection of the dipole moment along the z-axis. The third-order average ⟨cos3 θ⟩ (crucial to electrooptic materials) tends to closely track ⟨cos θ⟩.5 Furthermore, if the system is composed of spherical molecules, its thermodynamic state can be uniquely represented as a function of reduced dipole moment and reduced density14 μ*2 =
μ2 σLJ3kBT
ρ* = ρN σLJ3
(3)
Here, T is the absolute temperature, and σLJ is the LennardJones diameter of the molecule. These reduced parameters can be generalized20 to nonspherical molecules by replacing the single σLJ with the molecular semiaxes (a, b, and c), such that ρspheroid * = ρNabc. Additionally, the dipole density of any polar material can be reduced to a single dimensionless parameter33
y=
■
COMPUTATIONAL METHODS
Simulations of dipolar spheroids were run at nine different aspect ratios and five different dipole moments, holding all other thermodynamic parameters and terms in the intermolecular potential constant. Maximum rotational and translational move sizes were optimized for each system in order to maximize the RMS displacement per Monte Carlo step. Methods were substantially identical to those used in our 2010 paper.1 All simulations used self-consistent reaction field (SCRF) electrostatic boundary conditions. Run lengths were extended compared to our prior work1 to improve convergence of second-order properties such as fluctuations in the dipole moment and fluctuations in the energy.43 The simulation size was set at N = 1728 particles, and simulations were initially run for a length of 120 000 cycles, with N trial moves per cycle for a total of 2.07 × 108 configurations per simulation. The final 60 000 cycles of each simulation were used for calculating properties. Oblate spheroids with y ≥ 3 exhibited slow convergence, and the simulation length for these systems was extended to 300 000 cycles, using the last 100 000 cycles to calculate properties. The run length was further extended if simulations were still not converged to an energy difference of i ⎢ ⎠ ⎩ ⎣ ⎝
calculated self-consistently based on the most recent reaction field update. The dielectric constant was calculated using eq 2 in the calculations with an external field. In calculations without an external field, it was calculated using a Kirkwood fluctuation expansion33,44
⎤ 1 1 2(εS − 1) ⎥ μi ·μj ⎥ + 3 (μi ·μj − 3(μi · riĵ )(riĵ ·μj )) − 3 ε 2 1 + rij R cut S ⎥⎦ ⎫ ⎪ 2εS + ε μi · E⎬ − 2εS + 1 ⎪ ⎭
⟨Mz⟩E0,z =
Here, εLJ is the Lennard-Jones energy (set to 217 K), σeff is the effective Lennard-Jones contact distance calculated using the method of Perram and Wertheim,13 Rcut is the cutoff distance for calculating electrostatic interactions (27.96 Å), εS is the reaction field dielectric constant, which is self-consistently updated every cycle,1 and E is a homogeneous external field applied along the z-axis of the simulation cell. Lennard-Jones interactions were truncated at a cutoff distance of 18.64 Å, corresponding to 4σ, for reasons of computational efficiency. Each of the nine spheroid geometries was run at all five dipole densities and a constant number density of ρ* = 1, intended to be representative of common organic liquids. The NVT ensemble was chosen to ensure constant dipole densities at all aspect ratios. The thermodynamic parameters μ* and ρ* were redefined to reflect the spheroidal nature of the systems as μ*2 = μ2/(σaσ2b)kBT and ρ* = ρN(σaσ2b), where σa is the length of the dipole axis of the spheroid and σb is the length of the nondipole axes. The electrostatic parameters for each y are in Table 1, and the spheroid geometries for the different aspect ratios are listed in Table 2.
T (K)
μ*2
ρ*
2 3 4 5 6
2.42 2.96 3.42 3.82 4.19
293 293 293 293 293
1.43 2.14 2.86 3.57 4.29
1.0 1.0 1.0 1.0 1.0
A (R ) =
σa (Å)
aspect ratio
V (Å3)
1:2 4:7 2:3 4:5 1:1 5:4 3:2 7:4 2:1
5.88 5.62 5.34 5.02 4.66 4.32 4.08 3.86 3.70
2.94 3.20 3.56 4.02 4.66 5.40 6.10 6.76 7.40
0.50 0.57 0.67 0.80 1.00 1.20 1.50 1.75 2.00
52.98 52.98 52.98 52.98 52.98 52.99 52.98 52.98 52.98
(9)
1 2
∫0
∞
dx (x + 1) (R2x + 1) 3/2
(10)
This approximation adds the electrostatic effects of molecular anisotropy (distortion of the dielectric cavity around the molecule) in the absence of intermolecular correlations. In addition to ordering resulting from polarization in response to an applied field, spontaneous order was calculated from simulations at zero field. As simulations at zero field do not have an external director axis that is known a priori, ordering is quantified by means of the de Gennes46 Q-tensor
Table 2. Simulation Parameters: Spheroid Geometries σb (Å)
(8)
where a is the major (dipole) semiaxis of the molecule and b is the minor semiaxis, such that the aspect ratio is R= a/b and
Q= eccentricity
(⟨M2⟩0 − ⟨M2⟩20 )
(ε − 1)(ε + A(a , b)(1 − ε)) =y 3ε
Table 1. Simulation Parameters: Dipole Strength μ (D)
3kBT
of the total dipole moment M of the system at zero applied field, assuming a linear response of M to E0 and substituting into eq 2 to obtain the dielectric response. The two types of dielectric calculations were then compared to determine whether or not the linear response approximation was valid under each condition. The Onsager model (eq 5) was used as an uncorrelated control model. Onsager’s predictions were calculated at each dipole density, neglecting electronic polarizability as per the assumptions of eq 2 such that the refractive index was n = 1. Additionally, the Onsager model can be generalized to spheroids (as initially suggested in Onsager’s derivation38) by means of a Dirichlet correction45 on the local field factor such that
(6)
y
E0, z
3⟨μi ⊗ μj ⟩ − I (11)
2
averaged over all μi ≠ μj. The dominant eigenvalue of Q is the nematic order parameter P2 with respect to the internal director axis d of the system, which is itself the dominant eigenvector of Q.34 The corresponding acentric order parameter P1 = M̂ ·d can be obtained by projection of the total dipole moment of the system on the director and is the zero-field analogue to ⟨cos θ⟩. Finally, the radial distribution function33,47,48 N
g (r ) =
N
V ∑i ∑ j > i nij(r , r + Δr ) 2πN (N − 1)r 2Δr
(12)
and first-order dipole−dipole (angular) correlation function1,49 N
In order to obtain both fluctuation-based and direct measurements of the dielectric constant, simulations were run under both zero applied field and an applied field of E = 50 V/ μm, with the local field E=
2εS + ε E 2εS + 1
gμ(r ) =
N
V ∑i = 1 ∑ j > i cos γij(r , r + Δr ) 2πN (N − 1)r 2Δr
(13)
were calculated, where nij(r,r+Δr) is the number of particles where the center−center distance rij is between r and r + Δr and cos γ = μ̂ i·μ̂ j. The former indicates the degree of positional correlation between atoms, with the function decaying to unity in a disordered structure such as an isotropic liquid or a glass
(7) C
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
Figure 1. Static dielectric constant of spheroids of different aspect ratios (geometries shown in (1a)), calculated using both dipole fluctuations (shown in (1b)) and in response to an applied external field of 50 V/μm (shown in (1c)). Due to a large range of values, the ordinate is plotted on a logarithmic scale. Calculations were run at 293 K and ρ* = 1 using SCRF boundary conditions at dipole densities in the range of y = 2−6. Predictions from the Onsager model (eq 9) are shown as solid lines. Dotted lines denote the maximum possible dielectric response (perfect order) under the applied field, and dashed lines are guides to the eye.
results, a table of dielectric constants is also available as Supporting Information. Results at the lowest dipole density (y = 2) were within the linear response regime, with strong agreement between the applied field and fluctuation results at all aspect ratios. The dielectric constant depended strongly on the aspect ratio, forming a Lorentzian-like function peaking at an aspect ratio of unity (the Stockmayer fluid). The peak value was almost four times larger than the lowest (2:1 prolate) result, consistent with literature values for similar systems.42 Spheres and nearly spherical systems were poorly described by the Onsager model (eq 9, shown as solid curves on the graphs). However, the more prolate spheroids were reasonably well described by the Onsager model, indicating a very low degree of net intermolecular correlation.
and exhibiting repeating peaks in a system with strong positional order such as a smectic LC. Equation 13 indicates the degree of angular correlation between dipoles, with negative values indicating net anticorrelation at a given distance and positive values indicating net correlation. In an isotropic material, the angular correlation will decay to zero (no correlation), but in a ferroelectric material, it will limit to a positive value.21,50
■
RESULTS AND DISCUSSION
Static Dielectric Constant. The static dielectric constant was calculated using both the applied field and fluctuation method at six dipole densities, ranging from y = 2 (lowest) to 6 (highest); results are shown in Figure 1 along with the geometries of the dipolar spheroids; due to the large range of D
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B While results for prolate spheroids at y = 3 followed similar trends to those of y = 2, the behavior of oblates differed dramatically. Mildly oblate (4:5) spheroids transitioned to a ferroelectric phase, and 2:3 spheroids showed evidence of ferroelectric domains. Furthermore, the linear response approximation in eq 8 began to break down, with substantial disagreement between applied field and fluctuation calculations for oblate (a < b) spheroids. While this could be due to slow oscillation of ⟨M2⟩ and/or difficulty in averaging out any net dipole moment,43,51,52 results were consistent among the trajectories, and calculated dielectric constants became smaller as simulation lengths increased. Prolate spheroids with an aspect ratio of 3:2 or higher were still well-approximated by the Onsager model. Behavior at y = 4 was qualitatively similar to that at y = 3 at all aspect ratios. However, the ferroelectric domains at the 2:3 aspect ratio disappeared, and the ferroelectric response at the 4:5 aspect ratio was enhanced. The prolate spheroids exhibited much lower dielectric constants than the Stockmayer fluid, and agreement between the fluctuation and applied field methods remained strong for all a > b. The range of aspect ratios for which the Onsager model was approximately valid decreased. By y = 5 and 6, the Stockmayer fluid exhibited a very high dielectric constant, consistent with our prior results.1 Trends were otherwise similar to those of y = 4. The 1:2 oblates at the highest dipole density showed a weak increase in dielectric response to an applied field, although this was not reflected in the dielectric response calculated from dipole fluctuations. Interestingly, even at the highest dipole densities, the applied field and fluctuation-based responses for prolates closely mirrored each other. The Onsager model remained a reasonable approximation for the 7:4 and 2:1 prolates. The close match between the Onsager model and Monte Carlo results for prolate spheroids over the wide range of dipole densities raises the possibility that the Onsager model’s ability to closely approximate experimental data for polar liquids is due to a fortuitous cancellation of error; nonspherical systems (especially prolate spheroids) exhibit far weaker intermolecular correlation than spheres, such that the Kirkwood g factor33,44 remains near unity for many realistic molecular geometries. However, the purely electrostatic contribution from the shape of the Onsager cavity around the molecule, which is accounted for in the extended Onsager model (eq 9), depends only weakly on the cavity shape. Two of the most polar small organic liquids (acetonitrile and nitromethane) are highly prolate, and others such as acetone and dimethylformamide are oblate but with an in-plane dipole instead of an axial one. As such, ordering of real molecules in response to a field is frustrated by intermolecular correlation1,39 due to molecular shape. The relatively small difference between results calculated with eqs 5 and 9 (see the Supporting Information) explains why the Onsager model, although constructed using spheres, is more effective at explaining dielectric trends in real liquids than as a standalone theory of simple liquids. Spontaneous Order. The order parameters ⟨P1⟩ and ⟨P2⟩ were calculated using eq 11 for the simulations at each dipole density in the absence of an applied external field. These results are shown in Figure 2. No spontaneous ordering was observed in prolate ellipsoids at any dipole density under the conditions studied. In general, prolates that form nematic phases at ρ* = 1 often have substantially larger aspect ratios than those considered here.3,19,50,53 Furthermore, no spontaneous ordering (centro-
Figure 2. Spontaneous acentric (2a) and centrosymmetric (2b) ordering as a function of aspect ratio for dipolar spheroids at six different dipole densities at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field, as obtained from the ordering matrix Q (eq 11).
symmetric or acentric) was observed for any aspect ratio at y = 2. At the other dipole densities, the more oblate systems with aspect ratios between 2:3 and 1:2 exhibited weak spontaneous ordering, with ⟨P2⟩ typically larger than ⟨P1⟩. These phases may have slight antiferroelectric character. Strong (⟨P1⟩ > 0.3) acentric ordering was observed for 4:5 oblates at y = 3, 4, and 5 and for the Stockmayer fluid at y = 5 and 6, consistent with the singularities in the dielectric constant seen at those dipole densities (see Figure 1). This strong ordering in the absence of a field is indicative of formation of a ferroelectric20,54 phase and is most pronounced for 4:5 oblates at y = 4, where ⟨P1⟩ approaches unity. For all of these ferroelectric systems, ⟨P1⟩ was larger than ⟨P2⟩. Only two conditions (4:5 oblates at y = 4 or 5) showed a degree of centrosymmetric ordering similar to the degree of acentric ordering observed. Strong centrosymmetric order in the absence of acentric order (pure nematic or smectic phases as observed in typical LCs) was not observed over the range of aspect ratios and dipole densities studied, though 2:3 oblates at y = 3 exhibited extensive centrosymmetric order with only a weak ferroelectric response. Thermodynamics. In addition to ordering behavior, the total potential energy (U) and the potential energy E
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B contribution55 to the constant-volume heat capacity (Cv) were examined for signs of phase transitions; results appear in Figure 3. Heat capacities were calculated from fluctuations in the total potential energy at constant temperature.
and is minimized for the systems with the weakest intermolecular interactions. The heat capacity varies weakly as a function of aspect ratio at all dipole densities for prolate spheroids and spheres, reaching a weak local minimum at an aspect ratio of 5:4. The trend in heat capacity for prolates and spheres is nearly independent of dipole density, and the change in heat capacity across the isotropic−ferroelectric phase transition is modest. A substantial shift in heat capacity (Cv) is also seen as the aspect ratio is decreased below unity. The heat capacity increases rapidly as a function of aspect ratio for oblates with y > 2. The highest heat capacity is observed for 1:2 oblates at the highest dipole densities. While the high heat capacity of oblate spheroids in some conditions (e.g., 4:5 aspect ratio, y = 3 or 4) can be explained by phase transitions, an overall increase combined with increased error is seen for all simulations in the glassy solid phase. While energy and dielectric response (fluctuation method) were well-converged for these simulations, it is likely that the number of configurations used was not adequate to sample over the large intermolecular phonon modes involving chains or clusters of spheroids in these strongly bound systems. Insufficient sampling could lead to an overestimation of the fluctuation magnitude and heat capacity; a reduction of calculated heat capacity was observed between the production simulations and shorter preliminary simulations (see the Supporting Information). In addition to the total energies, the ratio of electrostatic to Lennard-Jones (dispersion + steric repulsion) energy was examined. Results are shown in Figure 4.
Figure 3. The total potential energy (3a) and potential energy contribution to the constant-volume heat capacity (3b) of ensembles of dipolar spheroids at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field are shown as functions of the aspect ratio. The heat capacity was calculated from fluctuations in the total potential energy.48 Figure 4. Ratio of the electrostatic and Lennard-Jones energy (UES/ ULJ) in ensembles of dipolar spheroids at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field. The ratio varies little as a function of aspect ratio for prolate spheroids; a detail of this region is shown in the inset.
As with their dielectric properties, the thermodynamics of oblate spheroids differ greatly from spheroids with aspect ratios of unity or greater, particularly at high dipole densities. Spheres and prolate spheroids exhibit only small changes in their total potential energy or heat capacity as a function of aspect ratio. However, the total energy decreases dramatically with the aspect ratio for oblate spheroids, consistent with previously observed trends in the electrostatic energy with aspect ratio.20 As the reduction in aspect ratio allows for close, energetically favorable stacking of dipoles and as εLJ is held constant, the electrostatic term dominates this decrease in energy. The total energy also decreases with increasing dipole density, again due to stronger electrostatic interactions. The potential energy contribution to the heat capacity of the spheroids follows an inverse relationship with potential energy
At the lowest dipole density (y = 2), the Lennard-Jones contribution to the total potential energy is a factor of 2 larger than the electrostatic contribution at all of the aspect ratios examined and nearly constant as a function of aspect ratio. At all dipole densities with y > 2, the relationship between the electrostatic and Lennard-Jones contributions to the potential energy changes dramatically. At the most oblate aspect ratio (1:2), the electrostatic interactions dominate heavily at y ≥ 4. The strong electrostatic interactions pull the spheroids into F
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B close contact and cause the Lennard-Jones interactions to become energetically unfavorable (overwhelmed by the r−12 nuclear repulsion term). The cohesion of the system is maintained by head-to-tail stacking of the spheroids along their dipole axis. The ratio reaches a maximum for the 2:3 oblates at y = 6, with the favorable electrostatic contribution 75 times as large as the unfavorable Lennard-Jones contribution. The aspect ratio at which Lennard-Jones interactions become unfavorable increases as a function of y, indicating that stronger dipoles are pulling ellipsoids into a tighter packing arrangement. For all dipole densities, the ratio decays toward a small, constant value (0.37 to 1.83 for the dipole densities studied) as the aspect ratio increases, with all of the prolate spheroids exhibiting UES/ULJ ratios of less than 2. Consequently, the large increase in the magnitude of the total potential energy of oblates at high dipole densities is nearly entirely from dipole− dipole interactions. The formation of spontaneous order in spheroids with aspect ratios between 0.5 and 2 is strongly correlated with a high UES to ULJ ratio. Radial Distribution Function. Positional ordering was examined by means of the radial distribution function (eq 12) for five aspect ratios at y = 2 and y = 4; results appear in Figure 5. At y = 2, the most oblate (1:2) and prolate (2:1) spheroids show only weak positional correlation, which decays rapidly after two neighbor shells. Spheroids with intermediate aspect ratios show stronger positional correlation with uniform spacing between shells of spheroids and a longer persistence distance. However, these spheroids do not have any long-range positional order, as evinced by g(r) decaying to unity before the edge of the simulation cell is reached. The rapid decay of the radial distribution function is likely indicative of a system in the liquid phase. At y = 4, spheres and mild prolates showed similar behavior to y = 2. The 4:5 oblates showed weaker initial correlation, but the correlation persisted substantially longer than that for spheres or prolates. The 1:2 oblates, however, showed clear signs of short-range ordering, with narrow peaks and very strong correlation of the first neighboring shell and erratic correlation thereafter. The persistence length of the order was no longer than that of spheres or prolates and may be indicative of an amorphous (glassy) solid. Dipole Correlation Function. The persistence distance for dipole alignment was examined by means of the first-order dipole−dipole correlation function (eq 13) for the same conditions as the radial distribution function. Results appear in Figure 6. Oblate and prolate spheroids at y = 2 both exhibited an initial anticorrelated peak followed by a shell of aligned dipoles, with rapid decay of correlation beyond second-nearest neighbors and negligible correlation after 15 Å. This is indicative of centrosymmetric pairing of dipoles without long-range order. The initial anticorrelation is weakest for the 4:5 oblate spheroids; the peak is barely visible in Figure 6a. The Stockmayer fluid does not show any initial anticorrelation but exhibits positive initial correlation, indicative of short-range acentric dipole pairing but with little correlation of dipoles beyond nearest neighbors. The highest aspect ratio (2:1) spheroids exhibit the strongest short-range anticorrelation as the spheroids can approach each other more closely in a direction orthogonal to the dipole moment, favoring antiparallel stacking of dipoles.1 The positive correlation is strongest for the Stockmayer fluid. Even at this low dipole
Figure 5. Radial distribution functions (g(r)) at y = 2 (5a) and 4 (5b) at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field.
density, slight anisotropy leads to a qualitative difference in local ordering. At y = 4, the behavior of the prolate spheroids is relatively unchanged compared to that at y = 2, but the behavior of spheres and oblates is greatly altered. Instead of oscillating around zero, the curve for the Stockmayer fluid remains positive until it decays to zero, indicating net local correlation of dipoles (formation of ferroelectric domains) but not pervasive ferroelectric behavior.1 Similar local correlation of dipoles is seen for the 1:2 oblates, albeit to a greater extent and in the form of sharp peaks instead of a smooth, polarized region. The 4:5 oblates, however, are strongly ferroelectric, with net dipole correlation throughout the entire simulation cell. Phase Behavior. On the basis of dielectric, thermodynamic, and structural data, simulations at each aspect ratio and dipole density were assigned to one of the four phases21,25 shown in Table 3; representative visualizations of the four phases listed above are shown in Figure 7. In comparing Figure 7b and 7c, it can be seen that the ferroelectric columnar phase, which only appeared for a single set of conditions (4:5 oblates at y = 4), is distinct from the ferroelectric liquid phase and is consistent with prior work26 G
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
chain. The chains themselves were randomly ordered such that the system has little or no net order. Both the polymerization and reduced ordering/dielectric constant in oblate systems that formed these sorts of chains was consistent with prior simulation work by Rommel and Robinson,4 and such chains have also been observed in the Stockmayer fluid under a variety of temperatures and densities.56,57 The large net dipole moment of each chain causes its neighbors to order themselves antiparallel to it, hindering the development of net acentric order. Figure 8 shows the results from all 45 sets of simulations, combined into a phase diagram defined by their dipole density and aspect ratio. At the temperature and density examined, the ferroelectric phase was found to exist within a narrow band of aspect ratios, with strong acentric ordering occurring only for spheres and 4:5 oblates and only for y ≥ 3. The Stockmayer fluid formed pervasively ferroelectric phases, with domain sizes larger than the simulation cell size. The equilibrium domain size was difficult to determine for 4:5 oblates at the highest dipole densities due to a very shallow potential energy surface, though these systems demonstrated very strong net polarization. Oblates with a 2:3 aspect ratio at y = 3 exhibited ferroelectric domains but weaker net polarization and were classified as being borderline between a ferroelectric and glassy phase. Spheroids more oblate than those in the narrow ferroelectric band formed glassy polymer-like materials, and all of the prolate spheroids examined remained in the isotropic phase even at a 5:4 aspect ratio and y = 6. While ferroelectric phases have been observed for prolates in calculations by Groh and Dietrich,3 they formed at packing fractions and/or aspect ratios higher than those considered in this work. The highly ordered, ferroelectric columnar phase was confined to a single dipole density/aspect ratio combination and developed late in the simulations, which passed through a metastable ferroelectric nematic phase before converging to even higher order (see the Supporting Information). Similar behavior was not observed for surrounding dipole densities and aspect ratios, indicating that this phase is stable only within a narrow range of electrostatic and thermodynamic conditions.
Figure 6. Dipole correlation functions (gμ(r)) at y = 2 (6a) and 4 (6b) at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field.
■
indicating the possibility of ferroelectric phases with positional order (e.g., ferroelectric nematic) for the Stockmayer fluid. The ferroelectric columnar phase exhibits nearly one-dimensional order with particles consistently aligned along a single director axis even in the absence of a field, whereas the ordinary ferroelectric phase exhibits net polarization without strong centrosymmetric alignment. The ferroelectric columnar phase also exhibits net positional correlation between particles into long, regularly aligned chains, though without the layering that would be characteristic of a smectic phase. Simulations in the glassy phase exhibited vermiculated structures with the molecular dipoles aligned along each
CONCLUSION
Subtle changes in the aspect ratio of dipolar spheroids can have dramatic effects on their ordering behavior. At the lowest dipole density (y = 2), all of the spheroids examined remained in the isotropic phase. At higher dipole densities, ferroelectric behavior was observed in a very narrow range of aspect ratios close to unity. None of the prolate spheroids studied at ρ* = 1 and 293 K, even those with a 5:4 aspect ratio, exhibited ferroelectric behavior, indicating that even the slightest anisotropy is enough to disrupt dipole correlation and prevent
Table 3. Phase Characteristics for Dipolar Spheroids isotropic ε/εOnsager UES/ULJ P1 P2 g(r) gμ(r)
1−5 small, often ≤1 0 0 smooth, short persistence smooth, oscillates +/−, short persistence
ferroelectric liquid
ferroelectric columnar
>5 small but >1 >0.1 >0.1 smooth, short persistence always positive, long persistence H
large >1 >0.7 >0.7 smooth, long persistence always positive, infinite persistence
glassy solid
Dielectric and Phase Behavior of Dipolar Spheroids Lewis E. Johnson,* Stephanie J. Benight, Robin Barnes, and Bruce H. Robinson* Department of Chemistry, University of Washington, Seattle, Washington 98195, United States S Supporting Information *
ABSTRACT: The Stockmayer fluid, composed of dipolar spheres, has a well-known isotropic−ferroelectric phase transition at high dipole densities. However, there has been little investigation of the ferroelectric transition in nearly spherical fluids at dipole densities corresponding to those found in many polar solvents and in guest−host organic electro-optic materials. In this work, we examine the transition to ordered phases of low-aspect-ratio spheroids under both unperturbed and poled conditions, characterizing both the static dielectric response and thermodynamic properties of spheroidal systems. Spontaneous ferroelectric ordering was confined to a small region of aspect ratios about unity, indicating that subtle changes in sterics can have substantial influence on the behavior of coarse-grained liquid models. Our results demonstrate the importance of molecular shape in obtaining even qualitatively correct dielectric responses and provide an explanation for the success of the Onsager model as a phenomenological representation for the dielectric behavior of polar organic liquids.
■
INTRODUCTION Coarse-grained spheroidal models of molecules approximate the steric behavior of a wide variety of molecules, including small organic liquids,1−3 organic electro-optic chromophores,4−7 and liquid crystals (LCs)8 at low computational cost. Such models, typically implemented using the Gay−Berne potential,9 variants thereof,10,11 the Berne−Pechukas potential,12 or the method of Perram and Wertheim,13 introduce a directional dependence to the well depth (εLJ) and/or van der Waals diameter (σLJ) of the ubiquitous Lennard-Jones potential, breaking its spherical symmetry. The steric anisotropy can result in dramatically different behavior in comparison with purely spherical models such as the Stockmayer fluid14 (sphere containing a point dipole at the center). For uniaxial LCs, the molecular aspect ratio (length/width) is particularly crucial. Onsager demonstrated that smectic or nematic phases could be induced in rigid spherocylinders (cylinders with hemispherical ends) by increasing their aspect ratio in his pioneering 1949 paper15 on LC phases of colloidal materials. It was later shown that LC phases could also be formed by prolate spheroids16 and that steric anisotropy was a necessary condition but not a sufficient condition for forming LC phases.17,18 Sufficiently oblate (discotic) spheroids with axial dipoles were demonstrated to form LC phases by Frenkel,19 Patey,20,21 Zannoni,22 and others.23 Oblate spheroids with transverse dipoles24 have also been demonstrated. Further work by Baus,25 Groh,3,26 and Ayton21 demonstrated that dipolar spheroids could form ferroelectric LC phases. In addition, work by Berne2,8 and co-workers explored the effects of electronic polarizability on simple models of LC systems. Hall and co-workers have also examined the phase behavior of dipolar Janus particles.27 Development of simulation methods for these types of systems is crucial for the development of © XXXX American Chemical Society
coarse-grained protocols for simulating acentric ordering in electro-optic chromophores and other highly polar molecular systems. In addition to changes in centrosymmetric (LC) and acentric (ferroelectric) ordering at equilibrium, changes in the aspect ratio of a molecule can have dramatic effects on the induced acentric ordering of dipole moments with respect to an applied electric field. Both analytic models and Monte Carlo simulations have demonstrated that dipolar spheres tend to align more strongly with an electric field than prolate spheroids.4,28 The dramatic difference in ordering behavior has had significant implications for the design of chromophores for electro-optic applications; such chromophores are typically highly prolate due to their long π-systems.29 In electro-optics, the ordering is typically quantified by the acentric order parameter ⟨cos3 θ⟩, where θ is the angle between the molecular dipole moment and the poling field used to align the chromophores. This cubic average cannot currently be measured directly in amorphous materials and is generally estimated through simulation or measured indirectly based on electro-optic response30 or second harmonic generation.31 Acentric order may also be indirectly quantified through capacitance measurements32 of the dielectric constant ε, where
(ε − 1)E = 4π P
(1)
Here, E is the macroscopic electric field applied to the material, and P is the total polarization of the material. We further apply four simplifying conditions: (1) ε is isotropic (scalar). Received: January 1, 2015 Revised: March 23, 2015
A
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B (ε − n2)(2ε + n2) =y ε(n2 + 2)
(2) The material is composed of identical molecules with permanent dipole moment μ. (3) All polarization of the system is due to reorientation of permanent dipoles as opposed to formation of induced dipoles (negligible electronic polarizability). (4) The system is oriented such that the electric field is aligned along the z-axis of the lab frame. Under these conditions, eq 1 simplifies to ε−1=
4πρN μ⟨cos θ ⟩ 4πMz = VEz Ez
where n is the refractive index of the liquid, was developed to correct for the Debye37 model’s incorrect prediction of a ferroelectric transition at y = 1. However, while the Onsager model is capable of reasonably approximating the experimental relationship41 between the dipole moment and dielectric constant in liquids, it is qualitatively incorrect in modeling the many-body behavior of the system it was derived for (dipolar spheres) compared to results obtained from numerical simulations.14,42 The Onsager model does not predict a ferroelectric phase transition under any conditions. However, unlike the Stockmayer fluid, most real polar liquids exhibit little intermolecular correlation.36 The low intermolecular correlation in real liquids results in good phenomenological agreement between the Onsager model and real liquids, despite the Onsager model’s underestimation of correlation in the Stockmayer fluid. Here, we present work that examines the ordering behavior and degree of intermolecular correlation (as quantified by deviation from the Onsager model) for a range of dipolar spheroids. While previous work has often focused on wide ranges of aspect ratios (up to several orders of magnitude)17 and/or temperatures22 and densities,19 this work will focus on the effects of dipole moment and aspect ratio at constant temperature, density, and molecular volume. We consider a small range of different aspect ratios centered about unity and determine the degree of ordering of dipolar spheroids in a fluid compared to that of the Stockmayer fluid.
(2)
in the static (DC) limit. Here, V is the volume of the system, ρN is the number density of the dipoles, and θ is the angle between any individual dipole and the applied field, such that ⟨cos θ⟩ is the cosine of the angle between the total dipole moment M and the applied field. Within the assumptions of eq 2, ⟨cos θ⟩ is simply the normalized projection of the dipole moment along the z-axis. The third-order average ⟨cos3 θ⟩ (crucial to electrooptic materials) tends to closely track ⟨cos θ⟩.5 Furthermore, if the system is composed of spherical molecules, its thermodynamic state can be uniquely represented as a function of reduced dipole moment and reduced density14 μ*2 =
μ2 σLJ3kBT
ρ* = ρN σLJ3
(3)
Here, T is the absolute temperature, and σLJ is the LennardJones diameter of the molecule. These reduced parameters can be generalized20 to nonspherical molecules by replacing the single σLJ with the molecular semiaxes (a, b, and c), such that ρspheroid * = ρNabc. Additionally, the dipole density of any polar material can be reduced to a single dimensionless parameter33
y=
■
COMPUTATIONAL METHODS
Simulations of dipolar spheroids were run at nine different aspect ratios and five different dipole moments, holding all other thermodynamic parameters and terms in the intermolecular potential constant. Maximum rotational and translational move sizes were optimized for each system in order to maximize the RMS displacement per Monte Carlo step. Methods were substantially identical to those used in our 2010 paper.1 All simulations used self-consistent reaction field (SCRF) electrostatic boundary conditions. Run lengths were extended compared to our prior work1 to improve convergence of second-order properties such as fluctuations in the dipole moment and fluctuations in the energy.43 The simulation size was set at N = 1728 particles, and simulations were initially run for a length of 120 000 cycles, with N trial moves per cycle for a total of 2.07 × 108 configurations per simulation. The final 60 000 cycles of each simulation were used for calculating properties. Oblate spheroids with y ≥ 3 exhibited slow convergence, and the simulation length for these systems was extended to 300 000 cycles, using the last 100 000 cycles to calculate properties. The run length was further extended if simulations were still not converged to an energy difference of i ⎢ ⎠ ⎩ ⎣ ⎝
calculated self-consistently based on the most recent reaction field update. The dielectric constant was calculated using eq 2 in the calculations with an external field. In calculations without an external field, it was calculated using a Kirkwood fluctuation expansion33,44
⎤ 1 1 2(εS − 1) ⎥ μi ·μj ⎥ + 3 (μi ·μj − 3(μi · riĵ )(riĵ ·μj )) − 3 ε 2 1 + rij R cut S ⎥⎦ ⎫ ⎪ 2εS + ε μi · E⎬ − 2εS + 1 ⎪ ⎭
⟨Mz⟩E0,z =
Here, εLJ is the Lennard-Jones energy (set to 217 K), σeff is the effective Lennard-Jones contact distance calculated using the method of Perram and Wertheim,13 Rcut is the cutoff distance for calculating electrostatic interactions (27.96 Å), εS is the reaction field dielectric constant, which is self-consistently updated every cycle,1 and E is a homogeneous external field applied along the z-axis of the simulation cell. Lennard-Jones interactions were truncated at a cutoff distance of 18.64 Å, corresponding to 4σ, for reasons of computational efficiency. Each of the nine spheroid geometries was run at all five dipole densities and a constant number density of ρ* = 1, intended to be representative of common organic liquids. The NVT ensemble was chosen to ensure constant dipole densities at all aspect ratios. The thermodynamic parameters μ* and ρ* were redefined to reflect the spheroidal nature of the systems as μ*2 = μ2/(σaσ2b)kBT and ρ* = ρN(σaσ2b), where σa is the length of the dipole axis of the spheroid and σb is the length of the nondipole axes. The electrostatic parameters for each y are in Table 1, and the spheroid geometries for the different aspect ratios are listed in Table 2.
T (K)
μ*2
ρ*
2 3 4 5 6
2.42 2.96 3.42 3.82 4.19
293 293 293 293 293
1.43 2.14 2.86 3.57 4.29
1.0 1.0 1.0 1.0 1.0
A (R ) =
σa (Å)
aspect ratio
V (Å3)
1:2 4:7 2:3 4:5 1:1 5:4 3:2 7:4 2:1
5.88 5.62 5.34 5.02 4.66 4.32 4.08 3.86 3.70
2.94 3.20 3.56 4.02 4.66 5.40 6.10 6.76 7.40
0.50 0.57 0.67 0.80 1.00 1.20 1.50 1.75 2.00
52.98 52.98 52.98 52.98 52.98 52.99 52.98 52.98 52.98
(9)
1 2
∫0
∞
dx (x + 1) (R2x + 1) 3/2
(10)
This approximation adds the electrostatic effects of molecular anisotropy (distortion of the dielectric cavity around the molecule) in the absence of intermolecular correlations. In addition to ordering resulting from polarization in response to an applied field, spontaneous order was calculated from simulations at zero field. As simulations at zero field do not have an external director axis that is known a priori, ordering is quantified by means of the de Gennes46 Q-tensor
Table 2. Simulation Parameters: Spheroid Geometries σb (Å)
(8)
where a is the major (dipole) semiaxis of the molecule and b is the minor semiaxis, such that the aspect ratio is R= a/b and
Q= eccentricity
(⟨M2⟩0 − ⟨M2⟩20 )
(ε − 1)(ε + A(a , b)(1 − ε)) =y 3ε
Table 1. Simulation Parameters: Dipole Strength μ (D)
3kBT
of the total dipole moment M of the system at zero applied field, assuming a linear response of M to E0 and substituting into eq 2 to obtain the dielectric response. The two types of dielectric calculations were then compared to determine whether or not the linear response approximation was valid under each condition. The Onsager model (eq 5) was used as an uncorrelated control model. Onsager’s predictions were calculated at each dipole density, neglecting electronic polarizability as per the assumptions of eq 2 such that the refractive index was n = 1. Additionally, the Onsager model can be generalized to spheroids (as initially suggested in Onsager’s derivation38) by means of a Dirichlet correction45 on the local field factor such that
(6)
y
E0, z
3⟨μi ⊗ μj ⟩ − I (11)
2
averaged over all μi ≠ μj. The dominant eigenvalue of Q is the nematic order parameter P2 with respect to the internal director axis d of the system, which is itself the dominant eigenvector of Q.34 The corresponding acentric order parameter P1 = M̂ ·d can be obtained by projection of the total dipole moment of the system on the director and is the zero-field analogue to ⟨cos θ⟩. Finally, the radial distribution function33,47,48 N
g (r ) =
N
V ∑i ∑ j > i nij(r , r + Δr ) 2πN (N − 1)r 2Δr
(12)
and first-order dipole−dipole (angular) correlation function1,49 N
In order to obtain both fluctuation-based and direct measurements of the dielectric constant, simulations were run under both zero applied field and an applied field of E = 50 V/ μm, with the local field E=
2εS + ε E 2εS + 1
gμ(r ) =
N
V ∑i = 1 ∑ j > i cos γij(r , r + Δr ) 2πN (N − 1)r 2Δr
(13)
were calculated, where nij(r,r+Δr) is the number of particles where the center−center distance rij is between r and r + Δr and cos γ = μ̂ i·μ̂ j. The former indicates the degree of positional correlation between atoms, with the function decaying to unity in a disordered structure such as an isotropic liquid or a glass
(7) C
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
Figure 1. Static dielectric constant of spheroids of different aspect ratios (geometries shown in (1a)), calculated using both dipole fluctuations (shown in (1b)) and in response to an applied external field of 50 V/μm (shown in (1c)). Due to a large range of values, the ordinate is plotted on a logarithmic scale. Calculations were run at 293 K and ρ* = 1 using SCRF boundary conditions at dipole densities in the range of y = 2−6. Predictions from the Onsager model (eq 9) are shown as solid lines. Dotted lines denote the maximum possible dielectric response (perfect order) under the applied field, and dashed lines are guides to the eye.
results, a table of dielectric constants is also available as Supporting Information. Results at the lowest dipole density (y = 2) were within the linear response regime, with strong agreement between the applied field and fluctuation results at all aspect ratios. The dielectric constant depended strongly on the aspect ratio, forming a Lorentzian-like function peaking at an aspect ratio of unity (the Stockmayer fluid). The peak value was almost four times larger than the lowest (2:1 prolate) result, consistent with literature values for similar systems.42 Spheres and nearly spherical systems were poorly described by the Onsager model (eq 9, shown as solid curves on the graphs). However, the more prolate spheroids were reasonably well described by the Onsager model, indicating a very low degree of net intermolecular correlation.
and exhibiting repeating peaks in a system with strong positional order such as a smectic LC. Equation 13 indicates the degree of angular correlation between dipoles, with negative values indicating net anticorrelation at a given distance and positive values indicating net correlation. In an isotropic material, the angular correlation will decay to zero (no correlation), but in a ferroelectric material, it will limit to a positive value.21,50
■
RESULTS AND DISCUSSION
Static Dielectric Constant. The static dielectric constant was calculated using both the applied field and fluctuation method at six dipole densities, ranging from y = 2 (lowest) to 6 (highest); results are shown in Figure 1 along with the geometries of the dipolar spheroids; due to the large range of D
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B While results for prolate spheroids at y = 3 followed similar trends to those of y = 2, the behavior of oblates differed dramatically. Mildly oblate (4:5) spheroids transitioned to a ferroelectric phase, and 2:3 spheroids showed evidence of ferroelectric domains. Furthermore, the linear response approximation in eq 8 began to break down, with substantial disagreement between applied field and fluctuation calculations for oblate (a < b) spheroids. While this could be due to slow oscillation of ⟨M2⟩ and/or difficulty in averaging out any net dipole moment,43,51,52 results were consistent among the trajectories, and calculated dielectric constants became smaller as simulation lengths increased. Prolate spheroids with an aspect ratio of 3:2 or higher were still well-approximated by the Onsager model. Behavior at y = 4 was qualitatively similar to that at y = 3 at all aspect ratios. However, the ferroelectric domains at the 2:3 aspect ratio disappeared, and the ferroelectric response at the 4:5 aspect ratio was enhanced. The prolate spheroids exhibited much lower dielectric constants than the Stockmayer fluid, and agreement between the fluctuation and applied field methods remained strong for all a > b. The range of aspect ratios for which the Onsager model was approximately valid decreased. By y = 5 and 6, the Stockmayer fluid exhibited a very high dielectric constant, consistent with our prior results.1 Trends were otherwise similar to those of y = 4. The 1:2 oblates at the highest dipole density showed a weak increase in dielectric response to an applied field, although this was not reflected in the dielectric response calculated from dipole fluctuations. Interestingly, even at the highest dipole densities, the applied field and fluctuation-based responses for prolates closely mirrored each other. The Onsager model remained a reasonable approximation for the 7:4 and 2:1 prolates. The close match between the Onsager model and Monte Carlo results for prolate spheroids over the wide range of dipole densities raises the possibility that the Onsager model’s ability to closely approximate experimental data for polar liquids is due to a fortuitous cancellation of error; nonspherical systems (especially prolate spheroids) exhibit far weaker intermolecular correlation than spheres, such that the Kirkwood g factor33,44 remains near unity for many realistic molecular geometries. However, the purely electrostatic contribution from the shape of the Onsager cavity around the molecule, which is accounted for in the extended Onsager model (eq 9), depends only weakly on the cavity shape. Two of the most polar small organic liquids (acetonitrile and nitromethane) are highly prolate, and others such as acetone and dimethylformamide are oblate but with an in-plane dipole instead of an axial one. As such, ordering of real molecules in response to a field is frustrated by intermolecular correlation1,39 due to molecular shape. The relatively small difference between results calculated with eqs 5 and 9 (see the Supporting Information) explains why the Onsager model, although constructed using spheres, is more effective at explaining dielectric trends in real liquids than as a standalone theory of simple liquids. Spontaneous Order. The order parameters ⟨P1⟩ and ⟨P2⟩ were calculated using eq 11 for the simulations at each dipole density in the absence of an applied external field. These results are shown in Figure 2. No spontaneous ordering was observed in prolate ellipsoids at any dipole density under the conditions studied. In general, prolates that form nematic phases at ρ* = 1 often have substantially larger aspect ratios than those considered here.3,19,50,53 Furthermore, no spontaneous ordering (centro-
Figure 2. Spontaneous acentric (2a) and centrosymmetric (2b) ordering as a function of aspect ratio for dipolar spheroids at six different dipole densities at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field, as obtained from the ordering matrix Q (eq 11).
symmetric or acentric) was observed for any aspect ratio at y = 2. At the other dipole densities, the more oblate systems with aspect ratios between 2:3 and 1:2 exhibited weak spontaneous ordering, with ⟨P2⟩ typically larger than ⟨P1⟩. These phases may have slight antiferroelectric character. Strong (⟨P1⟩ > 0.3) acentric ordering was observed for 4:5 oblates at y = 3, 4, and 5 and for the Stockmayer fluid at y = 5 and 6, consistent with the singularities in the dielectric constant seen at those dipole densities (see Figure 1). This strong ordering in the absence of a field is indicative of formation of a ferroelectric20,54 phase and is most pronounced for 4:5 oblates at y = 4, where ⟨P1⟩ approaches unity. For all of these ferroelectric systems, ⟨P1⟩ was larger than ⟨P2⟩. Only two conditions (4:5 oblates at y = 4 or 5) showed a degree of centrosymmetric ordering similar to the degree of acentric ordering observed. Strong centrosymmetric order in the absence of acentric order (pure nematic or smectic phases as observed in typical LCs) was not observed over the range of aspect ratios and dipole densities studied, though 2:3 oblates at y = 3 exhibited extensive centrosymmetric order with only a weak ferroelectric response. Thermodynamics. In addition to ordering behavior, the total potential energy (U) and the potential energy E
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B contribution55 to the constant-volume heat capacity (Cv) were examined for signs of phase transitions; results appear in Figure 3. Heat capacities were calculated from fluctuations in the total potential energy at constant temperature.
and is minimized for the systems with the weakest intermolecular interactions. The heat capacity varies weakly as a function of aspect ratio at all dipole densities for prolate spheroids and spheres, reaching a weak local minimum at an aspect ratio of 5:4. The trend in heat capacity for prolates and spheres is nearly independent of dipole density, and the change in heat capacity across the isotropic−ferroelectric phase transition is modest. A substantial shift in heat capacity (Cv) is also seen as the aspect ratio is decreased below unity. The heat capacity increases rapidly as a function of aspect ratio for oblates with y > 2. The highest heat capacity is observed for 1:2 oblates at the highest dipole densities. While the high heat capacity of oblate spheroids in some conditions (e.g., 4:5 aspect ratio, y = 3 or 4) can be explained by phase transitions, an overall increase combined with increased error is seen for all simulations in the glassy solid phase. While energy and dielectric response (fluctuation method) were well-converged for these simulations, it is likely that the number of configurations used was not adequate to sample over the large intermolecular phonon modes involving chains or clusters of spheroids in these strongly bound systems. Insufficient sampling could lead to an overestimation of the fluctuation magnitude and heat capacity; a reduction of calculated heat capacity was observed between the production simulations and shorter preliminary simulations (see the Supporting Information). In addition to the total energies, the ratio of electrostatic to Lennard-Jones (dispersion + steric repulsion) energy was examined. Results are shown in Figure 4.
Figure 3. The total potential energy (3a) and potential energy contribution to the constant-volume heat capacity (3b) of ensembles of dipolar spheroids at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field are shown as functions of the aspect ratio. The heat capacity was calculated from fluctuations in the total potential energy.48 Figure 4. Ratio of the electrostatic and Lennard-Jones energy (UES/ ULJ) in ensembles of dipolar spheroids at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field. The ratio varies little as a function of aspect ratio for prolate spheroids; a detail of this region is shown in the inset.
As with their dielectric properties, the thermodynamics of oblate spheroids differ greatly from spheroids with aspect ratios of unity or greater, particularly at high dipole densities. Spheres and prolate spheroids exhibit only small changes in their total potential energy or heat capacity as a function of aspect ratio. However, the total energy decreases dramatically with the aspect ratio for oblate spheroids, consistent with previously observed trends in the electrostatic energy with aspect ratio.20 As the reduction in aspect ratio allows for close, energetically favorable stacking of dipoles and as εLJ is held constant, the electrostatic term dominates this decrease in energy. The total energy also decreases with increasing dipole density, again due to stronger electrostatic interactions. The potential energy contribution to the heat capacity of the spheroids follows an inverse relationship with potential energy
At the lowest dipole density (y = 2), the Lennard-Jones contribution to the total potential energy is a factor of 2 larger than the electrostatic contribution at all of the aspect ratios examined and nearly constant as a function of aspect ratio. At all dipole densities with y > 2, the relationship between the electrostatic and Lennard-Jones contributions to the potential energy changes dramatically. At the most oblate aspect ratio (1:2), the electrostatic interactions dominate heavily at y ≥ 4. The strong electrostatic interactions pull the spheroids into F
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B close contact and cause the Lennard-Jones interactions to become energetically unfavorable (overwhelmed by the r−12 nuclear repulsion term). The cohesion of the system is maintained by head-to-tail stacking of the spheroids along their dipole axis. The ratio reaches a maximum for the 2:3 oblates at y = 6, with the favorable electrostatic contribution 75 times as large as the unfavorable Lennard-Jones contribution. The aspect ratio at which Lennard-Jones interactions become unfavorable increases as a function of y, indicating that stronger dipoles are pulling ellipsoids into a tighter packing arrangement. For all dipole densities, the ratio decays toward a small, constant value (0.37 to 1.83 for the dipole densities studied) as the aspect ratio increases, with all of the prolate spheroids exhibiting UES/ULJ ratios of less than 2. Consequently, the large increase in the magnitude of the total potential energy of oblates at high dipole densities is nearly entirely from dipole− dipole interactions. The formation of spontaneous order in spheroids with aspect ratios between 0.5 and 2 is strongly correlated with a high UES to ULJ ratio. Radial Distribution Function. Positional ordering was examined by means of the radial distribution function (eq 12) for five aspect ratios at y = 2 and y = 4; results appear in Figure 5. At y = 2, the most oblate (1:2) and prolate (2:1) spheroids show only weak positional correlation, which decays rapidly after two neighbor shells. Spheroids with intermediate aspect ratios show stronger positional correlation with uniform spacing between shells of spheroids and a longer persistence distance. However, these spheroids do not have any long-range positional order, as evinced by g(r) decaying to unity before the edge of the simulation cell is reached. The rapid decay of the radial distribution function is likely indicative of a system in the liquid phase. At y = 4, spheres and mild prolates showed similar behavior to y = 2. The 4:5 oblates showed weaker initial correlation, but the correlation persisted substantially longer than that for spheres or prolates. The 1:2 oblates, however, showed clear signs of short-range ordering, with narrow peaks and very strong correlation of the first neighboring shell and erratic correlation thereafter. The persistence length of the order was no longer than that of spheres or prolates and may be indicative of an amorphous (glassy) solid. Dipole Correlation Function. The persistence distance for dipole alignment was examined by means of the first-order dipole−dipole correlation function (eq 13) for the same conditions as the radial distribution function. Results appear in Figure 6. Oblate and prolate spheroids at y = 2 both exhibited an initial anticorrelated peak followed by a shell of aligned dipoles, with rapid decay of correlation beyond second-nearest neighbors and negligible correlation after 15 Å. This is indicative of centrosymmetric pairing of dipoles without long-range order. The initial anticorrelation is weakest for the 4:5 oblate spheroids; the peak is barely visible in Figure 6a. The Stockmayer fluid does not show any initial anticorrelation but exhibits positive initial correlation, indicative of short-range acentric dipole pairing but with little correlation of dipoles beyond nearest neighbors. The highest aspect ratio (2:1) spheroids exhibit the strongest short-range anticorrelation as the spheroids can approach each other more closely in a direction orthogonal to the dipole moment, favoring antiparallel stacking of dipoles.1 The positive correlation is strongest for the Stockmayer fluid. Even at this low dipole
Figure 5. Radial distribution functions (g(r)) at y = 2 (5a) and 4 (5b) at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field.
density, slight anisotropy leads to a qualitative difference in local ordering. At y = 4, the behavior of the prolate spheroids is relatively unchanged compared to that at y = 2, but the behavior of spheres and oblates is greatly altered. Instead of oscillating around zero, the curve for the Stockmayer fluid remains positive until it decays to zero, indicating net local correlation of dipoles (formation of ferroelectric domains) but not pervasive ferroelectric behavior.1 Similar local correlation of dipoles is seen for the 1:2 oblates, albeit to a greater extent and in the form of sharp peaks instead of a smooth, polarized region. The 4:5 oblates, however, are strongly ferroelectric, with net dipole correlation throughout the entire simulation cell. Phase Behavior. On the basis of dielectric, thermodynamic, and structural data, simulations at each aspect ratio and dipole density were assigned to one of the four phases21,25 shown in Table 3; representative visualizations of the four phases listed above are shown in Figure 7. In comparing Figure 7b and 7c, it can be seen that the ferroelectric columnar phase, which only appeared for a single set of conditions (4:5 oblates at y = 4), is distinct from the ferroelectric liquid phase and is consistent with prior work26 G
DOI: 10.1021/acs.jpcb.5b00009 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
chain. The chains themselves were randomly ordered such that the system has little or no net order. Both the polymerization and reduced ordering/dielectric constant in oblate systems that formed these sorts of chains was consistent with prior simulation work by Rommel and Robinson,4 and such chains have also been observed in the Stockmayer fluid under a variety of temperatures and densities.56,57 The large net dipole moment of each chain causes its neighbors to order themselves antiparallel to it, hindering the development of net acentric order. Figure 8 shows the results from all 45 sets of simulations, combined into a phase diagram defined by their dipole density and aspect ratio. At the temperature and density examined, the ferroelectric phase was found to exist within a narrow band of aspect ratios, with strong acentric ordering occurring only for spheres and 4:5 oblates and only for y ≥ 3. The Stockmayer fluid formed pervasively ferroelectric phases, with domain sizes larger than the simulation cell size. The equilibrium domain size was difficult to determine for 4:5 oblates at the highest dipole densities due to a very shallow potential energy surface, though these systems demonstrated very strong net polarization. Oblates with a 2:3 aspect ratio at y = 3 exhibited ferroelectric domains but weaker net polarization and were classified as being borderline between a ferroelectric and glassy phase. Spheroids more oblate than those in the narrow ferroelectric band formed glassy polymer-like materials, and all of the prolate spheroids examined remained in the isotropic phase even at a 5:4 aspect ratio and y = 6. While ferroelectric phases have been observed for prolates in calculations by Groh and Dietrich,3 they formed at packing fractions and/or aspect ratios higher than those considered in this work. The highly ordered, ferroelectric columnar phase was confined to a single dipole density/aspect ratio combination and developed late in the simulations, which passed through a metastable ferroelectric nematic phase before converging to even higher order (see the Supporting Information). Similar behavior was not observed for surrounding dipole densities and aspect ratios, indicating that this phase is stable only within a narrow range of electrostatic and thermodynamic conditions.
Figure 6. Dipole correlation functions (gμ(r)) at y = 2 (6a) and 4 (6b) at 293 K and ρ* = 1 under SCRF boundary conditions in the absence of an applied field.
■
indicating the possibility of ferroelectric phases with positional order (e.g., ferroelectric nematic) for the Stockmayer fluid. The ferroelectric columnar phase exhibits nearly one-dimensional order with particles consistently aligned along a single director axis even in the absence of a field, whereas the ordinary ferroelectric phase exhibits net polarization without strong centrosymmetric alignment. The ferroelectric columnar phase also exhibits net positional correlation between particles into long, regularly aligned chains, though without the layering that would be characteristic of a smectic phase. Simulations in the glassy phase exhibited vermiculated structures with the molecular dipoles aligned along each
CONCLUSION
Subtle changes in the aspect ratio of dipolar spheroids can have dramatic effects on their ordering behavior. At the lowest dipole density (y = 2), all of the spheroids examined remained in the isotropic phase. At higher dipole densities, ferroelectric behavior was observed in a very narrow range of aspect ratios close to unity. None of the prolate spheroids studied at ρ* = 1 and 293 K, even those with a 5:4 aspect ratio, exhibited ferroelectric behavior, indicating that even the slightest anisotropy is enough to disrupt dipole correlation and prevent
Table 3. Phase Characteristics for Dipolar Spheroids isotropic ε/εOnsager UES/ULJ P1 P2 g(r) gμ(r)
1−5 small, often ≤1 0 0 smooth, short persistence smooth, oscillates +/−, short persistence
ferroelectric liquid
ferroelectric columnar
>5 small but >1 >0.1 >0.1 smooth, short persistence always positive, long persistence H
large >1 >0.7 >0.7 smooth, long persistence always positive, infinite persistence
glassy solid
