# Ferroelectric glass of spheroidal dipoles with impurities: Polar nanoregions, response to applied electric field, and ergodicity breakdown

## Abstract

Using molecular dynamics simulation, we study dipolar glass in crystals composed of slightly spheroidal, polar particles and spherical, apolar impurities between metal walls. We present physical pictures of ferroelectric glass, which have been observed in relaxors, mixed crystals (such as KCNKBr), and polymers. Our systems undergo a diffuse transition in a wide temperature range, where we visualize polar nanoregions (PNRs) surrounded by impurities. In our simulation, the impurities form clusters and their space distribution is heterogeneous. The polarization fluctuations are enhanced at relatively high depending on the size of the dipole moment. They then form frozen PNRs as is further lowered into the nonergodic regime. As a result, the dielectric permittivity exhibits the characteristic features of relaxor ferroelectrics. We also examine nonlinear response to cyclic applied electric field and nonergodic response to cyclic temperature changes (ZFCFC), where the polarization and the strain change collectively and heterogeneously. We also study antiferroelectric glass arising from molecular shape asymmetry. We use an Ewald scheme of calculating the dipolar interaction in applied electric field.

## I Introduction

Ferroelectric transitions have been attracting much attention in various systems. It is known that they can occur even in simple particle systems. For example, one-component spherical particles with a point dipole undergo a ferroelectric transition in crystal or liquid-crystal phases if the dipole interaction is sufficiently strong (5); (1); (2); (6); (3); (4); (9); (10); (8); (7). Such spherical dipoles form various noncubic crystals in ferroelectric phases (7); (8). Ferroelectriciity was also studied in positionally disordered dipolar solids(4). Recently, Johnson et al.(11); (12) have investigated a ferroelectric transition of spheroidal particles with a dipole moment parallel to the spheroidal axis. They found that the static dielectric constant increases up to with increasing if the aspect ratio is close to unity. In this paper, we examine ferroelectric transitions in mixtures of slightly spheroidal dipoles and spherical impurities.

In many solids, the polarization is induced by ion displacements within unit cells and the dielectric constant is very large. As a unique aspect, the ferroelectric transitions become diffuse with a sufficient amount of disorder(15); (13); (14), which take place over a wide temperature range without long-range dipolar order. Notable examples are relaxors (16); (19); (20); (22); (21); (13); (17); (23); (18) such as Pb(MgNb)O (PMN), where the random distribution of Mg and Nb at B sites yields quenched random fields at Pb sites (24); (25). In relaxors, temperature-dependence of the optic index of refraction suggested appearance of mesoscopic polarization heterogeneities (26), called polar nanoregions (PNRs). They are enhanced at relatively high as near-critical fluctuations and are frozen at lower . It is widely believed that these PNRs give rise to a broad peak in the dielectric permittivity as a function of (24); (16); (19); (20); (13); (17); (27); (21); (22); (23); (18); (28). They have been detected by neutron and x-ray scattering (33); (32); (30); (29); (31); (23) and visualized by transmission electron microscopy (35); (37); (36) and piezoresponse force microscopy (39); (40); (38). Strong correlations have also been found between the PNRs and the compositional heterogeneity of the B site cations (33); (34); (35); (24); (40); (37); (43); (44); (41); (42).

Relaxor behaviors also appear in other disordered dipolar systems(13); (14); (15). In particular, orientational glass has long been studied in mixed crystals such as KCNKBr or KLiTaO (15); (47); (48); (49); (50); (55); (45); (46); (51); (54); (52); (53), where the two mixed components have similar sizes and shapes. Upon cooling below melting, they first form a cubic crystal without long-range orientational order in the plastic crystal phase. At lower , an orientational phase transition occurs, where the crystal structure becomes noncubic. In nondilute mixtures, this transition is diffuse with slow relaxations, where the orientations and the strains are strongly coupled, both exhibiting nanoscale heterogeneities(55); (56); (58); (57). Some polymers also undergo ferroelectric transitions due to alignment of permanent dipoles (13); (59); (60); (61). In particular, poly(vinylidene fluoride-trifluoroethylene) copolymers(62); (63) exhibited large electrostriction and relaxor-like polarization responses after electron irradiation (which brings disorder in polymer crystals). We also mention strain glass in shape-memory alloys (64), where the dipolar interaction does not come into play but a diffuse ferroelastic transition occurs with strain heterogeneities. We now recognize the universal features of glass coupled with a phase transition, where the order parameter fluctuations are frozen at low .

In their molecular dynamics simulation of relaxors, Burton et al.(65); (66); (67) started with a first-principles Hamiltonian for atomic displacements in perovskite-type crystals. As a compositional distribution, they assumed nanoscale chemically ordered regions embedded in a chemically disordered matrix. On the other hand, we investigate general aspects of ferroelectric glass with a simple molecular model. In electrostatics, we use an Ewald scheme including image dipoles and applied electric field (69); (68), which has been used to study water between electrodes(70); (71); (72). To prepare a mixed crystal, we cool a liquid mixture from high ; then, our impurity distribution at low is naturally formed during crystallization (58); (57).

Our system consists of spheroidal dipoles and spherical apolar particles only. Nevertheless, we can realize enhanced polarization fluctuations forming PNRs and calculate the frequency-dependent dielectric permittivity. We can also calculate the responses to applied electric field and to ZFCFC (zero-field-cooling and field-cooling) temperature changes. In the latter, nonergodicity of glass is demonstrated, so its experiments have been performed in spin glass(73); (75); (74), relaxors(25); (36); (76), orientational glass(48); (53), relaxor-like polymers (63) , and strain glass(64).

The organization of this paper is as follows. In Sec. II, we will explain our theoretical scheme and numerical method. In Sec. III, we will explain a structural phase transition in a one-component system of dipolar spheroids. In Sec. IV, we will examine diffuse ferroelectric transitions with impurities. Furthermore, we will examine responses to cyclic applied field in Sec.V and to cyclic temperature changes in Sec.VI. Additionally, antiferroelectric glass will be briefly discussed In Sec.VII.

## Ii Theoretical background

We treat mixed crystals composed of spheroidal polar particles as the first species and spherical apolar particles (called impurities) as the second species. These particles have no electric charge. As in Fig.1(a), we suppose smooth metal walls at and to apply electric field to the dipoles. The periodic boundary condition is imposed along the and axes with period . Thus, the particles are in a cell with volume .

In terms of the impurity concentration , the particle numbers of the two species are written as

(1) |

where the total particle number is set equal to . Their positions are written as (). The long axes of the spheroidal particles are denoted by unit vectors ().

### ii.1 Potential energy

The total potential energy is expressed as

(2) |

Here, is the sum of modified Lennard-Jones potentials between particles and (),

(3) |

where , is the characteristic interparticle energy, and in terms of the particle lengths and . The factor depends on the angles between spheroid directions and as (57); (58)

(4) |

where in the first term, in the second term, and represents the molecular anisotropy. For , we have , which vanishes for for (and for ). We assume a relatively small size difference and mild anisotropy as

(5) |

Then, at density , our system forms a crystal without phase separation and isotropic-nematic phase transition (77); (57). For larger and , the latter processes may take place during slow quenching from liquid. Because is minimized at for fixed and , we regard the anisotropic particles as spheroids with aspect ratio . Notice that our potential is similar to the Gay-Berne potential for rodlike molecules (78).

The second term in Eq.(2) is the sum of strongly repulsive, wall potentials as (69)

(6) |

We set and to make the potentials hardcore-like. Then, the distances between the dipole centers and the walls become longer than .

### ii.2 Electrostatic energy and canonical distribution

We assume permanent dipolar moments along the spheroid direction () written as

(7) |

where is a constant dipole moment. There is no induced dipole moment. The electric potential can be defined away from the dipole positions . We impose the metallic boundary condition at and :

(8) |

where is the applied potential difference and is the applied electric field. In this paper, we perform simulation by controlling (or ). In our scheme, can be nonstationary.

The boundary condition (8) is realized by the surface charge densities at and (see Appendix A). As a mathematical convenience, we instead introduce image dipoles outside the cell for each dipole at in the cell. As in Fig.1(a), we first consider those at () with the same moment , where is the unit vector along the axis. Second, at , we consider those with the image moment given by

(9) |

where is the image position closest to the bottom wall. For , the real and image dipoles and the applied field yield the following potential,

(10) | |||||

where , , and with , and being integers. Here, the first term is periodic in three dimensions (3D). Along the axis the period is because of the summation over or over the image dipoles. We confirm that the first term in Eq.(10) vanishes at and with the aid of Eq.(9).

At fixed , the total electrostatic energy in Eq.(2) is now written in terms of and as(68); (69)

(11) |

Here, is the dipolar tensor with its component being . In the first term, the self-interaction contributions and ) are removed in . In the second term, we set . In the last term, is the component of the total polarization,

(12) |

For each dipole , the electrostatic force is given by and the local electric field by

(13) |

We can also obtain by subtracting the self contribution from in Eq.(10) as

(14) |

We consider the Hamiltonian , where is the total kinetic energy. In our model, the applied field appears linearly in in Eq.(11). Then, we find

(15) |

where is the Hamiltonian for . This form was assumed in the original linear response theory(79). For stationary , the equilibrium average, denoted by , is over the canonical distribution . Then, for any variable (independent of ), its equilibrium average changes as a function of as(80)

(16) |

where is fixed in the derivative and . For the average polarization , we consider the differential susceptibility . In equilibrium, it is related to the variance of as

(17) |

As , tends to the susceptibility in the linear regime. In this paper, we calculate the time averages of the physical quantities using data from a single simulation run. In our case, the ergodicity holds at relatively high , but we do not obtain Eq.(17) at low because of freezing of mesoscopic PNRs in our finite system (see Sec.IVC and Fig.7).

### ii.3 Kinetic energy and equation of motions

The total kinetic energy depends on the translational velocities () and the angular velocities () as

(18) |

where is the mass common to the two species, and is the moment of inertia. We set in this paper. The Newton equations for are given by

(19) |

where . On the other hand, the Newton equations for () are of the form(69); (57),

(20) |

where , is the unit tensor, and is the local orientating field on dipole . The left hand side of Eq.(20) is perpendicular to from . The right hand side vanishes if is parallel to . From Eqs.(19) and (20) the Hamiltonian changes as (without thermostats). Thus, is conserved for stationary .

At low , we have for most , where is nearly parallel to . From Eq.(2) we set

(21) |

where is the long-range dipolar part in Eq.(13) and is the short-range steric part from the orientation dependence of in Eq.(3). Some calculations give

(22) |

where main contributions arise from neighbors with . These neighbor impurities () yield local random pinning fields (see Fig.3(a)).

### ii.4 Simulation method

We integrated Eqs.(19) and (20) for . We used the 3D Ewald method on the basis of in Eq.(11) (70); (68); (69); (71); (72). To realize crystal, we slowly cooled the system from a liquid above the melting temperature () at density . In crystal, there is no translational diffusion. We attached Nosé-Hoover thermostats to the particles in the layer regions and . We fixed the cell volume at with mostly, but we slightly varied in time to obtain the field-induced strain in Sec.V.

In our system, the dielectric response strongly depends on the dipole moment in Eq.(7)(12), so we present our results for and 1.6 in units of . For example(11); (12), if K and , these values of are D and 2.10 D, respectively.

We measure space and time in units of and

(23) |

Units of , electric potential, and electric field are , , and , respectively. For K and , we have V, Vnm, and (elementary charge).

## Iii Ferroelectric transition for

We first examine a ferroelectric transition in crystal composed of dipolar spheroids with in Eq.(4) without impurities. See similar simulation by Johnson et al.(12) for the prolate case with aspect ratio 1.25.

It is convenient to introduce orientational order parameters defined for each dipole as

(24) |

where and . We sum over neighbor dipoles with , where is their number. Then, represents the local dipolar order and the local quadrupolar order(46); (45). These variables will be used also for ferroelectric glass with .

In Fig.2, we examine the transition by slowly lowering for , , and . In (a), we plot the averages and . Here, the transition is steep but gradual due to the finite-size effect imposed by the metal walls at and . In our system, the spheroidal particles form a fcc crystal in the plastic crystal phase (46); (45); (57) in the range . For lower , a polycrystal with eight rhombohedral variants appears, where the spheroid directions are along except those near the interfaces.

In the transition range , the system is composed of disordered and ordered regions with sharp interfaces. We give snapshots of relatively ordered regions with at (b) and (c) , where we pick up (b) 10 and (c) of the total dipoles. These patterns are stationary in our simulation time intervals. In (d), at , we give a snapshot of polycrystal state with eight variants.

The rhombohedral structure is characterized by the angles of its lozenge faces of a unit cell. At low , we find for but for . See Sec.VA for the reason of this dependence.

## Iv Ferroelectric transition for

### iv.1 Role of impurities

The impurities hinder the spheroid rotations and suppress long-range orientational order not affecting the crystal order. In our previous papers (58); (57), this gave rise to orientational glass without electrostatic interactions. In a mixture of nematogenic molecules and large spherical particles, surface anchoring of the former around the latter suppresses the long-range nematic order (83).

In Fig.3(a), we display the dipole directions for c=0.2 on a (111) plane at . Many of them tend to align in the directions parallel to the impurity surfaces or perpendicular to (), because for in Eq.(3). However, this anchoring is possible only partially, because the dipoles are on the lattice points and the impurities form clusters. This picture resembles those of PNRs on crystal surfaces of relaxors (38); (39); (40).

In Fig.3(b), we display all the impurities in the cell for , where clustering is significant. As guides of eye, we write bonds between pairs of impurities if their distance is smaller than 1.4. In this bond criterion, we find large clusters composed of many members including a big one percolating through the cell. These clusters were pinned during crystallization, so they depend on the potentials and the cooling rate. They strongly influence the shapes of PNRs (see Figs.9 and 10 also).

Correlated quenched disorder should also be relevant in real systems. For relaxors, much effort(33); (34); (35) has been made to determine the distribution of the B-site ions (Mg and Nb for PMN) using effective atom-atom interactions, while Burton et al. (65); (66); (67) demonstrated strong influence of compositional heterogeneity on the PNRs.

### iv.2 Diffuse transition toward ferroelectric glass

The dipole moment determines relative importance of the dipolar and steric parts, and , in the orientating field in Eq.(21), since they are proportional to and , respectively, for . For example, for and , we average (, ) over all the dipoles to obtain (, ) for and and (, ) for and . Thus, is more important for larger in the dipole orientations. Here, the amplitude of the local electric field is mostly of order , where . This large size of is realized within mesoscopic PNRs.

In Fig.4, we examine the transition with and for the two cases and , where the net polarization nearly vanishes. At each , we waited for a time . In (a) and (b) we show gradual dependence of . They take appreciable values in the presence of small PNRs. In (c) and (d), we also show their variances,

(25) |

The orientation fluctuations are frozen at large sizes at low . We also see that in (a) and in (c) exhibit small maxima at low , but they should disappear in the ensemble averages.

In Fig.5, we plot the time-correlation functions for one-body angle changes defined by

(26) |

where the average is taken over the initial time . In (a) and (b), the angle changes slow down with lowering . We define the reorientation time by

(27) |

where 0.1 is smaller than the usual choice since decays considerably in the initial thermal stage for not very low . The PNRs are broken on this timescale. In (c) and (d), we display vs . where can well be fitted to the Vogel-Fulcher form (14),

(28) |

Here, , , and are constants with being for and for .

### iv.3 Dielectric permittivity

We next examine the dielectric permittivity. We calculated its real part and imaginary part as functions of and the frequency by applying small ac field in the linear response regime (see Appendix B).

In Fig.6, we show and the ratio vs at several low frequencies for the two cases (left) and 1.6 (right). In (a) and (b), increases with decreasing and exhibit a broad maximum at a temperature for each . With decreasing , decreases (with weaker dependence for smaller ) and the peak height increases. For , we have , so tends to the linear dielectric constant . However, for , decreases to zero with lowering or increasing , where the response of the PNRs to small ac field decreases. On the other hand, exhibits a maximum for each and shifts to a lower temperature with lowering . These behaviors characterize ferroelectric glass (17); (27); (16); (19); (20); (62); (23); (21); (22); (18); (50); (49). Similar behaviors were found for the frequency-dependent magnetic susceptibilities in spin glass(73). Furthermore, in Appendix B, we will present analysis of and for at relatively high on the basis of the linear response theory(79).

We write the inverse relation of as

(29) |

leading to . Here, () holds for (). In (c) and (d) of Fig.5, we compare the inverse and in Eq.(27) for and 1.6. We find . Thus, represents a characteristic frequency of the dipole reorientations. Previously, for relaxors and spin glasses, Stringer et al.(28) nicely fitted to the Vogel-Fulcher form, which is in accord with (c) and (d) of Fig.5.

Our system is ergodic at relatively high , but becomes nonergodic as is lowered. The boundary between these two regimes weakly depends on the observation time. In Fig.7(a), evolves on a wide range of time scales in a time interval with width for , , and . At , its time average becomes small, but its fluctuations are large. In contrast, at , it remains negative around , on which smaller thermal fluctuations with faster time scales are superimposed. Note that the ensemble average of should vanish at any for .

For a single simulation run, we consider the time average of the normalized polarization variance, written as . To avoid confusion, we define it explicitly as

(30) |

We set with here) for any time-dependent variable . This averaging procedure has already been taken for the quantities in Figs.4 and 5. In the nonergodic range, remains nonvanishing even for , while arises from the (thermal) dynamical fluctuations and tends to zero as . In Fig.7(b), we plot numerical results of for and at as functions of . These two curves nearly coincide for yielding , but is considerably larger than for . In their simulation, Burton et al.(65); (66) calculated a dielectric constant from polarization fluctuations, which corresponds to in our case.

In Fig.7(b), and steeply grow as . From the curves of and in its inset, and can fairly be fitted to the Curie-Weiss form,

(31) |

with and at ) for . At , however, we find and . In experiments, the behavior (31) was found for orientational glass(45), but a marked deviation was detected close to for relaxors(18); (21); (27). Thus, if is somewhat above , our polarization fluctuations resemble the critical fluctuations in systems near their critical point (23); (18). In our disordered system, these near-critical fluctuations are slowed down and eventually frozen as is further lowered, as in relaxors. This can also be seen in (c) and (d) of Fig.4. Furthermore, for , there is a tendency of interface formation between adjacent PNRs for , which will be discussed in future.

For relaxors, Stock et al.(32) divided the scattering intensity into frozen and dynamic parts, where the former (latter) increases (decreases) with lowering . Similar arguments of nonergodicity were made for polymer gels(84); (85), where the fluctuations of the polymer density consist of frozen and dynamic parts. Moreover, if gelation takes place in a polymer solution close to its criticality, the critical concentration fluctuations are pinned at the network formation (85); (80).

We next confirm Eq.(17) by increasing at with and , where the observation time is much longer than . In Fig.7(c), we compare the differential formula and the fluctuation formula for the field-dependent dielectric constant. The former is calculated from the data in Fig.12(a) and the latter from Eq.(30), where these two curves are surely very close. At this , the polarization fluctuations are suppressed with increasing .

### iv.4 Polar nanoregions in diffuse transition

In our diffuse transition, the PNRs are relatively ordered regions consisting of aligned clusters enclosed by impurities. At relatively high , they have finite lifetimes (within observation times)(21); (13). This feature is illustrated in two snapshots in Fig.8, which were taken at two times separated by in the same simulation run. They display the dipoles with for , , and . These two patterns are very different, so their lifetime is shorter than . In fact, is of order at in Fig.5(c).

The PNRs are frozen with lowering . In Fig.9, we give examples for , 0.2, and 0.3 with , , and . The left panels display the particles on the boundaries (), while the right ones the relatively ordered dipoles with . The dipoles depicted in the latter amount to 37, 20, and 13 of the total dipoles from above. For , we can see well-defined ordered domains consisting of eight variants, whose interfaces are trapped at impurities(57); (58) (see Fig.10(a)). These domains are broken up into smaller PNRs with increasing . For , the PNRs mostly take compressed, plate-like shapes under the constraint of the spatially correlated impurities (see Fig.10 also). For , the dipole orientations are highly frustrated on the particle scale without well-defined interfaces.

To be quantitative, we define PNRs as follows. In each PNR, any member satisfies and for some within the same PNR. In Fig.9, the dipole number in a PNR is , and on the average from above. Thus, the connectivity of the PNRs sensitively depends on . In the following, we treat the case .

### iv.5 Single polar nanoregion and local electric field

We visualize individual PNRs frozen at low . When the system is composed of PNRs, the local electric field in Eq.(13) arises mainly from the dipoles within the same PNR in the bulk. Its amplitude is of order for not very large .

In Fig.10, we pick up (a) a single PNR for at the cell center and (b) another one for () near the upper wall, where and . We depict the impurities whose distance to some dipole in the PNR is shorter than 1.4. We find the numbers of the constituent dipoles and impurities as (a) and (b) using the definition of PNRs in Sec.IVD. Here, the dipoles tend to be parallel to the impurity surfaces, as discussed in Sec.IVA, and almost all the impurities are on the PNR boundaries, resulting in plate-like PNRs.

In Fig.10 (right), we choose a typical dipole in the PNR interior (not in contact with the impurities) and display its and , where they are nearly parallel. Here, we divide the dipolar part of into the contributions from those inside and outside the PNR, written as and . Then, Eq.(A9) in Appendix A gives

(32) |

where the last term arises from the surface charges. In (a), we find , which occurs mostly for the dipoles in the interior of PNRs in the bulk. In (b), on the other hand, we find , where is of the same order as and is much larger than . Here,